(C) PLOS One This story was originally published by PLOS One and is unaltered. . . . . . . . . . . Population heterogeneity in Plasmodium vivax relapse risk [1] ['Eva Stadler', 'The Kirinstitute', 'Unsw Sydney', 'Sydney', 'Deborah Cromer', 'Somya Mehra', 'School Of Mathematics', 'Statistics', 'University Of Melbourne', 'Melbourne'] Date: 2022-12 A key characteristic of Plasmodium vivax parasites is their ability to adopt a latent liver-stage form called hypnozoites, able to cause relapse of infection months or years after a primary infection. Relapses of infection through hypnozoite activation are a major contributor to blood-stage infections in P vivax endemic regions and are thought to be influenced by factors such as febrile infections which may cause temporary changes in hypnozoite activation leading to ‘temporal heterogeneity’ in reactivation risk. In addition, immunity and variation in exposure to infection may be longer-term characteristics of individuals that lead to ‘population heterogeneity’ in hypnozoite activation. We analyze data on risk of P vivax in two previously published data sets from Papua New Guinea and the Thailand-Myanmar border region. Modeling different mechanisms of reactivation risk, we find strong evidence for population heterogeneity, with 30% of patients having almost 70% of all P vivax infections. Model fitting and data analysis indicates that individual variation in relapse risk is a primary source of heterogeneity of P vivax risk of recurrences. Despite elimination efforts, malaria continues to be a public health burden world-wide. Partially due to its ability to remain dormant in the liver for weeks or months, the malaria parasite Plasmodium vivax has not responded well to elimination efforts. These dormant parasites may reactivate and thereby cause disease and contribute to further transmission of the disease. Though it is often assumed that reactivations of dormant P vivax parasites occur at a constant rate, it has also been proposed that there is a time of increased risk of reactivation (‘temporal heterogeneity’) and there may be differences in individual’s reactivation risks (‘population heterogeneity’). We created models for constant reactivations, temporal heterogeneity, and population heterogeneity which we use to analyse data of P vivax malaria events from the Thailand-Myanmar border region and Papua New Guinea. We find strong evidence for population heterogeneity as a major determinant of reactivation patterns. Further analysis of the data suggests that spatial heterogeneity in exposure to infectious mosquito bites is a potential contributor to this heterogeneity. Thus, we find that population heterogeneity plays an important role in the overall epidemiology of P vivax recurrences. Funding: This work was funded by the Australian Research Council (ARC; https://www.arc.gov.au ) (grants DP120100064 & DP180103875 (to DSK, MPD, DC) and DP200100747 (to JAF)) and the National Health and Medical Research Council (NHMRC; https://www.nhmrc.gov.au ) of Australia (grants 1082022 (to MPD, DC), 1173528 (to DC), 1141921 (to DSK), 1080001 and 1173027 (to MPD) and?Senior Principal Research Fellowship 1135820 (to NA)). LJR is supported by NHMRC ( https://www.nhmrc.gov.au ) Career Development Fellowship (GNT1161627). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript. Data Availability: The data set from the Thailand-Myanmar border region has been previous published and made available by Taylor et al. (2019) in an online repository ( https://github.com/jwatowatson/RecurrentVivax/blob/master/RData/TimingModel/Combined_Time_Event.RData ). The data set from Papua New Guinea was published by Robinson et al. (2015) and is available in an online repository ( https://datadryad.org/stash/dataset/doi:10.5061/dryad.m1n03 ). All data and code are also publicly available at https://github.com/InfectionAnalytics/Heterogeneity_Pvivax_Relapse . Copyright: © 2022 Stadler et al. This is an open access article distributed under the terms of the Creative Commons Attribution License , which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited. Despite evidence for heterogeneity in P vivax relapses, much of the mathematical modelling of P vivax recurrences to date has assumed either constant or periodic recurrence rates [ 7 , 8 , 15 , 25 , 26 ], with some exceptions. In particular, Taylor et al. have included inter-individual heterogeneity in their time-to-event modelling by including a random effect that determines the probability of the next event being a relapse or a new infection [ 8 ]. Further, other models have included age (a proxy for immunity) as a source of heterogeneity [ 4 , 14 ], and a recent study incorporated the concept of variable risk in the form of a ‘high susceptibility’ subpopulation that has a higher risk of recurrence [ 14 ]. Here, we seek direct evidence of temporal and population heterogeneity in P vivax recurrences in two previously published studies [ 4 , 8 , 23 , 27 ]. The first study by Robinson et al. [ 4 ] contains data from a randomized placebo-controlled trial of blood- plus liver-stage drugs in Papua New Guinean children. The second data set consists of two studies by Chu et al. [ 23 , 27 ] and was made available by Taylor et al. [ 8 ]. Chu et al. conducted randomized trials in patients from the Thailand-Myanmar border region with symptomatic P vivax malaria, treated them with various therapies, and recorded their P vivax recurrences over a one-year follow-up [ 23 , 27 ]. A key difference in these studies is that the Robinson et al. study recruited individuals from communities regardless of whether they were infected, whereas the Chu et al. studies focus on treatment of individuals with symptomatic P vivax malaria at enrolment. These datasets provide an excellent opportunity for modelling of P vivax recurrences in areas with a short latency P vivax phenotype. We compare whether a model that considers (i) only a constant risk in relapse, (ii) temporarily increased relapse risk, (iii) individual variation in relapse risk, or (iv) both temporal and population heterogeneity, best explains the pattern of relapses observed in these trials. The evidence suggests that population heterogeneity in relapse risk is a primary determinant of relapse patterns. We speculate that this population heterogeneity in relapse is likely mediated by different hypnozoite burdens in individuals due to variation in individual exposure to inoculation with sporozoites. This agrees also with the observation that in the Thailand-Myanmar data, almost 70% of all P vivax infections occur in only 30% of the patients who are treated only for blood-stage infections [ 6 , 28 ]. In addition to population heterogeneity, we find evidence that temporal heterogeneity may also contribute to the overall relapse kinetics in these studies, albeit to a lesser extent. Several factors are thought to influence the timing of hypnozoite activation leading to relapse. In particular, recent infection with P falciparum or another infectious agent, and other factors may cause a temporary elevation in the risk of relapse [ 8 , 10 , 11 ]. Besides these temporal factors influencing the risk of relapse (‘temporal heterogeneity’), other factors such as transmission intensity, variation in latency due to differences in P vivax strains (e.g. temperate and tropical latency phenotypes) [ 6 , 12 ], host immunity and age are also known to influence the pattern of blood-stage infections [ 4 , 6 , 10 , 13 , 14 ]. Further, across a patient population, individuals are likely to harbor different numbers of hypnozoites due to differences in infection risk and a skewed distribution of sporozoite numbers inoculated [ 6 , 15 , 16 ], use of primaquine for radical cure and CYP2D6 polymorphisms causing treatment failures in some individuals [ 17 ], or variation in infection susceptibility [ 18 ] (which for example may be due to G6PD deficiency [ 19 ]). Given these differences, we might expect that individuals within a population may vary significantly in the number of hypnozoites they harbor, which is expected to contribute directly to their risk of hypnozoite reactivation and P vivax relapse [ 20 ]. In particular, studies of P cynomolgi in Rhesus macaques have shown that there is a correlation between the sporozoite inoculation and relapse frequency [ 21 ]. Thus, more exposed individuals within a community are expected to be more likely to have a larger hypnozoite reservoir and to have more frequent relapses [ 6 , 16 ]. Indeed, heterogeneity in malaria transmission, infections, and P vivax recurrences (i.e., P vivax infections that are either due to a new, mosquito-borne infection or a relapse) have been described previously [ 13 , 14 , 18 , 22 ]. For example, Chu et al. found that all recurrences observed after day 35 post enrolment occurred in only 12% of individuals [ 23 ]. This is consistent with studies of P falciparum infection, in which Cooper et al. found evidence that 20% of the population accounted for around 80% of transmission events [ 24 ]. Plasmodium vivax is a major cause of clinical malaria with about 4.5 million cases in 2020 [ 1 ]. An important feature of P vivax malaria is the ability of the parasite to form latent liver-stage parasites (hypnozoites), which can later activate and initiate blood stage infection in the absence of new mosquito inoculation [ 2 – 6 ]. It has been estimated that 66% to 96% of P vivax blood-stage infections are relapses caused by the activation of hypnozoites [ 2 , 4 , 7 – 9 ]. These and other studies have highlighted that hypnozoite reactivation is a major source of observed blood-stage infections and presents a major barrier to effective control and eradication of P vivax malaria [ 4 , 5 , 7 , 9 ]. Although the drug primaquine can effectively clear hypnozoites, its use has been limited in part because it can induce severe haemolysis in glucose-6-phosphate dehydrogenase (G6PD) deficient individuals [ 2 , 4 , 5 , 9 ]. Methods Data: Papua New Guinea data The Papua New Guinea data set contains the data from a randomized placebo-controlled trial of liver-stage treatment using primaquine. The data was made publicly available by Robinson et al. [4] (where the details of the study can be found). The trial was conducted between August 2009 and 20 May 2010 in 5 village clusters in the Maprik District, East Sepik Province, Papua New Guinea. Children were enrolled and treated for blood-stage infections with chloroquine or artemether-lumefantrine and with either primaquine for liver-stage infections or with a placebo. Primaquine was administered at a dose of 0.5 mg/kg/day for 20 days (10 mg/kg total dose). After initial treatment, there were fortnightly active surveillance visits for 32 weeks and passive surveillance throughout the trial. The data contains (amongst others) the time from enrollment to P vivax infection by PCR and microscopy from 504 children as well as the children’s age (between 4 and 10 years), village cluster, and their treatment (placebo or primaquine). For the model fitting, we used the time to P vivax infection by PCR data and the village cluster information. Data: Thailand-Myanmar border region data We consider the combined data from two randomized trials conducted in the Thailand-Myanmar border region, the VivaX History trial (VHX) [27] and the Best Primaquine Dose trial (BPD) [23] as published by Taylor et al. and made available online [8]. Details of the studies have been published previously [8,23,27]. In short, the VHX study was conducted between May 2010 and October 2012, it included 644 patients who were enrolled with uncomplicated P vivax malaria and randomized to either artesunate monotherapy, chloroquine monotherapy, or chloroquine with high dose primaquine. The BPD study was conducted between February 2012 and July 2014, it included 655 patients with symptomatic P vivax malaria, randomized to a primaquine (at one of two different doses) with either chloroquine or dihydroartemisinin-piperaquine (Table I in S1 Tables). Patients received either 0.5 mg/kg/day of primaquine for 14 days or 1 mg/kg/day of primaquine for 7 days, i.e., the total primaquine dose was 7mg/kg for all patients. In both studies, infections were diagnosed using a malaria smear and antimalarial treatment was allocated based on a factorial design. In the VHX study, patients were excluded from primaquine treatment if they were found to be glucose-6-phosphate dehydrogenase (G6PD) deficient and in the BPD study G6PD deficiency was an exclusion criterion. Individuals were followed up for one year from enrolment and recurring P vivax infections were treated with the same antimalarial treatment as the first P vivax infection (VHX study) or with the standard chloroquine and primaquine regimen (BPD study). The data includes patient id, episode number, antimalarial treatment, time from last event to current event, censored variable, study, and overall follow-up time. One individual was excluded from the data as the data indicates censoring at the time of enrolment. This leaves 1298 individuals (446 patients who received blood-stage treatment and 852 patients who received both primaquine and blood-stage treatment). Data was analyzed in R (version 3.6.0) [29] using the survival [30, 31] and survminer [32] packages. Models We considered four different mathematical models for the risk of relapse and recurrence in the cohorts. All models include a prophylactic effect of the antimalarial before drug washout and a constant rate of new infections. Patients are protected at enrolment due to the prophylactic effect of the antimalarial treatment, the time to drug washout is assumed to be lognormal distributed, and after drug washout, individuals become susceptible to new infection and relapses (see Fig A in S1 Methods for a model scheme). The models differ in the construction of the relapse rate. Model 1: constant relapse rate In model 1, the recurrence rate is the sum of the constant rate of new infections and the constant relapse rate. Thus, the dynamics of the fraction of susceptible individuals S(t) at time t is given by: where w(t;μ,σ) is the distribution in drug washout times across the population and is assumed to be described by the probability density function of the lognormal distribution with parameters μ and σ, r is the constant relapse rate, and n is the constant rate of new infections. Model 2: temporal heterogeneity Model 2 considers a relapse rate that is time dependent. Since our data analysis suggests that the relapse rate decreases after an initial peak (Fig 1), we chose a decreasing relapse rate such that the prophylactic effect of the antimalarial treatment in combination with the decreasing relapse rate can capture the observed change in the relapse rate (Fig 1). We assume that there is an initial relapse rate (I) and that the relapse rate decreases exponentially over time (with rate d). Thus, the relapse rate is given by r(t) = Ie−dt and the dynamics for the fraction of susceptible individuals at time S(t) is given by: where w(t;μ,σ) is the probability density function of the lognormal distribution with parameters μ and σ, r(t) = Ie−dt is the time-dependent relapse rate, and n is the constant rate of new infections. PPT PowerPoint slide PNG larger image TIFF original image Download: Fig 1. Time to first recurrence and weekly incidence per patients at risk in the Thailand-Myanmar and Papua New Guinea data. (A, B) Time from enrolment to the first recurrence for patients who received blood-stage treatment only (red) and patients who received primaquine and blood-stage treatment (blue). In the Papua New Guinea data, all patients (n = 504) received blood-stage treatment and either a placebo (red) or primaquine (blue). For the time to first recurrence in the Thailand-Myanmar data (n = 1298) by treatment and study see Fig H in S1 Figs (C, D) Weekly incidence per patients at risk for patients treated with primaquine and blood-stage treatment (blue dots) and blood-stage treatment only (red dots). The weekly incidence per patients at risk is the number of patients that had a recurrence within the current week divided by the number of patients who were at risk (i.e., the patients who have not yet had a recurrence) at the beginning of the week. The curves are fitted to the weekly incidence per patients at risk (using splines). https://doi.org/10.1371/journal.pntd.0010990.g001 Model 3: population heterogeneity In model 3, the relapse rate is constant but drawn from a lognormal distribution to model population heterogeneity. To simplify the numerical computations, we group the population in k different risk groups of equal size. The relapse rate for each risk group is the median relapse rate for this group that is computed using the percentiles of the lognormal distribution of relapse rates (see Fig B in S1 Methods). Thus, model 3 is given by: where S i (t) is the fraction of susceptible individuals that are in risk group i (i∈{1,2,…,k}) at time t, k is the number of relapse risk groups, r i is the median relapse rate of group i, and n is the constant mosquito-borne infection rate. The overall fraction of susceptible individuals at time t, S(t), is the sum of the fractions of all susceptible individuals in the different risk groups: Model 4: temporal and population heterogeneity Finally, we considered a model that considers both temporal and population heterogeneity and is a combination and extension of models 2 and 3. As for model 3, we group the population in k different relapse risk groups of equal size and as for model 2 this relapse risk decreases over time, i.e., r i (t) = I i e−dt where I i is the initial relapse risk for relapse risk group i. The initial relapse risk (I i ) is the median relapse rate for this group that is computed using the percentiles of the lognormal distribution of relapse rates (see Fig B in S1 Methods). Model 4 is given by: where S i (t) is the fraction of susceptible individuals that are in risk group i (i∈{1,2,…,k}) at time t, k is the number of relapse risk groups, r i (t) = I i e−dt is the time-dependent relapse rate of group i, and n is the constant mosquito-borne infection rate. As for model 3, the overall fraction of susceptible individuals at time t, S(t), is the sum of the fractions of susceptible individuals in the different risk groups. To fit the models not only to the first recurrence after enrolment but to the first and second recurrence, we extended the models to take two recurrences into account by adding additional compartments to the model (see Fig C in S1 Methods). After a susceptible individual has a P vivax recurrence, the patient is again protected with the same drug washout time distribution as at enrolment because patients are treated with the same antimalarials at each recurrence. As for the first recurrence, patients become susceptible to mosquito-borne infections and relapses after drug washout. For all models and all simulations shown, we assume that the drug washout time depends on the antimalarial treatment, the rate of new infections differs between the two different studies in the Thailand-Myanmar data as incidence rates changed between the time of the two studies (8) (Table A and Fig A in S1 Model comparisons), and we assumed that patients treated with primaquine have no relapses, i.e., we assume they are subject only to new infections (as previously; see Model fitting below and Supplement). Note that since only the Thailand-Myanmar data contains multiple recurrence times, we used only the Thailand-Myanmar data for fitting to the first and second recurrence time and simulating multiple recurrences under the temporal and population heterogeneity models. The Papua New Guinea data on the other hand contains spatial information (villages) which we used to study the association between new infections and relapses. Model fitting to recurrence times Each model was implemented as an ordinary differential equation (ODE) or a system of ODEs and solved numerically using the function ode15s in MATLAB [33]. We constructed a likelihood function using this numerical solution to fit our models to the recurrence time data (including censored time intervals) as described below and obtain maximum likelihood estimates for the model parameters. From the numerical solution of the model equations, we obtained the fraction of uninfected individuals at time t, U(t), as the sum of the fraction of individuals who are still protected and the fraction of susceptible individuals: where w(t;μ,σ) is the probability density function of the lognormal distribution of drug washout times with parameters μ and σ. We interpret U(t) as the probability of remaining uninfected until time t. The probability of having an infection at the follow-up visit on day t (G(t)) is then the probability of being uninfected on the previous follow-up visit but infected on day t, i.e., where t—Δ is the time of the last follow-up visit before day t. Using U(t), the probability of being uninfected until time t, and G(t), the probability of having an infection on day t, we used the following loglikelihood function to fit the models to the first recurrence times: where p are the parameters of the model (see S1 Methods for the parameters of the different models), D is the data, i.e., the recurrence times, N 0 and N 1 are the numbers of individuals with zero and at least one recurrence, respectively, U(t) is the probability of being uninfected until time t, t i is the time of the first recurrence of individual i, G(t) is the probability of having an infection at the follow-up visit on day t, and T i is the overall follow-up time of individual i (for censored time intervals). For fitting the data to the first and second recurrence time simultaneously (in the Thailand-Myanmar data), we constructed a different loglikelihood function. We use both the first and the second recurrence information (as well as censoring times) from the data in our likelihood function and also incorporate the different relapse risk groups in models 3 and 4 to allow for heterogeneity in relapse rates across the population. The loglikelihood function for the parameter set p and the data D is given by: where N 0 , N 1 , and N 2 are the numbers of individuals with 0, 1, or at least 2 recurrences, respectively, r is the number of risk groups, U j (t) is the probability of remaining uninfected to time t for individuals in risk group j, T i is the follow-up time of individual i (i.e., the time of censoring), G j (t) is the probability of having an infection on the follow-up visit on day t for individuals in risk group j, and t i,1 , t i,2 are the times from the start of the study to the first recurrence and from the first to the second recurrence of individual i, respectively. Note this loglikelihood function can be simplified for models 1 and 2 as there is only one risk group in these models (see S1 Methods for the derivation and more details on the loglikelihood function). The model was fit by minimizing the negative loglikelihood function with the function fmincon in MATLAB [33]. We did 100 minimizations of the negative loglikelihood function with random initial parameter values to find the Maximum Likelihood Estimates (MLEs) for the parameters. Additionally, when fitting the extended model to the first and second recurrence data, we used the parameters of the best fit to the first recurrence times as initial parameter values. Confidence intervals were computed using bootstrapping and the percentile method (see S1 Methods). We fit models 1 to 3 using one mosquito-borne infection rate for the Papua New Guinea data and two different infection rates for the two studies in the Thailand-Myanmar data (Table A and Fig A in S1 Model comparisons). We also assessed the importance of assumptions regarding the follow-up schemes (Table B, Figs B, C, and D in S1 Model comparisons), different numbers of relapse risk groups in model 3 (Table C and Fig E in S1 Model comparisons), and the relapse risk distribution in the population heterogeneity model (Table D, Figs F and G in S1 Model comparisons) in the Thailand-Myanmar data. For the final model comparisons and simulations both for fitting to the first recurrence and for fitting to the first and second recurrence simultaneously, we chose two different infection rates, daily follow-ups, and 10 relapse risk groups for all models. [END] --- [1] Url: https://journals.plos.org/plosntds/article?id=10.1371/journal.pntd.0010990 Published and (C) by PLOS One Content appears here under this condition or license: Creative Commons - Attribution BY 4.0. via Magical.Fish Gopher News Feeds: gopher://magical.fish/1/feeds/news/plosone/