hybrid = solve_dm(5/365, foiP3, Amax=1095)The Superinfection Problem
A hybrid variables approach
What is the waiting time to clearing malaria in a model with superinfection? We present a closed-form solution to a classical mathematical problem in malaria epidemiology.
The Problem
Immunity to malaria does not necessarily block new infections, infections can take a long time to clear, and exposure to malaria can be quite intense. Humans can thus be reinfected with malaria multiple times, which was called superinfection. This gave rise to a question about how to model clearance rates with superinfection.
History
In 1950, George Maccdonald published a model of malaria superinfection [1]. He was trying to understand how reexposure and reinfection would affect the prevalence of infection in a population. In particular, if people are continuously becoming reinfected, then how long would it take until a person would clear all of the parasites to become uninfected again?
Let \(h\) denote the force of infection, and \(r\) the rate that a simple infection would clear. Let \(x\) denote the fraction of a population that is infected. Ross’s model was equivalent to: \[
\frac{dx}{dt} = h(1-x) - r x
\] Macdonald’s model [1] was
\[
\frac{dx}{dt} = h(1-x) - F(h,r) x
\] where \[
F(h,r) =
\begin{cases}
h-r & \mbox {if } h<r \\
0 & \mbox {if } h \geq r \\
\end{cases}
\] In the text of the paper, Macdonald had described a model where infections would clear independtly of one another. The model did not match the description of the process. In 1975, Paul Fine wrote an essay about it that we recommended highly [2].
In 1974, Dietz presented slightly different solution to the problem of clearance with superinfection in a mathematical model developed for the Garki Project [3]; the clearance rate was a function of the force of infection (\(h\)) and the clearance rate (\(r\)):
\[F(h,r) = \frac{h}{e^{h/r} -1}\]
Dietz’s function was based on a steady state solution to a model from queuing theory, called \(M/M/\infty\).
Queuing Theory
The model for parasite clearance described by Macdonald is called \(M/M/\infty\). In this model:
each new infection increases the MoI by one;
the transition rate to a higher MoI is the force of infection (FoI), denoted \(h\);
Each parasite can clear at some rate, \(r\).
The following diagram shows how the fraction of the population with a given MoI, denoted \(\zeta_i=i\), changes over time.
\[\begin{equation} \begin{array}{ccccccccc} \zeta_0 & {h\atop \longrightarrow} \atop {\longleftarrow \atop r} & \zeta_1 & {h\atop \longrightarrow} \atop {\longleftarrow \atop {2r}} & \zeta_2 & {h \atop \longrightarrow} \atop {\longleftarrow \atop {3r}} & \zeta_3 & {h \atop \longrightarrow} \atop {\longleftarrow \atop {4r}}& \ldots \end{array} \end{equation}\]
Changes in the fraction uninfected (MoI=0) are described by:
\[\frac{d \zeta_0}{dt} = r \zeta _1 - h \zeta_0.\]
For those who are infected (with MoI\(=i\)), the dynamics are:
\[\frac{d \zeta_i}{dt} = h \zeta_{i-1} + r(i+1) \zeta_{i+1} - (r+h) \zeta_i\] The software in ramp.falciparum includes functions to set up and solve the queuing model \(M/M/\infty.\)
While \(M/M/\infty\) is an infinite system of differential equations, the function solveMMinfinity solves a finite system of differential equations. The maximum MoI is chosen to be large enough that it doesn’t affect the results.
Hybrid Models
An elegant solution to the problem of superinfection was proposed by Nåsell, who developed a hybrid modeling approach [4].
While the infinite system of differential equations can be formulated and solved numerically, Nåsell showed that the system can be described by a simple equation. Let \(m\) denote the mean MoI.
\[ \frac{dm}{dt} = h-rm \]
He also showed that the distribution of the MoI would converge to a Poisson distribution. If the initial distribution was Poisson, then it would remain Poisson forever.
with(hybrid, plot(time, m, type = "l", ylab = expression(m[tau](a)), xlab = "a - cohort age (in days)"))Superinfection Dynamics
With all this information at hand, we are now in a position to formulate a model for superinfection.
We note that since \(m\) is the mean of a Poisson, the prevalence of infection can be computed directly from \(m\). It is the complement of the zero term from a Poisson:
\[x(t) = 1-\mbox{Pois}(\zeta =0; m(t)) = 1 - e^{-m(t)}\]
It is, nevertheless, sometimes useful to know the clearance rate, a transition from MoI = 1 to MoI=0. Clearance from the infected class can only occur if a person is infected with an MoI of one, and if that infection clears:
\[ R(m) = r \frac{\mbox{Pr}(\zeta=1; m(t))}{\mbox{Pr}(\zeta>0; m(t))} = r \frac{m e^{-m}}{1-e^{-m}} = r \frac{m}{e^m-1} \]
The change in prevalence is thus described by the equation:
\[ \frac{dx}{dt} = h(1-x)-R(m) x \]
While this is interesting, and while it solves the problem of modeling clearance under superinfection, the equation depends on \(m\). Since we can compute \(x(t)\) directly from \(m(t)\): \[m = -\ln (1-x)\] So we can express the per-capita clearance rate as a function of \(x.\)
\[ \frac{dx}{dt} = h(1-x) - F_r(x) x \] where
\[ F_r (x) = \begin{cases} -r \ln (1-x) \frac{1-x}{x} & \mbox {if } x > 0 \\ r & \mbox {if } x = 0 \\ \end{cases} \]
So if \(x>0\), then can rewrite the equation:
\[ \frac{dx}{dt} = h(1-x) + r (1-x) \ln (1-x) \]