The Potential Elimination of Plasmodium vivax Malaria by Relapse Treatment: Insights from a Transmission Model and Surveillance Data from NW India

Background With over a hundred million annual infections and rising morbidity and mortality, Plasmodium vivax malaria remains largely a neglected disease. In particular, the dependence of this malaria species on relapses and the potential significance of the dormant stage as a therapeutic target, are poorly understood. Methodology/Principal Findings To quantify relapse parameters and assess the population-wide consequences of anti-relapse treatment, we formulated a transmission model for P. vivax suitable for parameter inference with a recently developed statistical method based on routine surveillance data. A low-endemic region in NW India, whose strong seasonality demarcates the transmission season, provides an opportunity to apply this modeling approach. Our model gives maximum likelihood estimates of 7.1 months for the mean latency and 31% for the relapse rate, in close agreement with regression estimates and clinical evaluation studies in the area. With a baseline of prevailing treatment practices, the model predicts that an effective anti-relapse treatment of 65% of those infected would result in elimination within a decade, and that periodic mass treatment would dramatically reduce the burden of the disease in a few years. Conclusion/Significance The striking dependence of P. vivax on relapses for survival reinforces the urgency to develop more effective anti-relapse treatments to replace Primaquine (PQ), the only available drug for the last fifty years. Our methods can provide alternative and simple means to estimate latency times and relapse frequency using routine epidemiological data, and to evaluate the population-wide impact of relapse treatment in areas similar to our study area.

tions, because in epidemic regions mosquito abundances and malaria cases vary significantly.
Therefore we model numbers, and not fractions, of infections.
The basic SEIH n QS model (without relapse treatment) is a special case of the full model that includes treatment, and so we present below equations for the latter (see Table S1 for parameter definitions): dS/dt = (δP + dP/dt) + µ IS I + µ QS Q + (aµ IH I + bµ EI E) − µ SE S − δS, dI/dt = (1 − b)µ EI E + nµ HI H n − (µ IH + µ IS + µ IQ )I − δI, (S1) The birth rate for the S class is set to ensure that S(t)+E(t)+I(t)+Q(t)+ k H k (t) = P (t), where P (t) is Kutch population size at time t obtained from interpolated census data. The advantage of using multiple H i (i = 1, · · · n) classes, which partitions the single dormant stage into a sequence of identical stages (Fig.1A), is to introduce a flexible Gamma-distribution for the relapse intervals s (Fig.3F) given by the formula G(s) = (nµ HI ) n s n−1 exp (−nµ HI ) /Γ(n) [7], where nµ HI denotes the transition rate from one H stage to the next. Relapse treatment reduces, by a fraction a, the rate µ IH at which the infected population enters dormancy via I-to-H 1 transition, and the treated humans join non-relapsing infecteds S1 in moving to the Q class. Equations for the control (untreated) SEIH 3 QS model, with n = 3, thus readily follow by setting a = 0.
The role of mosquitoes in the transmission from infected to susceptible humans is represented implicitly through their effect in generating a distributed time delay between the current rate of transmission experienced by a susceptible human at time t, µ SE (t), and the force of infection resulting from levels of infection in the human population at all previous times [1] . Because this force of infection, λ, determines the number of infected vectors, but does not affect transmission until the parasite completes its development (sporogony) within the surviving mosquitoes over this time, we refer to λ as the "latent" force of infection. Thus, the inoculation rate µ SE (t) can be expressed as the integral over time of this force of infection, weighted by a probability that describes the delay resulting from parasite development, as follows: For the delay probability function γ(s) we choose a gamma distribution (with mean τ and variance τ 2 /k), because it allows us to implement a flexible shape, in particular one with a characteristic time scale of parasite development [7], in contrast to an exponential distribution.
Moreover, the latent force of infection is given by where q denotes a reduced infection risk from humans in the Q class, and the transmission term β(t) includes three extrinsic drivers corresponding, respectively, to (1)

t). Stochasticity is incorporated by a
Gamma distribution Γ(t) with intensity σ 2 pro (see [1,2] for the details of implementing s i and Γ). We construct the rainfall function C(t) by accumulating rain over the preceding 4 months for each reported month, and then setting January-May and November-December rain to zero. Thus, C(t) = t t−4/12 r(s)ds between July and October and = 0 for other times, where r(s) is a spline interpolation of the discrete (monthly) rainfall data. This particular choice follows from observed correlation properties between monthly vivax cases and rainfall (Fig.S1). We then normalize C(t) to zero mean and unit variance before using it in (S3).
For numerical convenience we do not implement the integral in (S2), and instead use the equivalent representation in terms of sequential transitions through a chain of identical stages, whose only purpose is to generate the desired distributed delay between λ and µ SE [7]: κ 1 , · · · , κ k (≡ µ SE ) as follows [7]: We use k = 2 as in [1,2], which reduces (S4) to equations for κ 1 and µ SE .
Although humans belonging to the E, I and Q classes can all carry liver-stage hypnozoites, we only consider the I-to-H-to-I transition to denote relapse (Fig.1A). This assumption was based on the following considerations: 1) clinical studies designed to determine important relapse patterns, such as relapse rate and frequency, typically track patients through relapsing episodes [8,9,10,11,12,13], and can be directly modeled with an I-to-H-to-I type loop; 2) model parsimony dictated by our inference goal limits the amount of unobserved details that time-series data can support [14]; and, 3) if relapse is primarily triggered by physiological stress from illness, as recently suggested [13], the E-to-H and Q-to-H transitions will mostly give non-relapsing humans who have no effect on transmission, and can therefore be safely ignored in the model. Note that this relapse loop includes the possibility of multiple relapses (albeit at a constant rate) as is common in India [9,15].
The equation for the measurement model, which couples the continuous-time dynamics

S3
of model (S1) with the discrete-time sequence y 1 , · · · , y N of monthly reported case data at times t 1 , · · · , t N , is given by, "Negbin(a, b)" is the negative binomial distribution with mean a and variance a + a 2 b (see [1] for more details). We carried out all numerical simulations in the R computing environment [16], and used the R package "Pomp" [17] to implement the algorithm for statistical inference, which is detailed elsewhere (e.g. see Supplement of [1] for a summary of the algorithms steps, and also [18]).
Equations of the SEIHQS model, with a single H class, can be obtained by removing the H k=2···n equations from the untreated version of (S1) and substituting n = 1 in (S1) and (S5). The equations for the two non-relapse models can be obtained from (S1) as follows: model SEIQS by setting µ IH = µ HI = 0, and SEIRS by setting q = µ IH = µ HI = µ IS = 0, Q≡R.