We have derived mathematical formulas that describe the dynamics of malaria infections as random variables in cohorts of humans as they age (Henry JM, et al., 2024)1 This R package is the computational companion.

In the following, we review the mathematical formulas and the functions in ramp.falciparum that compute the multiplicity of infection (MoI), the age of infection (AoI), and the age of the youngest infection (AoY).

Formulas and Computation

  1. Let zτ(α,a)z_\tau(\alpha,a) denote the density for infections of age α\alpha in a host cohort of age aa born on day τ\tau.

    • The dynamics of z(α,a)z(\alpha,a) are described by the following: za+zα=rz \frac{\partial z}{\partial a} + \frac{\partial z}{\partial \alpha} = - r z with the boundary condition: zτ(0,a)=hτ(a). \begin{equation} z_\tau(0,a) = h_\tau(a). \end{equation} Its solutions are given by: zτ(α,a)=hτ(aα)erα z_\tau(\alpha,a) = h_\tau(a-\alpha) e^{-r \alpha}

    • The function zda computes zτ(α,a).z_\tau(\alpha, a).

  2. The mean MoI is given by the formula: mτ(a)=0azτ(α,a)dα m_\tau(a) = \int_0^a z_\tau(\alpha, a) d \alpha

    • The distribution of the MoI is Poisson (Nåsell I, 1985)2 with mean mτ(a).m_\tau(a).

    • The function meanMoI computes mτ(a)m_\tau(a) using zda

  3. The true prevalence is the pτ(a)=1emτ(a) p_\tau(a) = 1 - e^{-m_\tau(a)} The function truePR computes the true prevalence using meanMoI

  4. The density function for the age of infection (AoI) is Aτ(a)fA(α,a,τ)=zτ(α,a)mτ(a) A_\tau(a) \sim f_A(\alpha, a, \tau) = \frac{z_\tau(\alpha,a)}{m_\tau(a)} and its moments are xn(a,τ|h)=0aαnfA(α,a,τ)dα x_n(a, \tau | h) = \int_0^a \alpha^n f_A(\alpha, a, \tau) d \alpha

  5. The age of the youngest infection (AoY) is defined as: Yτ(a)fY(α,a,τ)=minζMτ(a){αi}i=1ζ,αiAτ(a) Y_\tau(a) \sim f_Y(\alpha, a, \tau) = \min_{\zeta \sim M_\tau(a)} \left\{ \alpha_i \right\}_{i=1}^\zeta, \alpha_i \sim A_\tau(a)

    • The density function can be expressed in terms of the density and distribution functions of the AoI and MoI. fY(α;a,τ)=fA(α,a,τ)emτ(a)FA(α,a,τ)mτ(a)pτ(a). f_Y(\alpha; a, \tau) = f_A(\alpha, a,\tau) e^{-m_\tau(a) F_A(\alpha, a,\tau)} \frac{m_\tau(a)}{p_\tau(a)}.

    • The distribution function for the AoY is: FY(a)1emτ(a)FA(α,a,τ)1emτ(a)=1emτ(a)FA(α,a,τ)pτ(a) F_Y(a) \sim \frac{1-e^{-m_\tau(a)F_A(\alpha, a,\tau)}}{1-e^{-m_\tau(a)}} = \frac{1-e^{-m_\tau(a)F_A(\alpha, a,\tau)}}{p_\tau(a)} \label{FY}

    • Its moments are: yn(a,τ|h)=0aαnfY(α|a,τ,h)dα y_n(a, \tau | h) = \int_0^a \alpha^n f_Y(\alpha | a, \tau, h) d \alpha

  6. We also developed functions to compute the age of the youngest of NN infections, called YoN Nτ(a)minN{αi}i=1N where αiAτ(a) N_\tau(a) \sim \min_{N} \left\{ \alpha_i \right\}_{i=1}^N \mbox { where } \alpha_i \sim A_\tau(a)

    • The distribution function for YoN, Nτ(a)N_\tau(a), is FN(α,a,t)1(1FA(α,a,τ))NF_N(\alpha, a, t) \sim 1- (1-F_A(\alpha, a, \tau))^N The following is a summary table of functions to compute the MoI, AoI, AoY, and all their moments.

    • The density function for YoN is found by differentiating:

fN(α,a,t)N(1FA(α,a,τ))N1fA(α,a,τ)mτ(a)f_N(\alpha, a, t) \sim N (1-F_A(\alpha, a, \tau))^{N-1}\frac{f_A(\alpha, a, \tau)}{m_\tau(a)}

Quick Reference

The following is a summary table of functions to compute the MoI, AoI, AoY, YoN, and all their moments.

MoI AoI AoY YoN
ζ\zeta α\alpha α\alpha α\alpha
ζ0\zeta \geq 0 0αa0 \leq \alpha \leq a 0αa0 \leq \alpha \leq a 0αa0 \leq \alpha \leq a
Random Variable Mτ(ζ,a,h)M_\tau(\zeta, a, h) Aτ(α,a,h)A_\tau(\alpha, a, h) Yτ(α,a,h)Y_\tau(\alpha, a , h) Nτ(α,a,h)N_\tau(\alpha, a, h)
Density Function fM(ζ,a,h)f_M(\zeta, a, h) fA(ζ,a,h)f_A(\zeta, a, h) fY(ζ,a,h)f_Y(\zeta, a, h) fN(ζ,a,h)f_N(\zeta, a, h)
dpois dAoI dAoY dYoN
Distribution Function FM(ζ,a,h)F_M(\zeta, a, h) FA(ζ,a,h)F_A(\zeta, a, h) FY(ζ,a,h)F_Y(\zeta, a, h) FN(ζ,a,h)F_N(\zeta, a, h)
ppois pAoI pAoY pYoN
Random Numbers M̂τ(ζ,a,h)\hat M_\tau(\zeta, a, h) Âτ(α,a,h)\hat A_\tau(\alpha, a, h) Ŷτ(α,a,h)\hat Y_\tau(\alpha, a , h) N̂τ(α,a,h)\hat N_\tau(\alpha, a , h)
rpois rAoI rAoY rYoN
Moments mτ(a,h)m_\tau(a, h) xn(a,τ,h)x_n(a, \tau, h) yn(a,τ,h)y_n(a, \tau, h)
meanMoI momentAoI momentAoY

Demonstration

Force of Infection (FoI)

To compute anything, we must first set up a function to describe exposure (see the FoI vignette). We define functions that plot the FoI for a cohort as it ages (in red), but we can also compute the population average FoI (in black). Different cohorts would experience different histories of exposure.

Computing zda

The function ramp.falciparum::zda(alpha, a, FoIpar, ...) uses the formula in Eq. 1 to compute the density of parasite infections in a cohort of humans as it ages.

Using zda, we can compute the density of parasites in a cohort of any age without solving a full system of equations. Given a function describing the FoI in the population, h(t)h(t), and the cohort birthday, τ.\tau.

alpha = 60
a = 6*365
zda(60, 6*365, foiP3) 
## [1] 0.001601196

The following computes the density of infections of every age in a cohort of age 3.

zz = zda(a3years, max(a3years), foiP3)

When we plot zτ(α,a)z_\tau(\alpha, a), we note that as α\alpha grows larger, the parasite cohort gets older. When we plot parasite cohorts by age, time is going backwards on the x-axis.

Now, we can imagine what zda would look like for several different host cohorts at age three, but who were born at different times. In effect, we are taking a snapshot of the cohorts at the same age, but at different times.

The curves are different because the hosts were born at different months, and they thus experienced different levels of exposure over the first two years of life. Here the annual FoI is 5 infections, per person, per year (h=5/365\bar h = 5/365):

Multiplicity of Infection (MoI)

We define a random variable MM describing the multiplicity of infection (MoI). The distribution of the MoI is Poisson (see the MoI vignette).

Mτ(a)fM(ζ;a,τ)=Pois(mτ(a))M_\tau(a) \sim f_M(\zeta; a, \tau) = \mbox{Pois}(m_\tau(a))

Since zτ(α,a)z_\tau(\alpha, a) describes the density of all infections of age α\alpha in a cohort of age aa, the density of all infections must be the MoI. Since 0α<a0 \leq \alpha < a, it must be true that:

mτ(a)=0azτ(α,a)dα\begin{equation} \tag{2} m_\tau(a) = \int_0^a z_\tau(\alpha, a) d \alpha \end{equation}

The function that computes mτ(a)m_\tau(a) is called meanMoI.

mm = meanMoI(a3years, foiP3, hhat=5/365)

Here, we plot the average MoI in the host cohort as it ages:

Age of Infection (AoI)

We define a random variable Aτ(a)A_\tau(a) that describes the age of infection (AoI), which is given by the formula

Aτ(a)fA(α;a,τ)=zτ(α,a)mτ(a) A_\tau(a) \sim f_A(\alpha; a, \tau) = \frac{z_\tau(\alpha,a)}{m_\tau(a)}

The Density Function, dAoI

We can compute Aτ(a)A_\tau(a) using the density function dAoI:

f_A = dAoI(a3years, max(a3years), foiP3)

Now, as we plot the distribution of the AoI in cohorts at age two, born at different months (as we did above), we notice that the distributions have changed shapes:

The Distribution Function, pAoI

The distribution function for Aτ(a)A_\tau(a) is:

FA(a)0αfA(α;a,τ)dα F_A(a) \sim \int_0^\alpha f_A(\alpha; a, \tau) d\alpha

F_A = pAoI(a3years, max(a3years), foiP3)

If our functions work correctly, then we should get approximately the same answer from computing the cumulative sum of dAoI.

F_A_alt = cumsum(f_A)

We shouldn’t expect the answers to be exactly the same, but they should be close, with the pAoI in black.

par(mar = c(5,4,1,1))

plot(a3years, F_A, type = "l", 
     xlab = "Parasite Cohort Age", 
     ylab = expression(1-F[X](alpha, a, tau)), lwd=3)

lines(a3years, F_A_alt, col = "red", lwd=2, lty =2)

Random Numbers, rAoI

The function rAoI uses pAoI to generate random numbers from FA(α)F_A(\alpha)

rhx = rAoI(10000, 3*365, foiP3)

A simple visual check computes the empirical CDF for the random variates against FA(α)F_A(\alpha) computed using pAoI

par(mar = c(5,4,1,2))
plot(stats::ecdf(rhx), xlim = c(0,1095), cex=0.2, main = "", 
     xlab = expression(list(alpha, paste("Parasite Age (in Days)"))), 
     ylab = expression(list(F[A](alpha), paste("ecdf"))))
lines(a3years, F_A, col = "red", lty = 2, lwd=2)

We can also plot the distribution functions.

AoI Moments

Let xx denote the first moment of of Aτ(a)A_\tau(a): xτ(a)=Aτ(a)=0αzτ(α,a)mτ(a)x_\tau(a) = \left< A_\tau(a) \right> = \int_0^\infty \alpha \frac{z_\tau(\alpha, a)} {m_\tau(a)}

Similarly, we let xτ(a)[n]x_\tau(a)[n] denote the higher order moments of Aτ(a)A_\tau(a): x[n](a,τ)=0αnzτ(α,a)mτ(a)x_{[n]}(a, \tau) = \int_0^\infty \alpha^n \frac{z_\tau(\alpha, a)} {m_\tau(a)}

moment1 = momentAoI(a3years, foiP3)
moment2 = momentAoI(a3years, foiP3, n=2)
moment3 = momentAoI(a3years, foiP3, n=3)

The first three moments of the AoY plotted over time. In the top plot, we’ve also plotted the nthn^{th} root of the nthn^{th} moment.

Age of the Youngest Infection (AoY)

We have derived a random variable Yτ(a)Y_\tau(a) describing the age of the youngest infection (AoY). The density function for the AoY is:

Yτ(a)fY(α;a,τ)=fA(α,a,τ)emτ(a)FA(α,a,τ)mτ(a)pτ(a)Y_\tau(a) \sim f_Y(\alpha; a, \tau) = f_A(\alpha, a,\tau) e^{-m_\tau(a) F_A(\alpha, a,\tau)} \frac{m_\tau(a)}{p_\tau(a)} The distribution function is:

FY(a)1emτ(a)FA(α,a,τ)1emτ(a)=1emτ(a)FA(α,a,τ)pτ(a)F_Y(a) \sim \frac{1-e^{-m_\tau(a)F_A(\alpha, a,\tau)}}{1-e^{-m_\tau(a)}} = \frac{1-e^{-m_\tau(a)F_A(\alpha, a,\tau)}}{p_\tau(a)}

The derivations are found in a Suppplement to Henry JM, et al. (2024).

The mean AoY is:

Yτ(a)=0aαfY(α,a,τ)dα \left< Y_\tau(a) \right> = \int_0^a \alpha \; f_Y(\alpha, a, \tau) \; d\alpha

And the higher order moments for the AoY are:

Yτ(a)n=0nαnfy(α,a,τ)dα\left< Y_\tau(a)^n \right> = \int_0^n \alpha^n \; f_y(\alpha, a, \tau) \; d\alpha

AoY Density, dAoY

The density function is computed with the function dAoY.

f_Y = dAoY(a3years, 3*365, foiP3)

We can compare fY(α)f_Y(\alpha) (in black) to fA(α)f_A(\alpha) (in grey).

Random Variables, rAoY

raoy = rAoY(10^5, 3*365, foiP3)
hist(raoy, breaks=seq(0, 1095, by = 15), 
     right=F, probability=T, main = "", 
     xlab = expression(list(alpha, paste("Parasite Age (in Days)"))), 
     border = grey(0.5)) -> out
lines(a3years, f_Y, type = "l", col = "red") 

AoY Moments

aa = seq(5, 3*365, by = 5) 
moment1y = momentAoY(aa, foiP3)
moment2y = momentAoY(aa, foiP3, n=2)
moment3y = momentAoY(aa, foiP3, n=3)

The first three moments of the AoY plotted over time. In the top plot, we’ve also plotted the nthn^{th} root of the nthn^{th} moment.


Next:

  • In the vignette MoI, we show that the mean MoI computed using the formula in Eq. 2 gives the same answer as other approaches.

  1. Henry JM, Carter AR, Wu SL, Smith DL (in preparation). A Probabilistic Synthesis of Malaria Epidemiology: Exposure, Infection, Parasite Densities, and Detection.↩︎

  2. Nåsell I (1985). Hybrid Models of Tropical Infections, 1st edition. Springer-Verlag. https://doi.org/10.1007/978-3-662-01609-1↩︎