A Hybrid Model for the Effects of Treatment and Demography on Malaria Superinfection

As standard mathematical models for the transmission of vector-borne pathogens with weak or no apparent sterilizing immunity, Susceptible-Infected-Susceptible (SIS) systems such as the Ross-Macdonald equations are a useful starting point for modeling the impacts of interventions on prevalence for diseases that cannot superinfect their hosts. In particular, they are parameterizable from quantities we can estimate such as the force of infection (FOI), the rate of natural recovery from a single infection, the treatment rate, and the rate of demographic turnover. However, malaria para-sites can superinfect their host which has the effect of increasing the duration of infection before total recovery. Queueing theory has been applied to capture this behavior, but a problem with current queueing models is the exclusion of factors such as demographic turnover and treatment. These factors in particular can affect the entire shape of the distribution of MOI generated by the super-infection process, its transient dynamics, and the population mean recovery rate. Here we show the distribution of MOI can be described by an alternative hyper-Poisson distribution. We then couple our resulting equations to a simple vector transmission model, extending previous Ross-Macdonald theory.


22
Malaria continues to have a significant impact on the world, and is responsible for over 400,000 23 deaths in 2015 alone [4]. Over the past century a repertoire of interventions has been developed 24 to prevent its spread and reduce its impact on those afflicted, but the parasite continues to adapt 25 and persist. Because malaria dynamics and control are so complex, mathematical and statistical 26 models play a significant role in the strategic implementation of these existing tools to maximize 27 their impact. One of the historically difficult aspects of malaria and other macroparasitic diseases to 28 model is the process of superinfection, wherein multiple distinct cohorts of pathogens infect the host 29 simultaneously. The standard SIS system usually used to model similar disease transmission using 30 a constant recovery rate from a single infected class does not capture this type of heterogeneity, and 31 may therefore be a poor approximation to the real system. 32 There has been a long history of modeling a number of disease states to track the proportion 33 of individuals with n infections, each of which is still susceptible to further infection [2,3,8]. Care 34 must be taken in modeling the distribution of the multiplicity of infection (MOI), that is the fraction 35 of the population that has exactly n "broods" or cohorts of parasites. This distribution strongly 36 affects the population mean recovery rate, which can be seen as a weighted sum of the rate at 37 which individuals are recovering from n infections weighted by the proportion of individuals with n 38 infections. The shape of this distribution can be profoundly altered by factors such as seasonality 39 and case management, which in turn can have considerable effects on the mean rate of recovery for 40 the population and the prevalence of the disease. 41 Under George Macdonald's original assumption that infection status is no barrier to reinfection, 42 MOI can be modeled as an M/M/∞ queue [2,8]. Though Macdonald's written description of the 43 system was correct, he mistakenly presented it mathematically as an M/M/1 queue [3]. This system 44 vastly overestimates an individual's time spent infected, with the average becoming infinite if the 45 force of infection exceeds the rate of recovery. Macdonald's model and its description generated 46 some confusion, substantial discussion, and a series of approximating models [10,12]. 47 However all of this analysis has historically been focused on the original M/M/∞ model, whose 48 assumptions include a strong steady state condition on the population wherein individuals experi-49 ence no demographic turnover or external forces that influence their susceptibility or infection status. 50 Critically, finite lifespans and treatment will affect the distribution of MOI. Including their inter-51 actions ad hoc in ordinary differential equation (ODE) models may lead to inconsistent transition 52 rates, and even inconsistent equilibria, when compared to the original assumptions of the stochastic 53 model.

54
For models in which infections are assumed to clear independently, the mean population recovery 55 rate is dependent on the proportion of individuals with a single infection -therefore the distribution 56 of MOI strongly influences the recovery rate. Excluding demographics and allowing for variable 57 force of infection (FOI, here denoted λ) allows hybridization of the transmission system given the 58 distribution is Poisson, and it was modeled elegantly by Nåsell with two relatively simple ODEs 59 [9]. However, introducing demography and treatment requires a different approach as they can 60 fundamentally alter the process generating the distribution of MOI.

61
Here we consider the system with demography and treatment, with and without transient chemo-62 protection, and present a new way of modeling the statics and dynamics of the distribution of MOI 63 and its nonlinear effects on the population mean recovery rate and resulting prevalence. In addition, 64 we show that this distribution is stable in the sense that any initial distribution with finite mean 65 converges in distribution to the stationary distribution. We then show the impact of an extended 66 protective period post-treatment, which has the effect of inducing a zero-inflation in the MOI dis-67 tribution. Using these results and their resulting dynamics, we then extend the Ross-Macdonald 68 model to include this heterogeneity and explore the effects of mortality and treatment on prevalence 69 under these assumptions. Consider a Poisson process X t (λ) of new infections in an individual person with intensity λ, the 74 force of infection (FOI). We will assume λ may be time-varying, but it is independent of the current 75 number of infections as per Macdonald's original assumption. In addition, we will have a second 76 Poisson process Y t (rX t ) of recovery with an intensity rX t . Again, we make no assumption about 77 r aside from it being independent of X t or Y t . Then the M/M/∞ queue can be represented by a 78 third Poisson process, defined to be the difference of the two: Z t := X t − Y t . This process counts 79 the current number of infections in an individual, and its range is on the non-negative integers. To 80 add demography and treatment, which both act to decrease from the infected classes and increase 81 the susceptible class, we need additional Poisson processes D t (δ) and T t (τ ), where δ and τ are 82 respectively the mean birth/death rate (assumed equal to maintain constant population sizes) and 83 the mean treatment rate. Treatment only acts on those with at least one infection, and because births 84 are assumed to replace deaths, demography also has a net effect of removing infected individuals 85 and introducing susceptible individuals. Therefore for brevity we can replace the two processes with 86 a single Poisson process M t with rate µ, where µ := δ + τ is the expected rate at which either death 87 or treatment occurs.

88
The master equations for the system can be built by considering the state space in Figure 1: dx n dt = λx n−1 + r(n + 1)x n+1 − (λ + rn + τ + δ)x n which we can rewrite using the combined treatment and demographic process with parameter µ: As this is an upper triangular infinite system of linear ODEs, it can't be solved iteratively; instead 91 we will introduce a generating function of the form x n s n which encodes the infinite variables as coefficients in its power series, simplifying further analysis. 93 We can directly substitute: [−(µ + λ + rn)x n + λx n−1 + r(n + 1)x n+1 ]s n This simplifies to the following partial differential equation (PDE) for the probability generating 95 function (PGF) G: Because G is a probability generating function, a natural boundary condition is imposed wherein 97 lim s→1 − G(t, s) = 1 owing to the fact that the sum of the coefficients of the power series is equal to 98 the sum of all the probabilities. 99

100
Note that if we let µ = 0 then we obtain the equation whose solution is the generating function 101 of a Poisson distribution, aligning with known results [2]. Including demography or treatment forces 102 the stationary distribution to belong to a different class. To find this class, we will start by solving 103 for the stationary distribution whose generating function can be found by setting ∂G ∂t = 0, leading 104 to the following ODE: with the condition that lim s→1 − G = 1 which is true for any probability generating function. 107 In the static case, the stationary distribution is parameterized by two parameters: α = λ/r cor-108 responding to the expected number of infections per recovery, and β = µ/r which is the expected 109 number of demographic events or treatments per recovery. As the parameters are independent of s, 110 this equation can be solved with the application of an integrating factor: where γ(a, x) is the lower incomplete gamma function. Because any PGF has the property that 112 the left limit as s approaches 1 should equal 1 and the term next to the integration constant has a 113 singularity at 1, the constant must be zero. Therefore the PGF for the steady state distribution of 114 MOI is given by Replacing γ with its power series representation yields where y (n) is the rising factorial of y. This is the PGF of the alternative hyper-Poisson distribution, 118 AHP (α, β + 1) [7], whose probability mass function (PMF) is This distribution is a more general form of a Poisson distribution with β acting as a dispersal 120 parameter. As β > 0 for any model that includes demography or treatment, the distribution is 121 overdispersed compared to a Poisson distribution. A plot of the distribution for varying β is given 122 in Figure 2. Note that despite qualitative similarities, it is distinct from the negative binomial 123 distribution in both form and derivation -and in many cases due to its similarity to the Poisson 124 Figure 2: Plots of the PMF of the stationary distribution of MOI with varying values of treatment and demography, corresponding to the alternative hyper-Poisson (AHP) distribution. α, which is proportional to the force of infection, is held constant at 7. We vary values of β, defined here to be the ratio of the rate of treatment or death to the rate of recovery. β is varied from 0 to 2 in increments of .1. The bold line represents a Poisson distribution, β = 0. Increasing β suppresses MOI and increases the probability of having fewer infections, and it increasingly resembles a geometric distribution distribution, as it is in the same family, it is sometimes a more appropriate alternative overdispersed 125 distribution.

126
Knowing the distribution of MOI is sufficient to get an estimate of the expected prevalence. 127 Prevalence in this model is the complement of the proportion of individuals with zero infections, 128 which is equal to the PGF evaluated at zero. Figure 3 shows plot of the prevalence as a function of 129 the force of infection under varying treatment rates using this conversion. As expected, prevalence 130 increases with FOI and decreases with treatment. 131

132
Now that we know the stationary distribution, we can turn our attention to the dynamics of 133 MOI. In order to proceed, note that the transient distribution will always be AHP-distributed with α 134 and β as being (possibly time-varying) parameters for the distribution. As these parameters depend 135 on the three original parameters λ, r, and µ, they satisfy their own ODEs which will be derived here. 136 This invariance of the class of AHP distributions under the dynamics can be seen by noting the set of 137 PGFs corresponding to the class of AHP distributions constitute the only nontrivial solutions to the 138 Figure 3: Plot of the prevalence as a function of annual FOI, excluding short term protection from treatment. The curves represent the expected prevalence for increasing FOI with set treatment rates. Prevalence monotonically decreases with increasing treatment rate. From top to bottom, the curves represent treatment rates of 0, 1, 3, 6, 9, and 12 per year. PDE, and once the system has entered this class it will not leave it without external perturbation. 139 To begin we go back to our PDE for the PGF and take a partial derivative with respect to s: Plugging in s = 1 gives us the governing equation for the first factorial moment, also equal to the 141 first raw moment: we then repeat the previous process with a second partial derivative to find an equation for the 143 second factorial moment: ∂G ∂s from which we can derive the equation for the second raw moment using the relationship between the second raw moment and the second factorial moment, that is We then can match parameters to moments to determine the relationship between the transient 145 moments and the transient parameter values that codetermine the exact distribution.

146
For this we note that for a random variable X ∼ AHP (α, β + 1), Matching moments allows us to get a functional relationship between the dynamics of the moments 149 and the shape of the distribution. Doing so gives us solving the system in terms of α and β gives us In particular, if we assume the system reached the limiting PMF at some point in the past 153 without subsequent perturbation, then this combined with the PMF of the AHP distribution gives 154 us a method of modeling the prevalence with a much more manageable system of two ODEs for the 155 two moments necessary to parameterize the transient distribution: which we can convert back to prevalence x by noting that prevalence is the complement of the 157 probability of having zero infections: Note that we can also find x n (t), the proportion of individuals with n infections, by considering 159 the PMF of the stationary distribution: with the same definitions as before.

161
In addition, this means that if we are given an explicit time-dependent functional form for 162 the FOI λ, we would be able to find the closed form solution for prevalence over time by solving 163 the inhomogeneous system of two linear equations for the moments iteratively through the use of 164 integrating factors. That is, in the language of queueing theory, we could find the transient solution 165 for the PMF of states given the initial distribution is AHP-distributed. This assumption will be 166 relaxed and a more general argument for the behavior of solutions will be made in the next section. Now that we know the stationary distribution, we can check its stability in the set of distributions. 170 That is, if we don't start off initially with an AHP distribution or it is perturbed through a non-171 uniformly acting intervention such as a targeted treatment or vaccine, will the distribution of MOI 172 return to being AHP-distributed? If so, the dynamics above which assume AHP-distributed MOI 173 will be good long-term approximations for any initial distribution. To investigate this, consider the 174 PDE for the PGF: As this is a linear equation, we see we can arrange it to be in the following form: We will restrict the domain for L to the following: It's clear that this forms a complete subspace of C 1 [0, 1] under the L 2 norm.

178
Now we can define where G is the generating function for the transient distribution at time t and G * is the generating 180 function for the stationary distribution found before. In particular this means G * ∈ {G : ∂G ∂t = 0}. 181 As the left hand limit of each of the functions on the right at 1 is 1 because they are PGFs, the 182 left hand limit of h must be 0. Also, as long as the mean of G is finite, h will be once continuously 183 differentiable. Therefore, h ∈ D.

184
Now that we know h is in our domain, we can consider the action of ∂ ∂t : where the last equality holds because G * being stationary implies that LG * = −µ. Now we can 186 investigate the spectrum of L| D , that is the operator restricted to the subspace that h is in. Setting 187 up the eigenvalue problem, for any eigenfunction u ∈ D we get This implies uniform convergence of any initial PGF to that of an AHP distribution on the unit 202 interval. Therefore any initial distribution with finite mean converges in distribution to an AHP 203 distribution, so it is globally asymptotically stable in this domain as long as the parameters are 204 non-negative and asymptotically constant with r > 0 for all time. For this system µ here acts as a 205 sort of bifurcation parameter, where µ = 0 is the boundary between stability and instability. It is 206 still stable at this boundary because any function which is not in the null space of the linear operator 207 will have a strictly negative time derivative, which implies the asymptotic stability of the stationary 208 distribution holds for the standard M/M/∞ queue with its stationary Poisson distribution as well. Now we will consider the same system as before but we further assume that treatment confers 212 some temporary protection from subsequent infection -once someone is treated they enter a protected 213 state as the treatment will linger in the body protect against new infections for on average 1/ω units 214 of time before the individual returns to a susceptible state. The distribution of MOI conditioned 215 on not being protected will remain AHP-distributed because adding this compartment does not 216 change the queuing system -it adds a "waiting room" between bouts in the queue. Repeating the 217 analysis done before by introducing a PGF results in another AHP-distribution, but with a zero 218 inflation introduced by the proportion of individuals in the protected class. Treatment with waning 219 protection therefore has three major effects on the distribution of MOI: decreasing the expected 220 value, inducing overdispersion relative to a Poisson distribution, and zero-inflation. Each of these 221 effects can individually correspond with an overall decrease in prevalence. See figures 5 and 6 for 222 plots of the PMF of MOI and prevalence curves, respectively.
Note that x has a nice physical interpretation in this form -it is the probability that an individual has 227 at least 1 infection given that they are not protected, times the probability they are not protected. 228 Figure 6: Plot of the prevalence as a function of annual FOI, including short term protection from treatment. Average duration of protection is set to 30 days. The curves represent the expected prevalence for increasing FOI with set treatment rates. Prevalence monotonically decreases with increasing treatment, but at high treatment rates prevalence will drop with increasing FOI due to the increase to the protected class. From top to bottom, the curves represent treatment rates of 0, 1, 3, 6, 9, and 12 per year.

229
Now that we have transient behavior established, we can go one step further and include a 230 pseudoequilibrium analysis. That is, rather than having to deal with variables which parameterize 231 the distribution we can assume the distribution approaches equilibrium on timescales that are much 232 smaller than the rate at which the force of infection changes. This was implicitly made in the model 233 developed for the Garki project [2], wherein the rate of recovery was given as the reciprocal of the 234 expected duration of infection. If we denote the time spent infected as X, then the expression for 235 the expectation is as follows: 236 E[X] = e λ/r − 1 λ It's useful to compare this to the case without superinfection. Because each infection has an 237 exponentially-distributed time to clearance, without superinfection the expected time spent infected 238 is 1/r. In the case with superinfection, we can see first that it is monotonically increasing with 239 respect to the force of infection. In addition, for small FOI the Taylor expansion of the exponential 240 shows that the expected time is approximately 1/r: Now we can consider the same simplification with treatment and demography. In order to 242 consider this case, we need to know how this expectation changes when introducing exponentially-243 distributed treatment and demography. Although the distribution of the time to recovery starting 244 at a single infection is not known in closed form, its Laplace transform is -and we can exploit that 245 to determine the impact of these perturbations [5]. Denoting the time until treatment (or death) 246 with T , then the total duration of infection D is D = min{X, T }. That is, the infection will end if 247 an individual either recovers naturally or gets treated or dies, whichever comes first. Then we can 248 write down the following: where this step uses the independence of the natural recovery and treatment/demography processes. 250 As these are exactly the CDFs of the respective variables, we get Because treatment and death are exponentially distributed, the minimum of the two is also 252 exponentially distributed with a rate equal to the sum of their rates. We will again denote this new 253 rate as µ. Knowing this allows us to plug in the CDF of this exponential distribution for F T : Using the fact that the integral of the complement of the CDF is equal to the expectation for any 255 random variable, we get where here L denotes the Laplace transform. Finally, using the relationship between the Laplace 257 transform and the integral of a function gives us a formula for the expected duration of infection as 258 a function of the Laplace transform of X: Although this rate looks complicated it is consistent; the limit as µ goes to zero is exactly the rate 263 found previously and implemented in the Garki model [2].

264
Before plugging this in to a model of infection and recovery, we must also determine the proportion of people who were treated rather than died or recovered, as they will enter the protected class. This can by found again through a probabilistic argument: That is, the probability of being treated to end the infection is the product of the probabilities 265 of ending the treatment through either treating or dying given that an individual won't recover 266 naturally, times the probability that the individual won't recover naturally. The first probability is 267 trivial; it's the probability that one exponential random variable is less than another, and in this 268 case is τ /(τ + δ) = τ /µ. The second probability is more interesting, since again with superinfection 269 X is not exponentially-distributed. Again we can write this probability out: which shows that the Laplace transform of a random variable X can also be interpreted as the 271 probability that a random variable defined over the nonnegative reals is less than or equal to an 272 exponential random variable with rate µ. Note that this quantity as a probability is bounded between 273 0 and 1, and as µ increases to infinity this probability decreases to 0. Multiplying those probabilities, 274 we see Therefore, we can complete the simplified model by replacing the simple rate of recovery r from 276 infected to susceptible with the reciprocal of this duration of infection as follows: x − ωp which we found can be written as is the probability an infectious human will infect a susceptible mosquito during a blood meal, g is 306 the death rate of mosquitoes, and ν is the extrinsic incubation period. These ODEs are derived in 307 the same manner as before, this time tracking the extra compartment.  It's important to recognize that this will hold if the initial distribution is in the null space 316 of the linear operator described in the previous section, which corresponds to the class of AHP 317 distributions. However, because we proved any distribution asymptotically approaches the class of 318 AHP distributions, knowing malaria has been endemic in a region for a long period of time without 319 strong perturbation to the structure of the system is enough to make this a reasonable assumption 320 for modeling local transmission.

321
Finally, in the specific case where there is no lasting protection from treatment we can also couple 322 the pseudoequilibrium simplification to transmission dynamics in an analogous way to achieve a set 323 of three simpler, albeit still nonlinear, ODEs: where 325 E[e −µX ] = r r + µ Φ(µ/r; µ/r + 2; − Habz rM ) Φ(µ/r; µ/r + 1; − Habz rM ) The interpretation of this formulation is much more straightforward than the previous two versions 326 -we no longer need information about the entire distribution. In particular, this form makes it clear 327 that the distribution of MOI only affects the rate of recovery and it is therefore much more directly 328 comparable to other Ross-Macdonald models. However, it is important to remember this model only 329 holds when the transient distribution is exactly equal to its stationary distribution at every point in 330 time. 331

332
In this paper we introduced a model of malaria superinfection which properly accounts for 333 treatment and demography, two epidemiologically-relevant features of malaria which influence the 334 distribution of MOI and have strong effects on the mean recovery rate. As the proportion of humans 335 with a single infection are the only ones recovering naturally, these perturbations to the recovery rate 336 can propagate through to differentially affect the modeled effect size of drugs on transmission. In 337 particular, this effect will be highly context-dependent -a drug will have a stronger measured effect 338 on prevalence in areas with a high force of infection than predicted in a model without superinfection. 339 In addition, temporary protection conferred by drugs can have a major impact on the prevalence, 340 as shown in comparing figures 3 and 6. This "waiting room" for the infection queue has the effect 341 of pulling down the equilibrium prevalence, and nonlinearly impacting the influence of FOI on 342 prevalence. This has interesting and somewhat troubling implications for malaria control in areas 343 with high treatment rates -reducing the effective FOI through any intervention appears to has the 344 opposite intended effect on overall prevalence.

345
While this model represents a substantial advance and generalization of past results, it is im-346 portant to keep in mind the underlying assumptions when interpreting the model predictions in 347 real transmission systems. Because infection events are modeled as a simple Poisson process, an 348 infectious bite can only result in a single infection. This allows the original M/M/∞ model to be 349 "tridiagonal" in the sense that an individual can only increase or decrease the number of infections by 350 one at a time. However in areas of high transmission intensity, a bite from a superinfected mosquito 351 may establish multiple infections from the same blood feeding event [1]. Parasite diversity can also 352 contribute to transmission; for example, deletion of HRP-2, the target of common rapid diagnostic 353 tests, can decrease the rate of detection and subsequent treatment of infections [6]. There is also 354 evidence to suggest a host who becomes anemic due to periods of high parasitemia may be less prone 355 to infection due to the efflux of iron from the liver [11]. All of these break Macdonald's original 356 assumption of independence of the infections, and generalizations which account for these could have 357 important implications for malaria control and elimination strategies. Future explorations of these 358 factors is needed.

359
Despite these caveats, this model serves two main purposes. First, it provides a consistent 360 framework with respect to the underlying assumptions of the stochastic process for investigating the 361 impact of drugs and demographic turnover on the prevalence of any pathogen which can superinfect 362 its host. Second, it provides an analytic demonstration of how the distribution of MOI under the 363 effects of superinfection is inextricably connected to the structure of the rest of the model, as well 364 as some tools for investigating the impact of this kind of model complexity. This connection and 365 the dynamics that it generates could change the policy recommendation for intervention packages 366 in different regions, impacting the lives of all people living in endemic areas. Because of that, we 367 must be cognizant of the impact of modifications of this type of hybrid model on the structure of 368 the distribution of MOI as we include more detail and realism in our attempts to better understand 369 and prevent the spread of infectious disease.