Modelling the drivers of the spread of Plasmodium falciparum hrp2 gene deletions in subSaharan Africa
 Cited 0
 Views 815
 Annotations
Abstract
Rapid diagnostic tests (RDTs) have transformed malaria diagnosis. The most prevalent P. falciparum RDTs detect histidinerich protein 2 (PfHRP2). However, pfhrp2 gene deletions yielding falsenegative RDTs, first reported in South America in 2010, have been confirmed in Africa and Asia. We developed a mathematical model to explore the potential for RDTled diagnosis to drive selection of pfhrp2deleted parasites. Low malaria prevalence and high frequencies of people seeking treatment resulted in the greatest selection pressure. Calibrating our model against confirmed pfhrp2deletions in the Democratic Republic of Congo, we estimate a starting frequency of 6% pfhrp2deletion prior to RDT introduction. Furthermore, the patterns observed necessitate a degree of selection driven by the introduction of PfHRP2based RDTguided treatment. Combining this with parasite prevalence and treatment coverage estimates, we map the modelpredicted spread of pfhrp2deletion, and identify the geographic regions in which surveillance for pfhrp2deletion should be prioritised.
https://doi.org/10.7554/eLife.25008.001eLife digest
Since the turn of the millennium, a large increase in international funding has helped to reduce the public health impact of malaria. The introduction of rapid diagnostic tests has played a central role in these efforts, particularly in remote areas that are heavily affected by the disease. These tests analyse human blood samples for specific proteins that are produced by malaria parasites.
The most common rapid diagnostic tests for malaria detect a protein called HRP2, which is produced by the deadliest malaria parasite, Plasmodium falciparum. Recently, however, cases have emerged where the tests have failed to detect these malaria infections. The first occurred in South America, and were found to be because some malaria parasites no longer possessed the gene that produces HRP2. Since then, malaria parasites that lack this gene have been found in several locations in Africa. This raises the question of whether using the tests favours the survival and spread of parasites that cannot produce the HRP2 protein.
Using mathematical modelling techniques, Watson et al. now present evidence that suggests that the use of HRP2detecting rapid diagnostic tests over the past 10 years could have favoured the evolution of malaria parasites that lack this protein. Furthermore, the models suggest that the conditions that are most likely to cause such selection are places where malaria infections are not common but people seek treatment at high rates.
Using this information, Watson et al. created a map of 160 locations in Africa most at risk of rapid diagnostic testdriven selection against the gene that produces HRP2. Public health authorities could use these maps to determine where they should more closely monitor malaria parasites to see if they lack this gene.
Future genetic investigations will be required in the highrisk areas to confirm and refine the predictions. The development of rapid diagnostic tests that detect other malaria proteins will also be essential if malaria parasites that lack HRP2 continue to spread.
https://doi.org/10.7554/eLife.25008.002Introduction
Efforts to control malaria globally have made substantial progress in the last 15 years (World Health Organization, 2015a). This progress reflects the impact made by reinvigorated political commitment that has yielded a twentyfold increase in international funding for the control and elimination of malaria (World Health Organization, 2015a). The World Health Organisation (WHO) Global Technical Strategy for Malaria 2016–2030 sets ambitious goals to further reduce incidence and mortality rates by 90% by 2030 (World Health Organization, 2015b). Central to achieving these goals is the need to test, treat and track all malaria (World Health Organization, 2010).
In subSaharan Africa (SSA), diagnostic testing of suspected malaria cases has risen from 36% to 60% between 2005 and 2014 (World Health Organization, 2015a). Microscopy was historically the most common method for diagnosis; however rapid diagnostic tests (RDTs) accounted for 71% of all diagnostic testing of suspected cases in 2014 (World Health Organization, 2015a). The most widely used RDTs target histidinerich protein 2 (HRP2), which is expressed by the Plasmodium falciparum (Pf) specific gene pfhrp2, with over 85% of RDTs tested in the WHO Foundation for Innovative New Diagnostics (FIND) Malaria RDT Evaluation Programme targeting PfHRP2 (World Health Organization, 2012a).
Falsenegative RDT results due to a partial or complete deletion of the pfhrp2 gene have been reported in areas of South America since 2010, resulting in the recommendation against the use of PfHRP2based RDTs in these areas (Akinyi et al., 2013; Abdallah et al., 2015; Cheng et al., 2014). These pfhrp2deleted mutants may still possess a functioning pfhrp3 gene; however, the cross reactivity between PfHRP2based RDT antibodies and PfHRP3 epitopes is such that a positive result may only occur at very high parasitaemia (Baker et al., 2005). Confirmed pfhrp2deleted mutants are rarer in Africa, with the first cases reported in Mali in 2012 (Koita et al., 2012). However, recently confirmed occurrences in Ghana, (Amoah et al., 2016) Zambia, (Laban et al., 2015) the Democratic Republic of Congo (DRC), (Parr et al., 2016) Rwanda (Kozycki et al., 2017) and Eritrea (Berhane et al., 2017) (Table 1) have prompted the WHO to host Technical Consultations on pfhrp2/3deletions and to issue interim guidance for malaria control programs (World Health Organization, 2017; World Health Organization, 2016a; World Health Organization, 2016b). This raises the concern that pfhrp2deleted mutants may be selected for by RDTguided treatment decisions – which if confirmed would be one of the first example of selection of a pathogen through diagnostic testing.
Here we use mathematical modelling to characterise the impact of introducing PfHRP2based RDTs on the emergence and spread of pfhrp2deleted parasites. We adapt a previously published transmission model (Griffin et al., 2016), incorporating the transmission of pfhrp2deleted mutants and the contribution of PfHRP3 crossreactivity to identify settings in which the selective pressure favouring pfhrp2deleted strains is greatest. In addition, we conduct sensitivity analyses to characterise the influence of assumptions within our model concerning adherence to RDTguided treatment decisions, the use of microscopybased diagnostic testing, fitness costs associated with the mutant parasite and the impact of nonmalarial fevers upon the selective advantage of pfhrp2 gene deletions. We continue by using a nationally representative crosssectional study of pfrhp2deletion in the DRC (Parr et al., 2016) to estimate the prevalence of pfhrp2deleted mutants prior to RDT introduction. This, in turn, allows us to map predicted geographical regions across SSA where pfhrp2deletion surveillance should be focused. These mapped predictions are explored across a range of estimates of the prevalence of pfhrp2deleted mutants prior to RDT introduction.
Results
Using our newly adapted model incorporating the transmission of pfhrp2deleted mutants, we first explored the potential for RDTguided treatment decisions to exert an evolutionary pressure on the prevalence of the mutant. Figure 1 shows the predicted proportion of strains that are pfhrp2deleted within the population after 10 years.
Within all settings that explored different transmission intensities and starting frequencies of pfhrp2deletion, RDT introduction is predicted to increase the proportion of pfhrp2deleted mutants. The strength of selection is predicted to be greatest at low PfPR (Figure 1a); however, a selective pressure is still predicted at both high PfPR and at low starting pfhrp2deletion frequencies (Figure 1b). The variance in the selection pressure exerted by RDTs is also predicted to be greatest at low PfPR (Figure 1c). A more gradual but analogous trend is predicted in the proportion of the population that were only infected with pfhrp2deleted mutants (Figure 1d). The prevalence of malaria within Figure 1a was also observed to increase after RDT introduction (Figure 1—figure supplement 1), with the greatest increase in lower transmission settings resulting from untreated infections due to falsenegative RDT results.
Within the sensitivity analyses, a selective pressure is observed to exist at comparative fitness costs of greater than 90% (see Figure 1—figure supplement 2), however below this the pfhrp2deletion allele is quickly lost. Both the introduction of additional diagnosis with microscopybased methods and nonadherence to RDT results decreased the selective pressure, slowing the rate of pfhrp2deletion emergence (see Figure 1—figure supplement 3). The introduction of nonmalarial fevers, however, increased the rate of pfhrp2deletion emergence (see Figure 1—figure supplement 4), even at 25% below the mean estimated rate of nonmalarial fever. When these opposing factors were combined, RDT introduction is still predicted to increase the proportion of pfhrp2deleted mutants (Figure 1—figure supplement 5).
The proportion of clinical cases seeking treatment (assumed here to be treated on the basis of an RDT result) is also predicted to exert a strong selection pressure for pfhrp2deletion (Figure 2). A consistent relationship was seen across comparable PfPR ranges, with the lowest treatment seeking rates (f_{T} = 0.2) yielding the slowest increase in the proportion of infections due to only pfhrp2deleted mutants. Again, the lower PfPR categories show the greatest selection pressures for pfhrp2deletion, with treatment seeking rates >30% and PfPR <25% leading to 20% of infections due to only pfhrp2deleted mutants in fewer than five years (Figure 2a).
The selection pressure favouring pfhrp2deletions is predicted to be weaker when PfHRP3 epitopes are assumed to cause positive RDT results (Figure 2b). In settings where PfHRP3 epitopes are assumed to cause a positive RDT result in 25% of cases (ε = 0.25), there are four fewer prevalence categories that reach 20% of infections due to only pfhrp2deleted mutants in fewer than five years. A similar effect is observed in the mean final frequency of pfhrp2deletion, with 64% frequency recorded after 20 years when no PfHRP3 epitope effect is assumed in comparison to 56% when ε is equal to 0.25 (Figure 2—figure supplement 1).
To estimate the starting frequency of pfhrp2deleted mutants, we used estimates of the proportion of pfhrp2deleted mutants from a national study in DRC (Parr et al., 2016) to calibrate the model. The calibration incorporated both the PfPR levels and estimates of the treatment rates in the 26 Divisions Provinciales de la Santé (DPS) that would drive selection of the mutant. We estimate a starting frequency of pfhrp2deleted P. falciparum of 6% in the DRC prior to any introduction of RDTs. The observed relationship between the proportion of infections due to pfhrp2deleted mutants and PCR PfPR among children 6–59 months of age (Figure 3a) displays a similar trend to the simulations, however with a notably steeper increase at lower prevalence. Of note, the same relationship was not predicted in the absence of selection pressure due to RDTbased treatments (i.e. purely on the basis of the variation in monoclonal infections) (Figure 3b).
Finally, using the baseline frequency estimate of 6% prior to RDT introduction, we explored 1000 different prevalence and treatment seeking rates spanning the range of estimates of the PfPR (Bhatt et al., 2015) and treatment levels across subSaharan Africa (SSA) in 2010 (Cohen et al., 2012) (Figure 4—figure supplement 2). The model output was aligned with these estimates by first administrative units (Figure 4—figure supplement 1), which enabled us to project the potential increase of the mutant strain and its impact on RDTguided treatment (Video 1). Our results suggest that 160 of 850 firstadministrative regions may have over 20% of all infections due to only pfhrp2deleted mutants by 2016 (Figure 4c). These areas, which we term of ‘high HRP2 concern’, are largely located in areas where PfPR_{210} in 2010 was less than 25% (Figure 4a). A number of other regions, classified as ‘moderate HRP2 concern’ have high treatment rates, and hence potential selective pressure, despite having comparatively higher transmission (Figure 4b). Our results also illustrate that regions with low transmission may have low HRP2 concern if the frequency of people seeking treatment is very low.
Discussion
Our results demonstrate that the key drivers of pfhrp2deletion selection are low malaria transmission and a high frequency of people seeking treatment and being correctly treated on the basis of diagnosis with a PfHRP2based RDT. Based on Africawide estimates of parasite prevalence and treatmentseeking behaviour at the time of RDTintroduction, we identified 160 firstadministrative units which we classify as ‘high HRP2 concern’. These are areas where the pfhrp2deleted strain is expected to increase in frequency over a relatively short timescale, and hence where further surveillance efforts should be concentrated.
Our results are based on calibration to a large representative survey of malaria across DRC. Due to its size and location in the centre of SSA, the DRC is arguably one of the most representative countries for endemic malaria in Africa. That the model was able to predict the observed relationship in the DRC, despite variability at a province level, provides support for the hypothesis that the variability in pfhrp2deletion frequency with transmission is driven by selection. However, in contrast to other reported surveys, the samples in this survey were primarily drawn from asymptomatic infections, and hence may not be representative of other reports of pfhrp2deletion in symptomatic cases with higher parasite density. However, it is interesting to note that our results show broad agreement with published data sets from Zambia (Laban et al., 2015) and Ghana (Amoah et al., 2016) (Table 1). In particular, our predictions confirm that the HRP2 concern would be greater in Ghana than in Southern Zambia. However, one study in Senegal found a lower prevalence of pfhrp2deletion than we predict (Wurtz et al., 2013).
A key uncertainty in predicting the potential spread of pfhrp2deletion due to RDTinduced selective pressure is the extent of use of, and adherence to, RDT results and the availability of appropriate treatment. On the one hand, if adherence to RDT results is poor (for example, with patients who show continued clinical symptoms of malaria in the absence of a positive test) or additional microscopybased detection is used (Figure 1—figure supplement 3), if appropriate treatment is not available (for example, due to stockouts), or if treatment is not fully curative (for example, due to patient nonadherence, drug resistance or fake drugs) then the spread of these deletions will be slower than predicted. On the other hand, in areas in which active case detection occurs, or in which treatment is sought for nonmalaria fevers (Figure 1—figure supplement 4), RDTbased treatment may also selectively clear nondeleted asymptomatic infections and hence increase the rate of spread of the deletion. However, when these factors, along with potential fitness costs associated with pfhrp2deletion, were investigated together we still observed an increase in pfhrp2deletion (Figure 1—figure supplement 5), which showed a similar rate of increase to that predicted by our model. Further data on RDT usage and adherence, as well as on nonmalarial fevers and the precise fitness cost of pfhrp2deletion, however, could help to refine mapping of areas of HRP2 concern.
A second uncertainty relates to the underlying prevalence of the mutation prior to RDT introduction. There is considerable variability in the estimates that have been measured both before and after RDT introduction, and it is entirely plausible that the presence of mutations could vary geographically at a range of spatial scales. However, estimating this variation is difficult given the lack of a sampling framework in reports mainly based on clinical cases and given the relatively small sample sizes. Thus our results should be interpreted not as predictions of the absolute levels of the gene deletion, but rather indicative of geographical areas in which surveillance should be focused. Similarly, these results should not be interpreted as predictions of the precise negative impact on malaria prevalence as a result of increased gene deletions (Figure 1—figure supplement 1), but illustrative of the potential impact of falsenegative test results upon malaria prevalence and the importance of alternative diagnostic methods (Figure 1—figure supplement 3). At the same time, further data collated in the coming months and years can be incorporated to iteratively update and refine our projections.
As with any modelling exercise, there are a number of important limitations. Firstly, we did not capture seasonality or any fitness cost associated with pfhrp2deletion. At a given transmission level, highly seasonal locations are likely to have a lower frequency of pfhrp2deletion in comparison to regions with perennial transmission. Seasonality could however cause substantial bottlenecks which may result in repetitive founder effects that could affect selection, resulting in either a decreased chance of pfhrp2fixation or an accelerated fixation if it occurred (Aguilée et al., 2009). In simulations incorporating a fitness cost the selection pressure was found to be considerably weaker (Figure 1—figure supplement 2). The exact fitness cost, despite being unknown, is likely subtle as our modelled fitness penalty would cause the strain to be eliminated at less than 90% comparative fitness. In addition, current theories concerning the role of PfHRP2 indicate a more minor role in heme detoxification than previously thought. Strains lacking PfHRP2 have been shown to be viable (Papalexis et al., 2001), with heme detoxification more dependent on the recently characterised haem detoxification protein (HDP) (Jani et al., 2008). Furthermore, in South America the first cases of pfhrp2deleted P. falciparum were confirmed prior to the introduction of RDTs (Gamboa et al., 2010). This suggests that these mutants may possess sufficiently high fitness such that the frequency of pfhrp2deletion is maintained in the absence of a selective advantage exerted through the use of HRP2based RDTs.
Secondly, our results depend on assumptions made regarding the contribution of PfHRP3 epitope crossreactivity and the potential for falsepositive RDT results. We found that increased crossreactivity with PfHRP3 epitopes decreases selection for pfhrp2deletion and was investigated due to confirmed observations of PfHRP2based RDTs detecting PfHRP3 epitopes at high parasitaemia (Baker et al., 2010). In simulations with no epitope effect, the model predicts the pattern in DRC well (Figure 3—figure supplement 1) and predicts a higher overall estimate of HRP2 concern (Figure 4—figure supplement 1), although the same regional patterns are identified (Figure 4—figure supplement 2). Furthermore, falsepositive RDT results would decrease the strength of the selection pressure, with pfhrp2monoclonal infections being treated. However, falsepositivity rates observed within round 6 of WHO RDT product testing were found to be low, with the median falsepositive rate on both clean negative samples and samples containing other infectious agents equal to 0%, and the overall falsepositive rate on samples containing immunological factors equal to 0.9% (World Health Organization, 2015c).
Thirdly, in the absence of systematic country introduction data, we assumed introduction of RDTs in all countries from 2010 in accordance with the WHO recommendation of testing in 2010 (World Health Organization, 2010). The precise date from region to region is undoubtedly more complex, however 2010 is a sensible estimate given the reported years at which RDTs were available at the community level in SSA by the WHO (Figure 4—figure supplement 5) (World Health Organization, 2012b). However, the ratio of testing via microscopy versus RDT is likely to have decreased over this period, and hence our estimate of RDT use (which our model assumes is 100% from introduction) is likely too high. The sensitivity of the output to this parameter is demonstrated in the data from the DRC, in which higher levels of pfhrp2deletion are observed in Kivu, an area in which RDT introduction likely occurred earlier than elsewhere in the country (Médecins Sans Frontières, 2007). Fourthly, extrapolating the starting frequency of pfhrp2deletion strains from the DRC across the rest of SSA is a clear oversimplification; however, in the absence of similar datasets, we feel it provides a reasonable first estimate. To assess the implications of this estimate, we also considered how the pattern of geographical areas that we have recommended for priority surveillance changes under different assumed starting frequencies of pfhrp2deletion (Figure 4—figure supplement 4). Despite the expected changes in the final frequency of pfhrp2deletion in these settings, the overall pattern of areas with the highest selectiondriven increase in pfhrp2deletion remains the same. A final limitation is that we assumed that treatment rates and transmission of malaria remain constant from 2010. This is clearly not the case, with 30 countries in SSA reporting a decline in prevalence from 2010 to 2015 (Bhatt et al., 2015). These combined effects, however, would presumably cause an increase in monoclonal infections and subsequent falsenegative RDTs due to pfhrp2deleted parasites.
In summary, our modelling predicts that an increased emergence of pfhrp2deleted mutants may be explained by the introduction of testing by PfHRP2based RDTs. If this is indeed the case, this would be, to our knowledge, one of the first examples of the emergence of resistance of a pathogen to a diagnostic test. The use of these RDTs will result in the greatest selection pressure in regions that have low malaria transmission and a high frequency of RDTbased treatment of clinical cases. Rapid and accurate diagnosis of P. falciparum infection, however, is essential for continued reduction in malaria transmission. In light of this, it may be sensible for public health bodies who are responding to reports of pfhrp2 gene deletions to focus surveillance in the regions we have identified as having a high HRP2 concern. This work should proceed alongside further improvement of nonHRP2based RDTs, such as those that detect lactate dehydrogenase, and the development of new alternative diagnostics.
Materials and methods
P. falciparum transmission model
An individuallevel stochastic model was developed to simulate the transmission dynamics of Plasmodium falciparum. The model is based upon previous modelling efforts (Griffin et al., 2016; Griffin et al., 2014; Griffin et al., 2015), and is described in full here before describing the extensions made with regards to PfHRP2 dynamics, and defining the parameters used and their sources. The model is implemented as stochastic individualbased model with a fixed daily time step, incorporating the necessary delay terms where mentioned, which is described in greater detail later. In overview, the transmission model considers people to exist in one of six infection states (Figure 5): susceptible (S), clinical disease (D), clinically diseased and receiving treatment (T), asymptomatic infection (A), protective state of prophylaxis (P), and subpatent infection (U).
Individuals begin life susceptible to infection (state S). At birth, individuals possess a level of maternal immunity that decays exponentially over the first 6 months. Each day individual i is probabilistically exposed to infectious bites governed by their individual force of infection (${\Lambda}_{i}$). ${\Lambda}_{i}$ is dependent on their preerythrocytic immunity, biting rate (dependent on both their age and their individual relative biting rate due to heterogeneous biting patterns in mosquitoes) and the mosquito population’s size and infectivity. Infected individuals, after a latent period of 12 days (${d}_{E}$), develop either clinical disease (state D) or asymptomatic infection (state A). This outcome is determined by their probability of acquiring clinical disease ($\varphi $), which is dependent on their clinical immunity. Individuals that develop disease have a fixed probability (f_{T}) of seeking treatment (state T), and a variable probability ($\epsilon}_{i$) that the clinical case yields a positive diagnostic result and subsequently receives treatment. ${\epsilon}_{i}$ is dependent on the assumed role of PfHRP3 epitopes, and the strain profile of infected individual i with respect to pfhrp2deleted mutants. Treated individuals are assumed to always recover, i.e. fullycurative treatment, and then enter a protective state of prophylaxis (state P) at rate r_{T}, before returning to susceptible at rate r_{S}. Individuals that did not receive treatment recover to a state of asymptomatic infection at rate r_{D}. Asymptomatic individuals progress to a subpatent infection (stage U) at rate r_{A}, before clearing infection and returning to susceptible at rate r_{U}. Additionally, superinfection is possible for all individuals in states D, A and U. Superinfected individuals who receive treatment will move to state T. Individuals who are superinfected but do not receive treatment in response to the superinfection will either develop clinical disease, thus moving to state D, or develop an asymptomatic infection and move to state A (except for individuals who were previously in state D, who will remain in state D).
The introduction of a fixed time step translates the waiting times at which individuals move from one infection state to another into a daily probability that this event occurs, with the probability drawn from the related exponential distribution. Thus the probability of a transition from state A to state B with hazard rate λ is given by:
The set of state transitions for individuals and their associated hazard rates are given below.
Process  Transition  Hazard 

Infection  S → I  ${\Lambda}_{i}(t{d}_{E})$ 
Progression of untreated disease to asymptomatic infection  D → A  ${r}_{D}=\frac{1}{{d}_{D}}$ 
Progression of asymptomatic infection to subpatent infection  A → U  ${r}_{A}=\frac{1}{{d}_{A}}$ 
Progression of subpatent infection to susceptible  U → S  ${r}_{U}=\frac{1}{{d}_{U}}$ 
Progression of treated disease to uninfected prophylactic period  T → P  ${r}_{T}=\frac{1}{{d}_{T}}$ 
Progression from uninfected prophylactic period to susceptible  P → S  ${r}_{P}=\frac{1}{{d}_{P}}$ 
Superinfection from untreated clinical disease, asymptomatic infection or subpatent infection  D → I A → I U → I  ${\Lambda}_{i}(t{d}_{E})$ 
Here state I denotes an infection state which is not tracked but which leads to either clinical disease (D), treated clinical disease (T) or asymptomatic infection (A). In the original model the probability of entering these states is determined by drawing a sequence of Bernoulli trials for each infected individual as:
For our model here, in which treatment is guided by RDTbased diagnostics, we introduce…
We assume that each person has a unique biting rate, which is the product of their relative age dependent biting rate, ${\psi}_{i}$, given by
and an assumed heterogeneity in biting patterns of mosquitoes, ζ_{i}, which we assume persists throughout their lifetime and is drawn from a lognormal distribution with a mean of 1,
where $1\rho $ is the relative biting rate at birth when compared to adults and ${a}_{0}$ represents the timescale at which the biting rate increases with age. The product of these biting rates is subsequently used to calculate an individual’s entomological inoculation rate, ${h}_{i}$, and subsequently their force of infection, which are given by
where ${\alpha}_{k}$ is the daily rate at which a mosquito takes a blood meal, ${I}_{M}$ is the size of the infected mosquito population, and ${b}_{i}$ is the probability of infection given an infectious mosquito bite.
The human population was assumed to have a maximum possible age of 100 years, with an average age of 21 years within the population yielding an approximately exponential age distribution typical of subSaharan countries. When an individual dies, they are replaced with a newborn individual whose individual biting rate due to heterogeneity in biting patterns is drawn again from a lognormal distribution with a mean of 1.
Immunity and detection functions
We consider three stages at which immunity may impact transmission:
Preerythrocytic immunity,${I}_{B}$; reduction in the probability of infection given an infectious mosquito bite.
Acquired and Maternal Clinical Immunity, ${I}_{CA}$ and ${I}_{CM}$ respectively; reduction in the probability of clinical disease given an infection due to the effects of blood stage immunity.
Detection immunity, ${I}_{D}$; reduction in the probability of detection and a reduction in the onward infectivity towards mosquitoes due to the effects of blood stage immunity.
Maternal clinical immunity is assumed to be at birth a proportion, ${P}_{M}$, of the acquired immunity of a 20 yearold and to decay at rate $\frac{1}{{d}_{M}}$.
The probabilities of infection, detection and clinical disease are subsequently created by transforming each immunity function by Hill functions. An individual’s probability of infection, ${b}_{i},$ is given by
where ${b}_{0}$ is the maximum probability due to no immunity, ${b}_{0}{b}_{1}$ is the minimum probability and ${I}_{B0}$ and ${\kappa}_{B}$ are scale and shape parameters respectively.
An individual’s probability of clinical disease, ${\varphi}_{i}$, is given by
where ${\varphi}_{0}$ is the maximum probability due to no immunity, ${\varphi}_{1}{\varphi}_{0}$ is the minimum probability and ${I}_{C0}$ and ${\kappa}_{C}$ are scale and shape parameters respectively.
An individual’s probability of being detected by microscopy when asymptomatic, ${q}_{i},$is given by
where ${d}_{1}$ is the minimum probability due to maximum immunity, and ${I}_{D0}$ and ${\kappa}_{D}$ are scale and shape parameters respectively. ${f}_{D}$ is dependent only on an individual’s age is given by
where ${f}_{D0}$ represents the timescale at which immunity changes with age, and ${a}_{D}$ and ${\gamma}_{D}$ are scale and shape parameters respectively. Lastly, ${\alpha}_{A}$ and ${\alpha}_{U}$ are parameters that determine the probability that an individual in states A and U are detectable by PCR, which are given by ${q}^{{\alpha}_{A}}$ and ${q}^{{\alpha}_{U}}$ respectively.
The contribution made by each infected individual towards the overall infectiousness of the human population towards mosquitoes is proportional to both their infectious state and their probability of detection, with a lower probability of detection assumed to correlate with a lower parasite density. Individuals who are in state D (clinically diseased), state U (subpatent infection) and state T (receiving treatment) contribute ${c}_{D}$, ${c}_{U}$ and ${c}_{T}$. In state A, infectious contribution,${c}_{A},$ is given by ${c}_{U}+\left({c}_{D}{C}_{U}\right){q}^{{\gamma}_{I}}$ where $q$ is the probability of being detected by microscopy when asymptomatic, and ${\gamma}_{I}$ is a parameter that controls how quickly infectiousness falls within the asymptomatic state.
Stochastic model equations
Given the definitions above, the full stochastic individualbased human component of the model can be formally described by its Kolmogorov forward equations. As before, let $i$index individuals in the population. Then the state of individual $i$ at time $t$ is given by $\left\{j,k,{t}_{k},l,{t}_{l},m,{t}_{m},cm,a,t\right\}$, where $a$ is age, $j$ represents infection status ($S,D,A,U,T$ or $P$), $k$ is the level of infectionblocking immunity and ${t}_{k}$ is the time at which infection blocking immunity was last boosted. Similarly, $l$ and ${t}_{l}$ denote the level and time of last boosting of clinical immunity, respectively, while $m$ and ${t}_{m}$ do likewise for parasite detection immunity, and $cm$ represents maternal immunity. Let ${\delta}_{p,q}$ denote the Kronecker delta (${\delta}_{p,q}=1$ if $p=q$ and 0 otherwise) and $\delta \left(x\right)$ denote the Dirac delta function. Defining ${P}_{i}$($j,k,{t}_{k},l,{t}_{l},m,{t}_{m},cm,a,t$) as the probability density function for individual $i$ being in state $\left\{j,k,{t}_{k},l,{t}_{l},m,{t}_{m},cm,a,t\right\}$ at time $t$, the time evolution of the system is governed by the following forward equation:
Here ${\mathcal{}\mathcal{\mathcal{O}}}_{b}$, ${\mathcal{}\mathcal{\mathcal{O}}}_{c}$ and ${\mathcal{}\mathcal{\mathcal{O}}}_{d}$ are commutative integral operators with the following action on an arbitrary density $f\left(j,k,{t}_{k},l,{t}_{l},m,{t}_{m},cm,a,t\right)$:
Finally, $\theta \left(x\right)$ is an indicator function such that $\theta \left(x\right)=1$ if $x<1$ and 0 otherwise. These functions allow the fixed periods of time in which immunities are not boosted after a previous boost to be included within the stochastic equations, while also allowing superinfection events to be incorporated.
For simulation, a discrete time approximation of this stochastic model was used, with a timestep of 1 day. For each individual $k$, $l$ and $m$ are set to zero at birth, while ${t}_{k}$, ${t}_{l}$ and ${t}_{m}$ are set to a large negative value ${T}_{big}$ (to represent never having been exposed or infected). Each immunity term increases by 1 for an individual whenever that individual receives an infectious bite ($k$), or is infected ($l$ and $m)$, if the previous boost to $k$, $l$ and $m$ occurred more than ${u}_{B}$, ${u}_{C}$ and ${u}_{D}$ days earlier, respectively. Immunity levels decay exponentially at rate ${r}_{B}$, ${r}_{CA}$ and ${r}_{ID}$, where ${r}_{B}$, ${r}_{CA}$ and $r}_{ID$ are equal to $\frac{1}{{d}_{B}}$,$\frac{1}{{d}_{CA}}$ and $\frac{1}{{d}_{ID}}$ respectively.
The stochastic model equations detailed above can be explained as follows. The first line is the total time derivative of ${P}_{i}$($j,k,{t}_{k},l,{t}_{l},m,{t}_{m},cm,a,t$). The next four lines describe the flows into states S, A, U and P due to progression through infection states.
The sixth line describes exposure to malaria that boosts preerythrocytic immunity but does not lead to an infection. The first term within the commutative integral operator ${\mathcal{}\mathcal{\mathcal{O}}}_{b}$ here considers the density of individuals who are in immunity class k1 and whose last boost to their preerythrocytic immunity was more than ${u}_{B}$ days earlier, and thus will be flowing into the considered density, k, from a lower preerythrocytic immunity. The second term in the integral will equate to 1 when considering individuals who are in immunity class k and whose last boost to their preerythrocytic immunity was less than ${u}_{B}$ days earlier and thus do not see their immunity boosted and hence remain in class k. This is needed to represent the current density of individuals in the considered density. There is no term for individuals in immunity class k whose last boost to their preerythrocytic immunity occurred more than ${u}_{B}$ days earlier as they would move out of the considered density (into class k + 1) and hence the indicator function will equate to 0 for these individuals.
The seventh line describes exposure events occurring to individuals in states S, A and U which do result in patent (bloodstage) infection, resulting in transition into states A, D or T. The force of infection acting on the density in state D is not included here but rather in the eighth line since these individuals may only move to states T or D and not A. In both the seventh and eighth lines, the commutative operators here function as described earlier. This tracks the density of individuals in immunity states one lower whose last boost to any of the three immunity types occurred a sufficient number of days earlier to flow into the considered density, while also considering the individuals already at the same immunity as the considered density to remain in their current combined infection/immunity state (when the indicator function equates to 1) or to move to a new infection/immunity state (when the indicator function equates to 0).
The ninth line (effectively a first order wave equation) represents deterministic exponential decay of the four different types of immunity. The tenth line represents the birth process. We assume a constant population size, so upon death, individuals flow into the state with no immunity and last immunity boosting times are set to ${T}_{big}$, chosen to be sufficiently early to allow immediate boosting upon exposure to infection (i.e. zero immunity other than maternal at birth). The last line shows the removal of individuals from the population through death, balancing the inflow from the previous line.
Mosquito population dynamics
The adult stage of mosquito development was modelled in a compartmental formulation. Susceptible adult mosquitoes (S_{M}) become infected at a rate which is proportional to the infectiousness of the human population lagged by ${d}_{g}$ days, which represents the delay from emergence of asexual bloodstage parasites to sexual gametocytes that contribute towards onward infectivity. The force of infection towards mosquitoes on a given day,${\Lambda}_{M},$ is represented by the sum of the contributions from each infected human, delayed by${d}_{g}$, towards the overall infectiousness of the human population, which is given by
Infected mosquitoes then pass through a latent stage (E_{M}) of duration${d}_{EM},$ before becoming infectious to humans (I_{M}). Infectious mosquitoes remain infectious until they die. The differential equations governing the adult stage of mosquitoes are given by
where ${\mu}_{M}$ is the daily death rate of adult mosquitoes, and M_{v} is the total mosquito population, i.e. $S}_{M}+{E}_{M}+{I}_{M$.
PfHRP2 dynamics
Individuals that are newly infected receive either a pfhrp2deleted mutant or a wild type, determined probabilistically by the ratio of pfhrp2deleted mutants in the contribution to onwards infectiousness governed by the human infectious population delayed by d_{EM}. An individual with clinical disease (state D), who possesses an equal number of wild type and pfhrp2deleted mutants will for example contribute ½c_{D} to both the wild type and the mutant profile of the onwards infectiousness to mosquitos. In simulations incorporating a fitness cost associated with pfhrp2deletion the contribution terms (c_{D}, c_{T}, c_{A} and c_{U}) required to calculate the contribution to the human infectious reservoir made by the deletion strains are comparatively decreased relevant to the wild type strains in order to represent an assumed decrease in parasitaemia and onward transmission. This effect would also capture if the pfhrp2deleted strain is comparatively less fit within the mosquito. This would be of importance when considering blood meals taken by mosquitoes feeding on polyclonally infected individuals, in which we would expect the fitter wildtype parasite to be probabilistically more likely to be onwardly transmitted.
If a newly infected individual is only infected with pfhrp2deleted mutants the probability that they enter the treated class is $\epsilon {f}_{T}$, where ε is equal to the cross reactivity contribution of PfHRP3 epitopes. If, however, they contain any wild type strains, ε is assumed to always equal 1. Additionally, if a subpatent individual is superinfected resulting in clinical disease, ε is assumed to equal the cross reactivity contribution of PfHRP3 epitopes if the acquired strain from superinfection is pfhrp2deleted. Analogously, ε is assumed to always equal 1 if the acquired strain was wild type. This is to reflect the inability of RDTs to detect any of the strains that were previously present within the subpatent individual.
Individuals that clear infection lose all strains, and infected individuals clear a random strain at rate $n{r}_{C}$, where $n$ is the total number of strains and ${r}_{C}$ is the rate at which one strain is cleared in a monoinfected individual, that is $\frac{1}{{d}_{A}+{d}_{U}}$. This introduces a carrying capacity on the number of strains an individual can be infected with, which scales with the transmission intensity. The multiplicity of infection and strain profile of an infected individual have no effect on the disease outcome except when the use of RDTs is introduced.
Model code availability and parameter values
The model code was developed using the R language (RRID:SCR_001905), (R Core Team, 2016) and is available with shape files and plotting scripts through an open source MIT license at https://github.com/OJWatson/hrp2malaRia (Watson, 2017). A copy is archived at https://github.com/elifesciencespublications/hrp2malaRia. The model is also written out in full as a pseudocode model using mathematical syntax (Supplementary file 1). Parameter estimates used within the model were taken from Griffin et al. (2014), (2015) and (Griffin et al., 2016) however have been included in Table 2 for clarity.
Characterising the epidemiological and clinical drivers of selection for pfhrp2deletion
The rate of pfhrp2deleted mutant emergence after the introduction of RDTs was examined across a range of malaria transmission intensities (10%, 25% and 60% parasite prevalence across all ages [PfPR]) and starting proportions of pfhrp2deleted mutants (2%, 8% and 12% mutants). For all simulations conducted, the proportion of clinically diseased cases seeking treatment was equal to 40% (f_{T} = 0.4). In all simulations ten stochastic realisations of 100,000 individuals were simulated for 60 years to reach equilibrium first, before exploring different parameter settings. RDTs were then introduced and the proportion of strains that are pfhrp2deleted recorded over the following 10 years. The effect of transmission intensities (0%–90% PfPR) was explored further by recording both the proportion of strains that are pfhrp2deleted and the proportion of individuals only infected with pfhrp2deleted mutants at 5 years after RDT introduction. In these simulations PfHRP3 epitopes were assumed to never yield a positive RDT result (ε = 0.0).
The rate of emergence was further examined under different assumptions about the proportion of people seeking treatment (f_{T} = 0.2–0.6), recording the time taken for the proportion of individuals only infected with pfhrp2deleted mutants to reach 20%. For each simulation we assumed 8% of strains carried pfhrp2deleted mutants prior to RDT introduction. We also considered the role that PfHRP3 antigens may have in the performance of PfHRP2based RDTs, assuming that 25% of individuals only infected with pfhrp2deleted mutants receive treatment due to the presence of PfHRP3 epitopes (ε = 0.25), based on an estimate of PfHRP2based RDT crossreactivity (Baker et al., 2005).
A series of analyses were additionally conducted to characterise the impact of a number of assumptions within the model. These sensitivity analyses were conducted at 20% PCR PfPR across all ages, with the proportion of clinically diseased cases seeking treatment equal to 40% (f_{T} = 0.4), and an assumed starting frequency of pfhrp2deletion equal to 10%. As before, in all simulations ten stochastic realisations of 100,000 individuals were simulated for 60 years to reach equilibrium first, before exploring different parameter settings.
We initially assessed the impact upon the strength of selection of a range of assumed comparative fitness costs associated with pfhrp2deletion (5%–100%). Secondly, we explored the impact on selection of introducing additional microscopybased diagnosis, occurring in 30% of cases in alignment with the 71% use of RDTs in 2014 (World Health Organization, 2015a), with and without an assumed nonadherence to RDT results in 10% of cases in alignment with estimated improving levels of adherence to RDT results (D’Acremont et al., 2013). Lastly, we investigated the impact of nonmalarial fevers (NMF), introducing an estimate for the annual NMF rate. This estimate was sourced by first finding household Demographic Health Surveys in SSA that surveyed whether individuals had been previously sick with a fever in the last 2 or 4 weeks and if and where they sought treatment for that fever. The resultant six surveys (Institut de Statistiques et d’Études Économiques du Burundi  ISTEEBU et al., 2013), (Institut de Statistiques et d’Études Économiques du Burundi  ISTEEBU et al., 2013) Liberia (2009 [National Malaria Control Program  NMCP/Liberia, Ministry of Health and Social Welfare/Liberia et al., 2012] and 2011 [National Malaria Control Program  NMCP/Liberia, Ministry of Health and Social Welfare et al., 2009]), Mali (2010) and Nigeria (2010 [National Population Commission  NPC/Nigeria et al., 2012] and 2015 [National Malaria Elimination Programme  NMEP/Nigeria et al., 2016]) were then subset by those that had sufficiently representative fever data across all ages, which yielded five surveys with Mali failing to be sufficiently representative at higher age ranges. The survey from Burundi was chosen for further analysis as both surveys from Liberia and Nigeria reported substantial treatment sought from drug peddlers and drug hawkers respectively, complicating inference on the clinical outcome of any treatment sought for the fever. From this survey an agebracketed annual rate of fever that led to treatment being sought was calculated, with smaller age brackets used at younger ages to capture the rapid change in fever rates at younger ages. This annual estimate was then scaled by 57% to represent the likely NMF rate, as estimated from a large scale estimate across Africa (Gething et al., 2010). We subsequently incorporated this rate to explore the impact of NMF upon selection. Within these simulations, we assume that individuals currently treated or in prophylaxis will not receive further antimalarial treatment when presenting with a NMF. Susceptible and subpatent individuals who seek treatment due to a NMF will only receive treatment due to nonadherence to test results. Lastly, asymptomatic and diseased individuals who seek treatment due to a NMF will always receive treatment, unless they are monoinfected with pfhrp2deleted parasites in which case they will only be treated due to potential PfHRP3 epitope contributions, nonadherence to RDT results or if they were diagnosed with microscopybased diagnosis.
Estimating the starting frequency and geographic spread of pfhrp2deletion
To estimate the current and future proportion of pfhrp2deleted mutants across SSA, we require a starting frequency of pfhrp2deletion. We used estimates of the proportion of pfhrp2deleted mutants (Parr et al., 2016) from the 2013–2014 DRC Demographic and Health Survey (DHS) (Meshnick et al., 2015) to infer the starting frequency before RDTs were introduced in 2010–2011, (Meshnick et al., 2015) using the weighted PCR prevalence of malaria in children aged 6–59 months (PCR PfPR 6–59 months) and the reported frequency of people seeking treatment in the 26 Divisions Provinciales de la Santé (DPS). The DHS survey was a nationally representative crosssectional study of 7137 children aged 6–59 months and 783 subjects with RDT/PCR+ results were tested using PCR assays to detect and confirm pfhrp2deletion.
We explored 50 starting frequencies between 0.1%–10%, with an assumed probability of a clinical case seeking treatment, who is only infected with pfhrp2deleted mutants, producing a positive RDT result ($\epsilon$) equal to 0.25. RDTs were assumed to be introduced in 2010 except for North and SouthKivu where the use of RDTs occurred from 2007 in the refugee camps. (Médecins Sans Frontières, 2007; United Nations High Commissioner for Refugees, 2013) For each starting frequency, ten stochastic realisations of 100,000 individuals were simulated for each DPS at malaria prevalence levels aligned to the observed weighted PCR prevalence of malaria in children aged 6–59 months for these provinces. These simulations were run for 60 years prior to the introduction of RDTs to ensure equilibrium was reached. The output from each set of simulations at a given starting frequency was smoothed using a local regression (LOESS) model, and the starting frequency identified as the set of simulations with the smallest residual sum of squares when compared to the recorded relationship from the DHS survey. It is important to highlight that due to the nonspatial nature of the model, each geographical region simulated occurs independently to neighbouring regions, i.e. there is no spatial spread of parasites between regions. Additionally, novel mutation emergence was not modelled explicitly and thus stochastic loss of the pfhrp2deletion genotype would always yield a final pfhrp2deletion frequency of 0%.
The estimated starting frequency was then used to simulate trends in the prevalence of pfhrp2deleted mutants across SSA, exploring a range of treatment coverages and transmission intensities, with ε = 0.25. These simulations considered populations of 100,000 individuals that were simulated for 20 years from 2010 to 2030, with the introduction of RDTs assumed across all regions in 2010. These outputs were matched to the mean microscopybased PfPR in 2–10 year olds (PfPR_{210}) in 2010 by first administrative unit and estimates of the proportion of cases seeking treatment from previously modelled estimates using the DHS and the Malaria Indicator Cluster Surveys (Cohen et al., 2012). The time taken for the proportion of infections due to only pfhrp2deleted mutants to reach 20% was recorded and classified to map areas of HRP2 concern under four qualitative classifications shown in Table 3.
References
 1

2
Allele fixation in a dynamic metapopulation: Founder effects vs refuge effectsTheoretical Population Biology 76:105–117.https://doi.org/10.1016/j.tpb.2009.05.003
 3
 4
 5

6
Genetic Diversity of Plasmodium falciparum Histidine‐Rich Protein 2 (PfHRP2) and Its Effect on the Performance of PfHRP2‐Based Rapid Diagnostic TestsThe Journal of Infectious Diseases 192:870–877.https://doi.org/10.1086/432010
 7
 8
 9
 10

11
WHO Informal Consultation on Fever Management in Peripheral Health Care Settings : A Global Review of Evidence and PracticeMalaria Policy Advisory Committee Meeting, 1315 March 2013.
 12
 13
 14
 15

16
Gradual acquisition of immunity to severe malaria with increasing exposureProceedings of the Royal Society B: Biological Sciences 282:20142657.https://doi.org/10.1098/rspb.2014.2657
 17
 18

19
FalseNegative Rapid Diagnostic Tests for Malaria and Deletion of the HistidineRich Repeat Region of the hrp2 Gene †The American Journal of Tropical Medicine and Hygiene 86:194–198.https://doi.org/10.4269/ajtmh.2012.100665
 20
 21

22
Improving Malaria Case Management in the Democratic Republic of the CongoWashington DC, United States: MalariaCare.

23
Demographic and Health Survey (DRCDHSII) 20132014 Supplemental Malaria ReportThe DHS Program.
 24

25
Liberia Malaria Indicator Survey 2009Monrovia, Liberia: NMCP, LISGIS, and ICF International.

26
Liberia Malaria Indicator Survey 2011Monrovia, Liberia: NMCP, LISGIS, and ICF International.

27
Nigeria Malaria Indicator Survey 2010Abuja, Nigeria: NPC/Nigeria, NMCP/Nigeria, and ICF International.
 28

29
Histidinerich protein 2 of the malaria parasite, Plasmodium falciparum, is involved in detoxification of the byproducts of haemoglobin degradationMolecular and Biochemical Parasitology 115:77–86.https://doi.org/10.1016/S01666851(01)002717

30
Pfhrp2 deleted Plasmodium falciparum parasites in the Democratic Republic of Congo: A national crosssectional surveyJournal of Infectious Diseases pp. jiw538–34.https://doi.org/10.1093/infdis/jiw538

31
R: A Language and Environment for Statistical ComputingR: A Language and Environment for Statistical Computing.
 32

33
http://www.unhcr.org/ukAccessed 2016.

34
hrp2malaRia, version 15ed5edGithub.
 35
 36
 37
 38
 39

40
Malaria Rapid Diagnostic Test Performance. Results of WHO Product Testing of Malaria RDTs: Round 6 (2014–2015), World Health OrganizationMalaria Rapid Diagnostic Test Performance. Results of WHO Product Testing of Malaria RDTs: Round 6 (2014–2015), World Health Organization.

41
FalseNegative RDT Results and Implications of New Reports of P. Falciparum HistidineRich Protein 2 / 3 Gene DeletionsWorld Health Organization.

42
WHO Malaria Policy Advisory Committee (MPAC) September 2016 MeetingWorld Health Organization.

43
WHO Malaria Policy Advisory Committee (MPAC) March 2017 MeetingWorld Health Organization.
 44
Decision letter

Ben CooperReviewing Editor; Mahidol Oxford Tropical Medicine Research Unit, Thailand
In the interests of transparency, eLife includes the editorial decision letter and accompanying author responses. A lightly edited version of the letter sent to the authors after peer review is shown, indicating the most substantive concerns; minor comments are not usually included.
Thank you for submitting your article "Modelling the drivers of the spread of Plasmodium falciparum hrp2 gene deletions in subSaharan Africa" for consideration by eLife. Your article has been reviewed by three peer reviewers, one of whom, Ben Cooper (Reviewer #1), is a member of our Board of Reviewing Editors, and the evaluation has been overseen by Prabhat Jha as the Senior Editor. The following individual involved in review of your submission has agreed to reveal their identity: Rachel Daniels (Reviewer #2).
The reviewers have discussed the reviews with one another and the Reviewing Editor has drafted this decision to help you prepare a revised submission.
Summary:
The authors adapt a previously published mathematical model of Plasmodium falciparum transmission to consider the impact of rapid diagnostic tests (RDTs) on the prevalence of pfhrp2 gene deleted mutants. Using data from the Democratic Republic of Congo they estimate the prevalence of pfhrp2deleted mutants prior to RDT introduction and identify the effects of malaria prevalence and frequency of people seeking treatment on the rate of spread on pfhrp2deleted mutants. Using parasite prevalence and treatment coverage estimates the model is used to identify areas of SubSaharan Africa of concern that should be prioritized for surveillance.
Essential revisions:
Broadly, the reviewers agreed this paper addressed an important subject, used an appropriate modelling framework and was wellpresented. However the methods were described with insufficient detail to enable a full understanding of what was done. It was also felt that the authors needed to do more to demonstrate that their conclusion about geographical areas of high concern are robust to alternative assumptions. We expand on these concerns below.
1) The Materials and methods section does not precisely define the actual model used. Instead an analogous deterministic model is described. This is problematic as there is not a single corresponding stochastic model. The authors need make it much clearer what they actually did by specifying precisely the stochastic model they actually used (there are no space limitations so this should be possible). Some aspects of the model also don't appear to be directly analogous to the defined d.e. model. For example, the equations in subsection “P. falciparum transmission model” do not allow for superinfection, but the text states that the actual model does, but there is no information on how this is done, what multiplicity of infection is allowed for etc. Another example: in the Materials and methods section it is stated that "each u term represents the time during which immunity cannot be boosted further after a previous boost ", which is straightforward to define for a stochastic model (assuming the authors mean a constant u – do they?), but doesn't have an obvious correspondence with the given PDEs. It is also not made clear in the Materials and methods section what the relationship is between this model and the models of Griffin et al., 2014,2015,2016. This should be clarified.
Additionally, assumptions made for the admin 1 predictions are not clear. For geographic predictions was epsilon 1 at all times when no RDTs are used? What about slowuptake of coverage of RDT use (or was it assumed that once RDTs are introduced there is immediate 100% coverage)? These assumptions have implications on the selection and on the conclusions about areas of concern for surveillance.
2) It would be a great help to nonmodellers and improve readability to outline the key assumptions of the model (without equations) in the Introduction.
3) In the supplement, it would add transparency if some level of validation was provided that the model produces output that matches the prev/incidence by country/admin 1 as previously used by this group
4) It appears that only one clinical case can occur per infection. Does this mean only fully curative treatment is assumed (no subcurative due to patient nonadherence, low dosing, lack of system compliance/resistance/fake drugs)? This should be clarified and the implications of this considered in the Discussion.
5) In the Results section and Discussion – "using the baseline frequency estimate of 6% prior to RDT" How do areas of concern change if this assumption was relaxed? The caveats around assuming a starting prevalence of 6% everywhere have been addressed in the Discussion, but in a somewhat simplified manner "Fourthly, extrapolating the starting frequency of pfhrp2deletion strains from the DRC across the rest of SSA is a clear oversimplification; however, in the absence of similar datasets, we feel it provides a reasonable first estimate". However, the work would be further strengthened if the authors undertook a simplified analysis to see if/how priority geographic areas for surveillance change if a lower or higher starting frequency of pfhrp2 was assumed. Along the lines of the other investigations concerning assumptions around epsilon etc. This seems particularly important if the authors agree with their own statement "Thus our results should be interpreted not as predictions of the absolute levels of the gene deletion, but rather indicative of geographical areas in which surveillance should be focused."
6) Discussion section: Concerning assumptions of nontreated RDT negative cases and coverage of RDT usage "However, the ratio of testing via microscopy versus RDT is likely to have decreased over this period, and hence our estimate of RDT use (which our model assumes is 100% from introduction) is likely too high". It is not clear that this, and adherence/lack of adherence to RDT treatment guidelines, would not have significant impact on the main conclusions. There is a strong argument for additional sensitivity analysis here.
7) Figure 4 and Discussion: the classification of areas into high, moderate, slight and marginal concern is based on the frequency of pfhrp2deleted mutants. But how does this frequency correlate with prevalence and ability to detect? From the results, it appears that low prevalence areas with high treatment rates are likely to increase in pfhrp2deleted mutant frequency first, but these areas may not correspond to areas that are easy to detect/sample from nor, more importantly, areas where clinical morbidity will increase as infections go untreated as they result in negative RDTs. Can the authors comment on this and address the impact on Pf prevalence in relative terms considering the dynamic feedback with treatment and transmission?
8) It appears that the selection pressure for deletion mutants arises solely as a result of differential treatment rates of patients who are clinically diseased as a result of malaria. However, in high prevalence areas, many of those treated for malaria may be parasite positive, have a clinical diagnosis of malaria, but have a nonmalarial cause of fever. In some settings these nonmalarial causes could account for the vast majority of patients treated for malaria (for example Crump JA, Morrissey AB, Nicholson WL, Massung RF, Stoddard RA, Galloway RL, et al. (2013) Etiology of Severe Nonmalaria Febrile Illness in Northern Tanzania: A Prospective Cohort Study. PLoS Negl Trop Dis 7(7): e2324. doi:10.1371/journal.pntd.0002324). Intuitively these effects could substantially change the selective pressures acting on the deletion mutants. If this process has been ignored (which is the reviewers' understanding) at the very least there should be a good argument explaining why it can be neglected. Additional analysis allowing for the fact that many patients treated for malaria may have other causes of fever would be the preferred option.
9) A related question is what happens to the patients who do have clinical disease resulting from malaria, but who don't get treated because they are infected only with the pfhrp2 deletion mutants (and don't have crossreactivity from pfhrp3). It appears from the model description these patients will never get treated for the current infection. If this is correct, how clinically realistic is this? In practice wouldn't a patient who continues to show clinical signs of malaria get treated for malaria even in the absence of a positive test result?
[Editors' note: further revisions were requested prior to acceptance, as described below.]
Thank you for resubmitting your work entitled "Modelling the drivers of the spread of Plasmodium falciparum hrp2 gene deletions in subSaharan Africa" for further consideration at eLife. Your revised article has been favorably evaluated by Prabhat Jha (Senior editor) and a Reviewing editor.
The manuscript has been improved but there are some remaining issues that need to be addressed before acceptance, as outlined below:
The authors have done a very good job in responding to the reviewers' comments, and the changes have greatly strengthened this work and have addressed the concerns raised.
The only outstanding issues for me are what I presume are typos in one or two equations (and in one case, a clarification is needed). These equations were not given in the original submission, but they are needed for a full description of the work.
Substantive concerns:
1) Something is wrong with the second equation in the Materials and methods section. First, this gives a biting rate that goes down with age, not up (contradicting the text which follows the third equation). Secondly it allows for negative biting rates (and will give them at default parameter values). Minus sign missing?
2) For the hazards for the first three transitions in subsection “Immunity and detection functions it wasn't clear to me why the period u_. (corresponding to the time when immunity cannot be further boosted after a previous boost) was accounted for by scaling the hazard rather than allowing it change with time. Wouldn't it be straightforward to specify a hazard of h_i if no previous boost in the previous u_B days for person i (or u_C or u_D) and 0 otherwise?
3) Two things that seem odd in the Kolmogorov equations in subsection “Stochastic Model Equations”: i) I can't make sense of the first term on the rhs (P_i(j,k,l,m,a,t). Is there a mu missing here? ii) on the fifth line of these equations shouldn't the $(1f_t) epsilon$ term be $(1f_t epsilon)$ ?
[Editors' note: further revisions were requested prior to acceptance, as described below.]
Thank you for resubmitting your work entitled "Modelling the drivers of the spread of Plasmodium falciparum hrp2 gene deletions in subSaharan Africa" for further consideration at eLife. Your revised article has been favorably evaluated by Prabhat Jha (Senior editor) and one reviewer, who is also a member of our Board of Reviewing Editors.
The manuscript has been improved and addressed the previous concerns but there are some remaining issues that need to be addressed before acceptance, as outlined below in the review.
Reviewer 1:
This revision has introduced a further change in the way the model is described. This is a good thing if it aids reproducibility and transparency. The one problem is that the new formulation is very difficult to make sense of, and the logic of the revised equation is not clear.
Some text should be added below these equations to explain in a less technical way what these terms on the r.h.s. represent. It may well be that this equation as it stands is correct. Or there could be one or more mistakes which need to be fixed.
Apart from this issue, I don't have any further concerns with the paper.
[Editors' note: further revisions were requested prior to acceptance, as described below.]
Thank you for resubmitting your work entitled "Modelling the drivers of the spread of Plasmodium falciparum hrp2 gene deletions in subSaharan Africa" for further consideration at eLife. Your revised article has been evaluated by Prabhat Jha (Senior editor) and a Reviewing editor.
Unfortunately, the revisions have not sufficiently addressed the concerns raised in the previous review and the issues raised there still need to be addressed before acceptance. Given the lengthy back and forth on this review, we must insist that this will be the last opportunity we will allow for a successful response to the concerns expressed by the Reviewing editor.
In particular, a precise, intelligible description of the model is an essential requirement. While the new mathematical description of the model is precise, it is not fully intelligible. This may be because there are still errors in the description, or it may because sufficient motivation for equations has not been provided. To take just one example, on the sixth line of the equations in subsection “Stochastic Model Equations” there is a term that seems to correspond to bites from infectious mosquitos which do not result in infection (the bites occur at rate h_i, the EIR, and (1b_i) is the probability of no infection given such a bite). Suppose a susceptible gets bitten, so j is S and t_k< t < t_k + u_b), so the function theta evaluates to 1. So in this case the operator b evaluates to the identify function, and plugging this back into the equation at the top of the page we learn that the rate P_i(S,…) is changing w.r.t time (keeping age constant) for susceptible hosts is (1b_i) * h_i (..) * P_i(S,…).
Why should this be the case? Why are these noninfecting bites leading to increases in the susceptible population? This is not clear and is merely the first place where the equations in this subsection cannot be readily understood.
We propose two possible ways forward:
i) a detailed explanation of the equations is added, explaining how they represent the model going into a similar level of detail as the text above.
ii) the simulation model is described in pseudocode rather than in equations. We appreciate that R code is already provided on Github. However, R will not be universally readable (either by other researchers now or in the future who are not familiar with R [and the R language is itself evolving so R code that is valid now may not run in the future]). If the model is described in pseudocode that should allow the twin aims of precision and intelligibility.
Option ii) might be easier, and of course (since space is not an issue) both i) and ii) could be done together.
https://doi.org/10.7554/eLife.25008.033Author response
Essential revisions:
Broadly, the reviewers agreed this paper addressed an important subject, used an appropriate modelling framework and was wellpresented. However the methods were described with insufficient detail to enable a full understanding of what was done. It was also felt that the authors needed to do more to demonstrate that their conclusion about geographical areas of high concern are robust to alternative assumptions. We expand on these concerns below.
1) The Materials and methods section does not precisely define the actual model used. Instead an analogous deterministic model is described. This is problematic as there is not a single corresponding stochastic model. The authors need make it much clearer what they actually did by specifying precisely the stochastic model they actually used (there are no space limitations so this should be possible).
We thank the reviewers for raising the need for further clarity pertaining to the precise numerical implementation of the stochastic model. We have redrafted the Materials and methods to fully specify the stochastic model in a way that we hope will still be accessible to those who do not follow the model equations. To do this, we have removed the deterministic descriptions and provided a list of transitions for the infection states in the human model and their associated hazards. In the section describing the acquisition and loss of naturallyacquired immunity, we have similarly provided a table with the transitions and associated hazards. We also explicitly include the standard equation to translate a continuoustime hazard rate into a daily transition probability since we implement the stochastic model on a daily time step. For mathematical clarity, we formalise this by presenting the Kolmogorov forward equation for the human infection states.
Some aspects of the model also don't appear to be directly analogous to the defined d.e. model. For example, the equations in subsection “P. falciparum transmission model” do not allow for superinfection, but the text states that the actual model does, but there is no information on how this is done, what multiplicity of infection is allowed for etc.
We thank the reviewers for highlighting the need for extra clarity and information here. In our revisions the superinfection process is clearly included in the equations and the tables of transitions.
In response to comments regarding how the multiplicity of infection is allowed for, we are deeply grateful for the reviewers bringing this to our attention as an additional paragraph detailing this was erroneously omitted in the earlier submission. This paragraph has been added in to the end of the PfHRP2 dynamics section of the transmission model, and reads as follows:
“Individuals that clear infection lose all strains, and infected individuals clear a random strain at rate nr_{C}, where n is the total number of strains and r_{C} is the rate at which one strain is cleared in a monoinfected individual, i.e. 1/(d_{A} + d_{U}). This introduces a carrying capacity on the number of strains an individual can be infected with, which scales with the transmission intensity. The multiplicity of infection and strain profile of an infected individual have no effect on the disease outcome except when the use of RDTs is introduced.”
We have not added the strains into the forward equations as this would make them very difficult to read but feel that this text description is now sufficiently clear to be reproducible.
Another example: in the Materials and methods section it is stated that "each u term represents the time during which immunity cannot be boosted further after a previous boost ", which is straightforward to define for a stochastic model (assuming the authors mean a constant u – do they?), but doesn't have an obvious correspondence with the given PDEs.
This was originally in the deterministic equivalent equations but has now been included in the forward equations. To make this clearer in the text this sentence now reads as follows: “Each u term represents the time during which immunity cannot be boosted further after a previous boost and each d term represents the duration of immunity. In the stochastic implementation of the model used here, each immunity term increases by 1 for an individual whenever that individual receives an infectious bite or is infected respectively, if the previous boost to I_{B}, I_{CA} and I_{D} occurred more than u_{B}, u_{C}, and u_{D} days earlier respectively.”
It is also not made clear in the Materials and methods section what the relationship is between this model and the models of Griffin et al., 2014,2015,2016. This should be clarified.
The basic disease model is the same as in Griffin et al., 2016, which is an amalgamation of previous iterations of this model in 2014 and 2015. The model is extended here to include the dynamics of pfhrp2deleted parasites, through tracking the infection history of individuals and introducing assumptions about the impact of these parasites upon HRP2based RDT guided treatment decisions. To make this clearer we have added the following to the beginning of the opening paragraph of the Materials and methods: “An individuallevel stochastic model was developed to simulate the transmission dynamics of Plasmodium falciparum. The model is based upon previous modelling efforts, and is described in full here before describing the extensions made with regards to PfHRP2 dynamics, and defining the parameters used and their sources.”
Additionally, assumptions made for the admin 1 predictions are not clear. For geographic predictions was epsilon 1 at all times when no RDTs are used? What about slowuptake of coverage of RDT use (or was it assumed that once RDTs are introduced there is immediate 100% coverage)? These assumptions have implications on the selection and on the conclusions about areas of concern for surveillance.
We thank the reviewers for these comments and have added the following to the Materials and methods to increase the clarity detailing the simulation assumptions for the Africa admin 1 predictions: “These simulations considered populations of 100,000 individuals that were simulated for 20 years from 2010 to 2030, with the introduction of RDTs assumed across all regions in 2010.” We fully agree with the implications related to the latter half of this comment, and have extended our sensitivity analysis accordingly and increased our discussion surrounding this within the Discussion section. Full details of this extension are made in response to comment 4 within the essential revisions where this point is raised again.
2) It would be a great help to nonmodellers and improve readability to outline the key assumptions of the model (without equations) in the Introduction.
We fully agree with these comments, and have extended the Introduction to highlight the assumptions detailing RDT nonadherence, microscopy use, fitness costs and the impact of nonmalarial fevers. This was added as follows:
“In addition, we conduct sensitivity analyses to characterise the influence of assumptions within our model concerning adherence to RDTguided treatment decisions, the use of microscopybased diagnostic testing, fitness costs associated with the mutant parasite and the impact of nonmalarial fevers upon the selective advantage of pfhrp2 gene deletions.”
3) In the supplement, it would add transparency if some level of validation was provided that the model produces output that matches the prev/incidence by country/admin 1 as previously used by this group
The model is set up to exactly match the prevalence at a given timepoint in any location through scaling the mosquito density M_{V}. For transparency we have now included an additional figure supplement (Figure 4—figure supplement 1), which shows the aligned model output compared to the estimates of prevalence within SSA produced by the Malaria Atlas Project. We have made reference to this extra figure in the final paragraph of the Results, the opening two sentences of which now read: “Finally, using the baseline frequency estimate of 6% prior to RDT introduction, we explored 1000 different prevalence and treatment seeking rates spanning the range of estimates of the PfPR and treatment levels across subSaharan Africa (SSA) in 2010 (see Figure 4—figure supplement 2). The model output was aligned with these estimates by first administrative units (see Figure 4—figure supplement 1), which enabled us to project the potential spread of the mutant strain and its impact on RDTguided treatment (see Video 1).”
The relationship between prevalence and clinical incidence is identical to our previous outputs since the model is an exact replica of this in the absence of RDT mutations.
4) It appears that only one clinical case can occur per infection. Does this mean only fully curative treatment is assumed (no subcurative due to patient nonadherence, low dosing, lack of system compliance/resistance/fake drugs)? This should be clarified and the implications of this considered in the discussion.
We agree that subcurative treatment is an important consideration but this is not captured in our current model. We have made clear our assumption about all treatment being fullycurative within the substantial rewrite of the second paragraph of the Materials and methods section. Subcurative treatment would likely slow the rate of selection, along with a number of other factors detailed in this section, which we have introduced into the 3^{rd} paragraph of the Discussion section.
5) In the Results section and Discussion – "using the baseline frequency estimate of 6% prior to RDT" How do areas of concern change if this assumption was relaxed? The caveats around assuming a starting prevalence of 6% everywhere have been addressed in the Discussion, but in a somewhat simplified manner "Fourthly, extrapolating the starting frequency of pfhrp2deletion strains from the DRC across the rest of SSA is a clear oversimplification; however, in the absence of similar datasets, we feel it provides a reasonable first estimate". However, the work would be further strengthened if the authors undertook a simplified analysis to see if/how priority geographic areas for surveillance change if a lower or higher starting frequency of pfhrp2 was assumed. Along the lines of the other investigations concerning assumptions around epsilon etc. This seems particularly important if the authors agree with their own statement "Thus our results should be interpreted not as predictions of the absolute levels of the gene deletion, but rather indicative of geographical areas in which surveillance should be focused."
We thank the reviewers for this useful comment. To address this we considered 3 additional different starting frequencies of pfhrp2deletion (3%, 4.5% and 7.5%) with an assumed probability of a clinical case seeking treatment, who is only infected with pfhrp2deleted mutants, producing a positive RDT result (ϵ)equal to 0.25. The final frequency of pfhrp2deletion across SSA was subsequently recorded and admin level 1 regions ranked accordingly, as shown in the new Figure 4—figure supplement 4. This extra sensitivity analysis is discussed within the penultimate paragraph of the Discussion as follows:
“Fourthly, extrapolating the starting frequency of pfhrp2deletion strains from the DRC across the rest of SSA is a clear oversimplification; however, in the absence of similar datasets, we feel it provides a reasonable first estimate. […]Despite the expected changes in the final frequency of pfhrp2deletion in these settings, the overall pattern of areas with the highest selectiondriven increase in pfhrp2deletion remains the same.”
6) Discussion section: Concerning assumptions of nontreated RDT negative cases and coverage of RDT usage "However, the ratio of testing via microscopy versus RDT is likely to have decreased over this period, and hence our estimate of RDT use (which our model assumes is 100% from introduction) is likely too high". It is not clear that this, and adherence/lack of adherence to RDT treatment guidelines, would not have significant impact on the main conclusions. There is a strong argument for additional sensitivity analysis here.
We agree fully with this comment, and have conducted a series of additional sensitivity analyses in response to this comment and comment 8 within the essential revisions in relation to the impact of nonmalarial fevers. With regards to nonadherence and the use of microscopybased diagnosis we conducted a similar suite of simulations to those in Figure 1, in which we incorporated a probability of diagnosis via microscopy occurring and a probability of nonadherence to RDT results occurring. The methodology for this is expanded within the new second paragraph within the Materials and methods, which also explores the introduction of nonmalarial fevers. The results of this are introduced in the new third paragraph of the Results section.
These factors, as expected, caused a decrease in the selective pressure (Figure 1—figure supplement 3). This decrease, along with the decrease due to potential fitness costs sit in opposition to an increase observed due to nonmalarial fever (discussed further in response to comment 8). As such we would want to stress our acknowledgment of these opposing factors, as discussed at greater length in the third paragraph of the Discussion, and highlight again a need for further data to allow our predictions of geographical areas of HRP2 Concern to be refined. However, with the data currently available we believe our model gives a good overall representation of how these factors might average out, as represented by the model’s ability to predict the relationship within DRC.
7) Figure 4 and Discussion: the classification of areas into high, moderate, slight and marginal concern is based on the frequency of pfhrp2deleted mutants. But how does this frequency correlate with prevalence and ability to detect? From the results, it appears that low prevalence areas with high treatment rates are likely to increase in pfhrp2deleted mutant frequency first, but these areas may not correspond to areas that are easy to detect/sample from nor, more importantly, areas where clinical morbidity will increase as infections go untreated as they result in negative RDTs. Can the authors comment on this and address the impact on Pf prevalence in relative terms considering the dynamic feedback with treatment and transmission?
We thank the reviewers for these comments, as the interaction between falsenegative test results, treatment and transmission intensity has important clinical implications resulting from changes in the frequency of pfhrp2deleted mutants. Firstly, to demonstrate this we looked at how the malaria prevalence in Figure 1A was changing as a result of increases in pfhrp2deletion due to the introduction of RDTguided treatment decisions (Figure 1—figure supplement 1). These simulations predict that as a result of increasing pfhrp2deletion, which will be greatest at the lowest transmission intensities, the prevalence of malaria increases due to an increased number of falsenegative test results. We have added these findings to the end of the second paragraph of the Results section as follows: “The prevalence of malaria within Figure 1A was also observed to increase after RDT introduction (Figure 1—figure supplement 1), with the greatest increase in lower transmission settings resulting from untreated infections due to falsenegative RDT results”. We have also incorporated these implications within the fourth paragraph of the Discussion, reiterating the need for careful interpretation of these results as follows: “Similarly, these results should not be interpreted as predictions of the precise negative impact on malaria prevalence as a result of increase gene deletions (Figure 1—figure supplement 1), but illustrative of the potential impact of falsenegative test results upon prevalence and the importance of alternative diagnostic methods (Figure 1—figure supplement 3).”
We also agree that low transmission areas where this effect is greatest may not be the easiest to sample from, for example due to the absence of adequate sampling infrastructure as a result of the areas historical low transmission. If these areas, however, have a high risk of pfhrp2deletion due to a high frequency of treatment it seems likely that the infrastructure in these settings would be sufficient. We also agree that areas with low prevalence and high treatment coverage will have fewer cases so the impact on clinical morbidity will likely be lower due to stronger health care systems and the ability to notice these falsenegative systems. For example, a recent report from Eritrea (Berhane A, Russom M, Bahta I, Hagos F, Ghirmai M, Uqubay S. 2017. Rapid diagnostic tests failing to detect Plasmodium falciparum infections in Eritrea: an investigation of reported false negative RDT results. Malar J 16: 105. doi:10.1186/s1293601717529), shows that although the effect on morbidity may not be severe, due to their extensive use of RDTs from 2007, within settings with a high treatment seeking rate, a substantial prevalence of pfhrp2deletion has now been observed. Conversely, settings with the highest clinical cases will also have the highest multiplicities of infection, so are less likely to experience negative RDT results as clinical cases are unlikely to be monoinfected.
8) It appears that the selection pressure for deletion mutants arises solely as a result of differential treatment rates of patients who are clinically diseased as a result of malaria. However, in high prevalence areas, many of those treated for malaria may be parasite positive, have a clinical diagnosis of malaria, but have a nonmalarial cause of fever. In some settings these nonmalarial causes could account for the vast majority of patients treated for malaria (for example Crump JA, Morrissey AB, Nicholson WL, Massung RF, Stoddard RA, Galloway RL, et al. (2013) Etiology of Severe Nonmalaria Febrile Illness in Northern Tanzania: A Prospective Cohort Study. PLoS Negl Trop Dis 7(7): e2324. doi:10.1371/journal.pntd.0002324). Intuitively these effects could substantially change the selective pressures acting on the deletion mutants. If this process has been ignored (which is the reviewers' understanding) at the very least there should be a good argument explaining why it can be neglected. Additional analysis allowing for the fact that many patients treated for malaria may have other causes of fever would be the preferred option.
We believe this a very important question, and are indebted to the reviewers for raising it. As mentioned in response to comment 6, we conducted additional sensitivity analyses in response to this comment. In brief, we attempted to estimate an annual rate of nonmalarial fever a different age ranges to capture the increased variability in fevers at lower ages. This was conducted by using all available Demographic Health Surveys that recorded data regarding previous fevers, and had sufficiently low rates of treatment sought from drug peddlers or hawkers. This was then averaged across discrete age brackets, with more age brackets at the lower ages and combined with the estimated proportion of fevers that are nonmalarial fevers from across Africa in 2007 (Gething PW, Kirui VC, Alegana VA, Okiro EA, Noor AM, Snow RW. 2010. Estimating the number of paediatric fevers associated with malaria infection presenting to Africa’s public health sector in 2007. PLoS Med 7. doi:10.1371/journal.pmed.1000301). We elected to use this estimate as we hoped it would be more representative of the range of nonmalarial fevers across Africa.
This inclusion of nonmalarial fevers created a substantial increase in the selective pressure (Figure 1—figure supplement 4), and if we had incorporated the proportion of fevers that are nonmalarial from the cited study from Tanzania, this increase would have been even greater. This increase was notably greater than individual decreases to the selective pressure from potential microscopy use or nonadherence to RDT results, the results of which are introduced in the new third paragraph of the Results section.
As mentioned earlier, we want to reiterate our acknowledgment of these opposing factors, and highlight that the presented model managed to capture the trend shown in DRC, and yielded a selection pressure that was very similar to that observed when considering the intermediate levels of the model assumptions together within our sensitivity analysis (Figure 1—figure supplement 5). We are aware, however, that in other regional settings these factors might lead to potentially a slower or faster rate of selectiondriven increase in pfhrp2deletion. An example of this can be seen in the recent report from Eritrea (Berhane A, Russom M, Bahta I, Hagos F, Ghirmai M, Uqubay S. 2017. Rapid diagnostic tests failing to detect Plasmodium falciparum infections in Eritrea: an investigation of reported false negative RDT results. Malar J 16: 105. doi:10.1186/s1293601717529), where RDTs were deployed in 2007 within settings with a high treatment seeking rate and a substantial prevalence of pfhrp2deletion has now been observed. We also have mentioned, however, that our model over predicted the likely prevalence of pfhrp2deletion in Senegal, which again highlights that our results should be interpreted not as predictions of the absolute levels of the gene deletion, but rather indicative of geographical areas in which surveillance should be focused. We have included a summary of these arguments within bot the third and fourth paragraphs of the Discussion, included the recent Eritrean study as well as a recent study from Rwanda within Table 1, and would like to thank the reviewers again for raising this important point.
9) A related question is what happens to the patients who do have clinical disease resulting from malaria, but who don't get treated because they are infected only with the pfhrp2 deletion mutants (and don't have crossreactivity from pfhrp3). It appears from the model description these patients will never get treated for the current infection. If this is correct, how clinically realistic is this? In practice wouldn't a patient who continues to show clinical signs of malaria get treated for malaria even in the absence of a positive test result?
We thank the reviewers for raising this point. In the model, a clinically diseased individual who does not receive treatment will move to an asymptomatic infection on average after 5 days, before potentially recovering via a subpatent infection, or being infected again at which point they may be treated. We agree that it is plausible that patients who show continued clinical symptoms of malaria in the absence of a positive test result may be treated, either through an alternative diagnostic method being used, or due to nonadherence to the test result. These potential routes to treatment are considered in response to earlier comments made by the reviewers. The extensions to the third paragraph of the Discussion raises this point, explaining how these routes may capture the clinical realism of patients being treated without a positive RDT result due to their continued presentation of clinical symptoms of malaria.
[Editors' note: further revisions were requested prior to acceptance, as described below.]
The only outstanding issues for me are what I presume are typos in one or two equations (and in one case, a clarification is needed). These equations were not given in the original submission, but they are needed for a full description of the work.
Substantive concerns:
1) Something is wrong with the second equation in the Materials and methods section. First, this gives a biting rate that goes down with age, not up (contradicting the text which follows the third equation). Secondly it allows for negative biting rates (and will give them at default parameter values). Minus sign missing?
Thank you for highlighting this omission. The minus sign was missing and has since been added.
2) For the hazards for the first three transitions in subsection “Immunity and detection functions it wasn't clear to me why the period u_. (corresponding to the time when immunity cannot be further boosted after a previous boost) was accounted for by scaling the hazard rather than allowing it change with time. Wouldn't it be straightforward to specify a hazard of h_i if no previous boost in the previous u_B days for person i (or u_C or u_D) and 0 otherwise?
These hazard rates have now been removed, as the extended Kolmogorov equations now better describe the acquisition of immunity using a series of commutative integral operators with indicator functions to represent the delay in further boosting.
3) Two things that seem odd in the Kolmogorov equations in subsection “Stochastic Model Equations”: i) I can't make sense of the first term on the rhs (P_i(j,k,l,m,a,t). Is there a mu missing here? ii) on the fifth line of these equations shouldn't the $(1f_t) epsilon$ term be $(1f_t epsilon)$ ?
Thank you for highlighting these. The first term on the rhs has been removed in the extended equations, and the epsilon term has been moved inside the bracket.
[Editors' note: further revisions were requested prior to acceptance, as described below.]
Reviewer 1:
This revision has introduced a further change in the way the model is described. This is a good thing if it aids reproducibility and transparency. The one problem is that the new formulation is very difficult to make sense of, and the logic of the revised equation is not clear.
Some text should be added below these equations to explain in a less technical way what these terms on the r.h.s. represent. It may well be that this equation as it stands is correct. Or there could be one or more mistakes which need to be fixed.
Apart from this issue, I don't have any further concerns with the paper.
We agree the extensions made to the stochastic equations could make the reproducibility of the model more difficult. However, we have included these equations in response to earlier requests made by the reviewers to include a single model that specified precisely the stochastic model. The extensions made in these revised equations now allow the period of time in which immunity is not boosted to be described formulaically, which addresses the reviewer’s earlier essential comments concerning the issue of describing an analogous deterministic model.
We do, however, believe the lines of the model mentioned by the reviewer are hard to interpret. Additional text has been added after these equations that now reads: “These functions allow the fixed periods of time in which immunities are not boosted after a previous boost to be included within the stochastic equations, while also allowing superinfection events to be incorporated.” We also want to mention again that the simulation model has been made publicly available within the Github repository listed within the Materials and methods to help reproducibility.
[Editors' note: further revisions were requested prior to acceptance, as described below.]
In particular, a precise, intelligible description of the model is an essential requirement. While the new mathematical description of the model is precise, it is not fully intelligible. This may be because there are still errors in the description, or it may because sufficient motivation for equations has not been provided. To take just one example, on the sixth line of the equations in subsection “Stochastic Model Equations” there is a term that seems to correspond to bites from infectious mosquitos which do not result in infection (the bites occur at rate h_i, the EIR, and (1b_i) is the probability of no infection given such a bite). Suppose a susceptible gets bitten, so j is S and t_k< t < t_k + u_b), so the function theta evaluates to 1. So in this case the operator b evaluates to the identify function, and plugging this back into the equation at the top of the page we learn that the rate P_i(S,…) is changing w.r.t time (keeping age constant) for susceptible hosts is (1b_i) * h_i (..) * P_i(S,…).
Why should this be the case? Why are these noninfecting bites leading to increases in the susceptible population? This is not clear and is merely the first place where the equations in this subsection cannot be readily understood.
We propose two possible ways forward:
i) a detailed explanation of the equations is added, explaining how they represent the model going into a similar level of detail as the text above.
ii) the simulation model is described in pseudocode rather than in equations. We appreciate that R code is already provided on Github. However, R will not be universally readable (either by other researchers now or in the future who are not familiar with R [and the R language is itself evolving so R code that is valid now may not run in the future]). If the model is described in pseudocode that should allow the twin aims of precision and intelligibility.
Option ii) might be easier, and of course (since space is not an issue) both i) and ii) could be done together.
In response to these valid concerns we have decided to carry out both suggested ways forward. In response to the request for model pseudocode we have written out the full model using a mathematical style syntax with detailed annotations that describing and explaining the rationale behind each line of code, including both a pseudocode syntax and colour key to aid understanding of the syntax. This has been included as a supplementary file (Supplementary file 1), and has been referenced to within the main text within the “Model code availability and parameter values” paragraph of the Materials and methods. We have added this as a supplementary file due to the considerable length of the pseudocode model, and to hopefully preserve the monospaced font and colour syntax, which aid the intelligibility of the pseudocode.
In addition we have added four new paragraphs following the stochastic model equations that provides a detailed explanation of the motivation for each line of equations. This new text reads as follows:
“The stochastic model equations detailed above can be explained as follows. The first line is the total time derivative of Pi (j,k,tk,l,tl,m,tm,cm,a,t). The next four lines describe the flows into states S, A, U and P due to progression through infection states. […] The last line shows the removal of individuals from the population through death, balancing the inflow from the previous line.”
https://doi.org/10.7554/eLife.25008.034Article and author information
Author details
Funding
Wellcome (109312/Z/15/Z)
 Oliver John Watson
National Institute of Allergy and Infectious Diseases (5R01AI107949)
 Steven R Meshnick
Imperial College London
 Hannah C Slater
Medical Research Council (MR/N01507X/1)
 Robert Verity
Department for International Development
 Azra C Ghani
National Institute of Allergy and Infectious Diseases (5T32AI007151)
 Jonathan B Parr
The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.
Reviewing Editor
 Ben Cooper, Reviewing Editor, Mahidol Oxford Tropical Medicine Research Unit, Thailand
Publication history
 Received: January 9, 2017
 Accepted: August 21, 2017
 Accepted Manuscript published: August 24, 2017 (version 1)
 Accepted Manuscript updated: August 29, 2017 (version 2)
 Version of Record published: September 18, 2017 (version 3)
 Version of Record updated: September 29, 2017 (version 4)
Copyright
© 2017, Watson et al.
This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.
Metrics

 815
 Page views

 146
 Downloads

 0
 Citations
Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.