# Probabilistic classification of late treatment failure in uncomplicated malaria

Somya Mehra, Aimee Taylor, Mallika Imwong, Nicholas J White, James A Watson

► To cite this version:

Somya Mehra, Aimee Taylor, Mallika Imwong, Nicholas J White, James A Watson. Probabilistic classification of late treatment failure in uncomplicated malaria. 2025. pasteur-05272916

HAL Id: pasteur-05272916

https://pasteur.hal.science/pasteur-05272916v1

Preprint submitted on 22 Sep 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Probabilistic classification of late treatment failure in uncomplicated malaria

Somya Mehra¹,²*, Aimee R Taylor³, Mallika Imwong¹,⁴, Nicholas J White¹,², James A Watson²,⁵

1. Mahidol Oxford Tropical Medicine Research Unit, Mahidol University, Bangkok, Thailand

2. Nuffield Department of Medicine, University of Oxford, Oxford, United Kingdom

3. Infectious Disease Epidemiology and Analytics G5 Unit, Institut Pasteur, Université Paris Cité, Paris, France

4. Department of Molecular Tropical Medicine and Genetics, Faculty of Tropical Medicine, Mahidol University, Bangkok, Thailand

5. Infectious Diseases Data Observatory, Big Data Institute, Oxford, United Kingdom

*Corresponding author: somya@tropmedres.ac

## Abstract

Distinguishing treatment failure (recrudescence) from reinfection in uncomplicated falciparum malaria is essential for characterising antimalarial treatment efficacy in malaria endemic areas. Classification of recrudescence versus reinfection is usually based on a comparison of parasite allelic calls derived from PCR amplification and electrophoresis of individual polymorphic markers in the acute and recurrent isolates. Match-counting methods (e.g. 3/3 or 2/3 matching alleles) have usually been applied, but these do not account for multiple comparisons per-marker when infections are polyclonal. We show that when infections are polyclonal, as is common in high transmission settings, currently used match-counting and model-based methods have unacceptably high false-discovery rates leading to overestimation of treatment failure. We develop the software *PfRecur* which provides analytical Bayesian posterior probabilities of treatment failure in recurrent falciparum malaria. We use data from a recent study in Angola to demonstrate the potential utility of our model in resolving complex polyclonal *P. falciparum* infections, thereby providing more accurate estimation of treatment failure rates.
Introduction

*Plasmodium falciparum* is estimated to cause around 250 million symptomatic malaria cases each year, the majority of which are in young children in sub-Saharan Africa [1]. The primary therapeutic goal of antimalarial treatment for uncomplicated malaria is to cure the infection. Effective treatment of uncomplicated malaria currently relies on artemisinin-based combination therapies (ACTs). The ACTs combine a rapidly acting but rapidly eliminated artemisinin derivative with a less active and slowly eliminated partner drug (for example, the combination of artemether and lumefantrine, AL). Artemisinin resistance is now widespread in Southeast Asia [2, 3] and recently has emerged independently in East Africa [4, 5]. Resistance to artemisinin results in an increased likelihood that the ACT fails to clear the blood stage infection and recrudesces. This amplifies the selective advantage of artemisinin resistant parasites and helps select for partner drug resistance. Characterising ACT failure rates is essential to guide policies and practices. This assessment is routinely performed through therapeutic efficacy surveillance studies. These are usually single arm observational studies which enrol symptomatic uncomplicated malaria patients and follow them for one to two months after treatment. The proportion of patients with a recrudescent infection (i.e., treatment failure) is then estimated [6, 7]. Treatment failure rates should not exceed 10% [1]. In areas of high transmission where the main burden of disease is in young children and where reinfection within one month is very common, estimating the antimalarial treatment failure rate relies on distinguishing recurrent bloodstream infections resulting from incomplete clearance of the incident infection (recrudescence) from new infections (reinfection) resulting from new mosquito bites [8, 9]. This relies on molecular correction methods which compare parasite genotypes between the baseline infection with those of the recurrent infection.

Classification of a recurrent infection as either a reinfection or a recrudescence is primarily done by PCR genotyping. Length polymorphisms within genes or microsatellite markers are compared in paired acute and recurrent isolates [9]. Classification is then based on the observed alleles. Most studies have used match-counting approaches, such as those recommended by the World Health Organisation (WHO) [9]. Match-counting algorithms usually take the presence of one or more shared alleles per marker in the paired isolates as evidence of recrudescence, and can either be strict (requiring matches at all markers) or relaxed (requiring matches at a subset of markers). However, this does not account sufficiently for multiple comparisons when there are multiclonal isolates (infections caused by multiple parasite clones), the allele specific sensitivity of the genotyping method [10], the relative frequency of the different alleles, and uncertainty in the allocation of alleles across individual haploid parasite clones within each isolate [8, 11]. Concerns have been raised around the reliability of match-counting algorithms. These highlight technical constraints associated with molecular genotyping [12–14] and have been evaluated in simulation models [15, 16]. Deep sequencing of longer amplicons from SNP-rich genomic regions (AmpSeq) addresses some of these concerns, but is not yet widely available. Statistical model-based approaches can address classification uncertainty correctly and take into account background allele frequencies, adjusting for multiple comparisons. The US CDC has developed a Bayesian classifier in which allelic states across clones are explicitly estimated using a Gibbs sampler and the likelihood of recrudescence is averaged over
each pairwise comparison of clones in the baseline and recurrent isolate [8]. An informal consultation convened by WHO in 2021 recommended study of “the impact of Bayesian analysis on recrudescence rates in areas of high transmission, given the trend in increased recrudescence rates using this method, to determine whether this is an artefact of the method, whether there is some reason for higher failure rates in these areas (e.g., a high MOI may be more challenging for antimalarial drugs to clear), or whether there is emergence of true antimalarial resistance that needs rigorous confirmation” [9].

In this work, we use microsatellite data provided from 70 paired *P. falciparum* recurrent infections from a study conducted in Angola in 2021 [17] to evaluate the CDC Bayesian classifier and explore the impact of model misspecification on recurrence classification. We formulate an alternative probabilistic Bayesian classifier (implemented in the R software *PfRecur*) for genotyping data from paired baseline and recurrent falciparum infections, allowing recurrent infections to be mixtures of recrudescent and newly-inoculated clones [9, 18]. We show that this formulation makes *PfRecur* robust to misspecification in the likelihood and solves combinatorial problems for multiclonal isolates, both for the pairwise comparisons between baseline and recurrent isolates [19], and the computation of sample allele frequencies [20–24]. Because the posterior is analytically tractable, *PfRecur* is fit to data analytically, precluding the need for an optimisation algorithm or numerical sampler. We compare the empirical performance of simple match-counting algorithms versus the CDC Bayesian classifier and our novel approach *PfRecur*. We show that in high transmission areas, where infections are often polyclonal, the CDC Bayesian classifier and match-counting approaches (as recommended by [16, 17, 25]) both give unacceptably high false positive rates. This can result in overestimation of treatment failure rates and thus overestimation of antimalarial drug resistance, potentially prompting expensive changes in treatment policy [26, 27].

# Results

## PfRecur: a probabilistic model of recrudescence

We developed a novel probabilistic model for recrudescence versus reinfection based on observed alleles from each marker in each baseline and paired recurrent infection; for simplicity we refer to this model as *PfRecur*. Rather than treating reinfection and recrudescence as mutually exclusive categories, each recurrent infection is modelled as a *mixture* of newly-inoculated (from reinfection) and recrudescent clones. This construction is motivated by a concern for robustness against model misspecification. Under the simplifying assumptions of marker-wise and clone-wise independence (within isolates and between non-recrudescent clones across isolates), we derive an analytic per-participant likelihood that averages over all allelic configurations that are compatible with the isolate MOI and the set of observed genotypes in each isolate. We accommodate undetected clones in baseline/recurrent isolates for the participant of interest only, whereby each clone is modelled to be genotyped (i.e, detected) with probability $\omega$ at each marker; this yields a marker-wise truncated binomial model for the number of clones that have contributed to the observed set of genotypes at each marker. Site specific allele frequencies for newly-inoculated clones are derived under a multinomial-Dirichlet model from the baseline isolates, assuming exchangeability among
study patients (excluding the patient of interest). These site specific allele frequencies are also used to impute allelic states at each marker for undetected clones in the baseline isolate for each patient. Allele frequencies in the recrudescent clones are based on the paired patient baseline clones following imputation. We adjust for genotyping error in baseline isolates, relative to the recurrent isolate of interest, through a non-parametric marker-agnostic model $\delta_{\ell}$ governing the probability that each allele called in a baseline isolate matches an allele called in the recurrent isolate; in this study, we consider a normalised geometric error model (with respect to allele repeat lengths) with error probability $\varepsilon$, adapted to length polymorphic markers. $\omega, \varepsilon$ are user-specified parameters in the model with default values of 0.9 and 0.05, respectively.

Under a Bayesian framework, with a symmetric beta binomial prior for the number of newly-inoculated vs recrudescent clones within a recurrent isolate, we derive an analytical posterior distribution for the number of recrudescent clones within each recurrent isolate. We consider two posterior summary metrics: the posterior probability of there being at least one recrudescent clone in the recurrent isolate (denoted M1); and the expected proportion of recrudescent clones in the recurrent isolate (M2). These metrics are identical for recurrent isolates with a single clone (i.e., MOI=1).

## Simulation study

We conducted a simulation study to evaluate the ability of our probabilistic classifier *PfRecur* to resolve mixtures of newly-inoculated and recrudescent clones. Simulated data (Appendix B) were well-specified relative to our classification model apart from three features: non-meiotic siblings (i.e. genetically related parasites derived from independent crosses of the same parental pair) could be present within isolates (thereby violating the assumption of clone-wise independence within isolates); the observed MOI of each isolate (i.e., the maximum observed cardinality across simulated markers) could be lower than the true MOI (because of undetected clones or allelic overlap between clones); and the detection of clones in baseline isolates (used to derive allele frequencies for newly-inoculated clones) was incomplete. We considered an idealised setting, where both genotyping error probabilities $\varepsilon$ and the per-clone marker-wise detection probability $\omega$ under which data were generated were known.

Using metric M2, our model *PfRecur* was able to recapitulate the proportion of recrudescent clones within simulated recurrences, albeit with decreasing precision for recurrences with high MOIs (Figure 1A). Metric M1 (i.e., the posterior probability of at least one recrudescent clone) was sensitive to simulated recurrences with $\ge 1$ recrudescent clone but was liable to call false positive recrudescences using a threshold of 0.5 (dashed blue line) at elevated MOIs (Figure 1B). Classification using a threshold of $1/m$ for recurrent isolates with elevated MOI $m$ may allow recrudescences to be called with greater specificity.

## Therapeutic efficacy study in Angola

We used open access data from a therapeutic efficacy study conducted in Angola in 2021 [17] to compare *PfRecur* with match-counting algorithms and the CDC model. The Angolan study evaluated artemether-lumefantrine (AL), artesunate-amodiaquine (ASAQ), dihydroartemisinin-piperaquine (DP),
<figure><img src="image_2.png" /><figcaption>**Figure 1:** Application of our probabilistic classifier *PfRecur* to simulated mixtures of newly-inoculated and recrudescent clones. Panel A shows the posterior probability of at least one recrudescent clone (M1) (annotated with classification sensitivity and specificity using a threshold of 0.5 to call recrudescence), while panel B shows the posterior proportion of recrudescent clones (M2), aggregated over 28773 simulated recurrences. Truth values are shown with vertical lines.</figcaption></figure>
<table><caption>**Table 1:** Summary of a therapeutic efficacy study conducted in Angola in 2021 by Dimbu *et al.* [17] (adapted from Tables 1 and 3 of [17]). ASPY=artesunate-pyronaridine, DP=dihydroartemisinin-piperaquine, AL=arthemether-lumenfantrine, ASAQ=artesunate-amodiaquine. *Genotypes are missing for one recurrent isolate from the ASPY arm in Benguela. †Based on the maximum observed cardinality across 7 neutral microsatellite markers for the day 0 isolate.</caption><thead><tr><th rowspan="2"></th><th colspan="2">Benguela</th><th colspan="2">Zaire</th><th colspan="2">Lunda Sul</th></tr><tr><th>ASPY</th><th>DP</th><th>AL</th><th>ASAQ</th><th>AL</th><th>ASAQ</th></tr></thead><tbody><tr><td>Study period</td><td>Mar-May 2021</td><td>May-Jul 2021</td><td>Feb-Mar 2021</td><td>Mar-Jun 2021</td><td>Feb-Apr 2021</td><td>Mar-Jun 2021</td></tr><tr><td>No. patients enrolled</td><td>104</td><td>105</td><td>104</td><td>105</td><td>104</td><td>100</td></tr><tr><td>No. patients finished follow-up</td><td>100</td><td>104</td><td>98</td><td>97</td><td>100</td><td>91</td></tr><tr><td>Median age, years (range)</td><td>6.9 (0.5-12)</td><td>7 (0.6-12)</td><td>2.5 (0.5-5)</td><td>2.5 (0.7-5)</td><td>2.6 (0.5-5)</td><td>2.7 (0.6-5)</td></tr><tr><td>% patients female</td><td>49%</td><td>47%</td><td>45%</td><td>56%</td><td>49%</td><td>46%</td></tr><tr><td>No. early treatment failures (day&lt;7)</td><td>0</td><td>0</td><td>1</td><td>3</td><td>0</td><td>0</td></tr><tr><td>No. late treatment failures (day≥7)</td><td>14* (3*)</td><td>6 (4)</td><td>22 (12)</td><td>16 (7)</td><td>13 (10)</td><td>0 (0)</td></tr><tr><td>(with baseline MOI&gt;1<sup>†</sup>)</td><td></td><td></td><td></td><td></td><td></td><td></td></tr></tbody></table>

and artesunate-pyronaridine (ASPY) across 3 study sites in 622 patients with uncomplicated malaria [17]. Recurrent parasitaemia was detected in 71 patients between 7 and 42 days of follow-up (Table 1). Paired baseline and recurrent *P. falciparum* isolates (70 pairs) were genotyped at 7 neutral microsatellite markers (*M2490*, *M313*, *M383*, *PfPK2*, *POLYA*, *TA1*, *TA109*). An additional 14 unpaired baseline isolates per study site were also genotyped to increase the precision of the estimation of population allele frequencies, stratified by study site. Transmission intensities varied across the three sites, with Zaire and Lunda Sul classified as high and moderate intensity transmission, respectively, and Benguela as low transmission intensity. In Zaire and Lunda Sul combined, there were 36/198 (18%) recurrent infections in the AL treated patients within 28 days of follow-up; and 19/188 (10%) recurrent infections in the artesunate-amodiaquine (ASAQ) treated patients within 28 days of follow-up. In Benguela, dihydroartemisinin-piperaquine (DP) and artesunate-pyronaridine (ASPY) were compared, with 14/100 and 6/104 recurrent infections within 42 days of follow-up, respectively. Under the CDC model, this translated into estimated day-28 PCR-corrected efficacies of 88% (95% CI 82-95) and 94.4% (95% CI 90-99) for AL; and 91% (95% CI 85-97) and 100% for ASAQ; and day-42 PCR-corrected efficacies of 99.6% (95% CI 99-100) for DP and 98.3% (95% CI 96-100) for ASPY. Of particular concern was the low efficacy estimate of AL in the Zaire site. AL is the most frequently used ACT in the world [28].
<figure><img src="image_3.png" /></figure>

*Expected value of 1-Nei's gene identity metric, assuming M equifrequent clones (where M is the maximum cardinality across all loci) exhibit A distinct alleles (where A is the locus cardinality)

**Figure 2:** (A) Parasitemia stratified by MOI (i.e., the maximum cardinality across 7 neutral microsatellite markers) for 181 *P. falciparum* isolates reported by Dimbu *et al.* [17]. (B) The within-host diversity of each microsatellite marker as a function of parasitemia.

## Artefacts at low parasite densities

Visual inspection of the genotypic data from [17] showed that some of the markers were much more likely to be polyallelic than others. Exploratory analysis demonstrated that the parasite density (log₁₀ parasitaemia) was strongly predictive of the observed MOI (maximum cardinality across the 7 markers; Figure 2A). Low densities (<1000 parasites per µL) were associated with MOIs of 2 or more. When assessing individual markers, it was apparent that the microsatellite *TA109* is problematic particularly at very low parasite densities (Figure 2B), yielding elevated apparent MOIs for recurrent parasitaemias below 1000/µL. This suggests that *TA109* multiplicity in low parasite density samples may be artefactual and result from methodological or laboratory artefacts (e.g. non-specific peaks). We conducted all analyses with and without *TA109* to test for robustness of the methodology relative to inclusion of an unreliable polymorphic marker.
## Estimation of treatment failure for isolates with baseline MOI>1

We estimated the probability of recrudescence for each patient with a recurrent infection in the Angola study, using the CDC model and *PfRecur*. We used metric M1 (probability of at least one recrudescent clone) for conceptual consistency with the CDC model [8]. For paired infections with a baseline MOI of 1, the probabilistic models do not differ substantially. Both models correct for chance allelic matches, which is an advantage over a simple match-counting algorithm. However, when the MOI of the baseline exceeds 1, issues of multiple comparisons become important. Figure 3 shows the model estimates for the 36 paired infections in the Angolan study with baseline MOI> 1, ordered by the posterior probability of recrudescence under the CDC model. For 13/36 (one third) recurrent infections, there was a non-negligible difference in the model estimates, with the CDC model consistently estimating higher recrudescence probabilities. Discrepancies were largely apparent for recurrences with intermediate posterior probabilities under the CDC model, and many of these differences persisted after exclusion of the problematic marker *TA109*. Visual inspection of the isolates where the models differ substantially suggest that the CDC model is over-calling recrudescence (Supplementary Figure 11).

## Estimation of false positive rates for recrudescence classification

We used a permutation method to estimate false positive rates for recrudescence classification. We generated 500 artificial ‘not-recrudescence’ datasets by shuffling the participant identifiers of the baseline isolates in [17], stratified by study site thus preserving population structure. Figure 4 shows the estimated false positive rates, calculated by averaging metric M1 of *PfRecur* and the posterior probability of recrudescence under the CDC model across recurrent isolates. Across the permuted datasets, the CDC model had median false discovery rates of 8.7% (95% confidence interval [CI] 1.2–19.1) in Benguela, 6.5% (95% CI 0.7–18.6) in Lunda Sul and 5.0% (95% CI 1.3–10.2) in Zaire, driven largely by permuted pairs with intermediate posterior probabilities of recrudescence (Supplementary Figure 10). In comparison, *PfRecur* had substantially lower median false discovery rates of 0.6% (95% CI 0.02–7.1) in Benguela, 0.1% (95% CI 0–2.8) in Lunda Sul and 0.4% (0.02–3.9) in Zaire, with false discovery rates remaining comparatively low even as the per-clone marker-wise probability of detection was relaxed from the default value $\omega = 0.9$ to $\omega = 0.75$ (assuming more clones evade detection tends to augment the posterior probability of recrudescence) (Supplementary Figure 9). The ‘≥ 4/7’ match-counting rule recommended by [16, 17, 25] yielded higher false discovery rates than both *PfRecur* and the CDC model.

We note that the false discovery rates are dependent in part on marker diversity: in settings with limited diversity, we would expect *PfRecur* to return the prior distribution over newly-inoculated vs recrudescent clones, and the CDC model to return a 0.5 posterior probability of recrudescence, yielding an elevated apparent false positive recrudescence rate using this method.
<figure><img src="image_4.png" /><figcaption>**Figure 3:** Summary of posterior estimates for *P. falciparum* recurrences reported by Dimbu *et al.* [17] with MOI>1 for the baseline isolate. We compare the probability of at least one recrudescent clone under *PfRecur* (with genotyping error probability ε = 0.05 under the normalised geometric model (29) and marker-wise per-clone detection probability ω = 0.9, applied to clones in the pair of baseline/recurrent isolates for the participant of interest) against the posterior probability of recrudescence based on the CDC model [8].</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>**Figure 4:** False positive recrudescence rates under *PfRecur* (averaging metric M1) vs the CDC model [8] (averaging the posterior probability of recrudescence) vs a match-counting approach (i.e., treating the presence of one or more shared alleles at 4 or more, or 5 or more, markers as evidence of recrudescence [16, 17, 25]; restricted to permuted pairs with ≥ 1 allele call at each of the 7 markers). Points indicate medians, bold error bars indicate the interquartile range and gray error bars indicate 95% confidence intervals over 500 permuted artificial 'not-recrudescence' datasets, each comprising 19 recurrent isolates from Benguela; 13 recurrent isolates from Lunda Sul and 38 recurrent isolates from Zaire.</figcaption></figure>

Discussion

Accurate characterisation of ACT failure rates in uncomplicated malaria is essential to guide policies and practices, especially now in Africa where artemisinin resistant *P. falciparum* is spreading and consequently ACTs are under threat [29]. Genotyping of polymorphic alleles has allowed evaluation of antimalarial therapeutic efficacy in malaria endemic areas where study participants are exposed to new infectious bites during follow-up. But in high transmission settings differentiating between reinfection and recrudescence remains a difficult problem. Many malaria infections are polyclonal, posing combinatorial problems since allele counts are not directly observable. The different PCR genotyping methods also have different sensitivities [10]. Limitations of match-counting methods – which do not adjust for multiple comparisons across polyclonal isolates, relative allele frequencies [8, 11] or the imperfect detectability of parasite clones [13, 14] – have been well-established in the literature [14–16]. In addition to having good laboratory techniques for accurate genotype calling, a robust statistical methodology is needed to assess probabilistically whether the recurrent infection is compatible with recrudescence (treatment failure). We show that when infections are polyclonal, match-counting methods and the methodology proposed by the CDC [8] may result in unacceptably high false positive recrudescence rates. In contrast, *PfRecur* has higher specificity whilst retaining good accuracy at identifying recrudescence. When applied to microsatellite data from a recent therapeutic efficacy assessment in Angola [17], our classification model outputs substantially different estimates for one third of the recurrent infections with baseline MOI> 1 [8]; however this does not greatly impact efficacy estimates across study arms (Supplementary Table 4). The systematic overestimation of recrudescence rates is likely to be greater in settings with greater parasite diversity, higher polyclonality, and more frequent reinfection.
The Bayesian model proposed by the CDC has three potential issues. Firstly, a per-clone unobserved allele “penalty” — a multiplicative factor applied to the likelihood when the imputed allele for a clone lies outside the observed set of alleles for a given isolate — is estimated simultaneously during classification (a relatively lax penalty in the range 0.25 to 0.45 was estimated for each study site in [17]). We suggest that this construction increases the estimated probability of recrudescence. This is because the unobserved alleles have little effect on the probability of reinfection (they only affect allele frequency estimates), but tend to augment the chance of allelic matches between the paired isolates and consequently the probability of classifying as a recrudescence. This bias is particularly pronounced for isolates with no genotyping data at a subset of markers, and isolates with unbalanced cardinality across loci. Markers with missing genotypes are included in the estimation procedure, and the per-clone unobserved allele penalty is applied irrespective of whether or not imputed alleles match those in the paired isolate. Isolates with unbalanced cardinality across loci result in a greater chance of introducing unobserved alleles, a feature that may be problematic when MOIs are inflated due to artificially high multiplicities at one marker (e.g. TA109 in the dataset analysed here). Secondly, estimation of the genotyping error appears to be unstable. In the dataset analysed, site specific estimates of genotyping error varied from 0.025 in Lunda Sul to 0.072 in Benguela (i.e. a three fold difference). However, all genotyping was presumably done by the same central laboratory. The difference in estimates suggests possible identifiability issues in the model. As demonstrated by the marker TA109, a major determinant of genotyping error is likely to be the parasite density in the filter paper. Thirdly, the reliability of the output model probabilities depend on convergence of the algorithm. To the best of our knowledge, there has been no in depth study of convergence of this algorithm, and the default parameters suggested (1000 iterations) appear insufficient. In contrast, in our proposed model *PfRecur*, we average analytically over compatible allelic configurations within each isolate to derive an analytically-tractable (discrete) posterior distribution for the number of newly-inoculated vs recrudescent clones within each recurrent isolate that can be evaluated directly, allowing for classification of recurrent isolates with multiplicity of infection up to 9 in the order of seconds.

A limitation of *PfRecur* in low transmission settings is the assumed lack of relatedness structure within or between isolates (the CDC model makes the same assumption). The classifier of Taylor *et al.* [30] for *P. vivax* recurrence explicitly accommodates sibling relationships within isolates (and between isolates, but only in support of relapse), and in addition, includes an optional “fudge-factor” for low-level background relatedness; but the transitive property of relatedness introduces substantial complexity. However, it is shown in [31] that sample allele frequencies, which are plugged into many models of malaria parasite genetic data, can be considered to encode implicitly average relatedness marker-wise.

In conclusion we present an analytically-tractable probabilistic classifier *PfRecur* for estimating recrudescence based on multi-allelic calls in recurrent falciparum malaria. In high transmission settings it is more accurate than current methods of analysis and therefore less likely to overestimate recrudescence rates.
# Methods

## Data from Dimbu et al. [17]

We re-analysed data from a six-arm therapeutic efficacy study conducted between February and July in 2021 across 3 provinces in Angola [17]. Outpatients presenting to urban clinics situated in provincial capitals were screened for inclusion in the study. Enrolment criteria included uncomplicated *P. falciparum* monoinfection and either a history of fever or an axillary temperature reading ≥ 37.5°C. To avoid confounding due to transmission intensities, some inclusion criteria varied by province. In Benguela (low to moderate transmission), children between 6 and 143 months of age with 1,000–100,000 parasites/µL at baseline were enrolled. In Lunda Sul (moderate to high transmission) and Zaire (moderate to very high transmission), children were enrolled with a narrower age range of 6 to 59 months, and higher baseline parasitemias of 2,000–200,000 parasites/µL. Parasitemia was quantified using microscopy. 622 patients were enrolled across the study arms.

Study participants were treated with 3-day regimens of either artemether-lumefantrine or artesunate-amodiaquine (Lunda Sul and Zaire), and dihydroartemisinin-piperaquine or artesunate-pyronaridine (Benguela). Dosing was determined by weight bands in accordance with manufacturer guidelines. Antimalarial treatment was largely supervised. Follow-up, entailing clinical examination and slide microscopy, occurred on days 1, 2, 3, 7, 14, 21 and 28 in addition to days 35 and 42 for patients treated with dihydroartemisinin-piperaquine and artesunate-pyronaridine. Recurrent infections, characterised by microscopy-detected asexual *P. falciparum* parasitemia from day 7 to the end of follow-up, were identified in 71 patients.

Genotyping was performed for 70 pairs of baseline and recurrent isolates, and an additional 42 baseline isolates. DNA was extracted from dried blood spots and *P. falciparum* diagnosis was confirmed by PCR. For a panel of 7 neutral microsatellite markers (*M2490*, *M313*, *M383*, *PfPK2*, *POLYA*, *TA1*, *TA109*), fragment lengths were then assessed using capillary electrophoresis. Classification of reinfection vs recrudescence was performed using the Bayesian CDC model [8], our novel probabilistic approach, *PfRecur*, and simple match-counting algorithms (at least 4/7 [16, 17, 25] or 5/7 matches, where a per-marker observation is a match if some alleles are the same at enrolment and recurrence), stratified by study site.

## Marker cardinality as a function of parasitemia in [17]

We explored the relationship between MOI (defined as the maximum cardinality across the genotyped panel of 7 neutral microsatellite panel) and parasitemia for 181 isolates with genotyping data available, using the Jonckheere-Terpstra test [32], with *p*-values derived under a normal approximation using the R function *PMCMRplus::jonckheereTest* [33]. We also examined the within-host diversity of each marker, as a function of parasitemia. We defined the within-host diversity of each marker for each isolate as the complement of Nei's gene identity metric [34], taking a uniform distribution over compatible allelic configurations at that locus (i.e., all possible ways of allocating alleles observed at that locus to the
number of clones given by the MOI). Nei's gene identity metric is formulated as the sum of squares of allele frequencies; the complement can be interpreted as the probability of differing alleles when a pair of clones is sampled (with replacement) from the isolate (see Appendix C.1).

## PfRecur

PfRecur classifies a recurrent infection as either a recrudescence or a reinfection based on either the posterior probability of there being at least one recrudescent clone (M1); or the posterior expected proportion of recrudescent clones (M2) (Figure 5). A complete description of the framework is provided in Appendix A; below, we provide a brief outline, using the terms marker and locus interchangeably. Classification of recurrent infections is done for each pair in a group of patients, where the grouping is determined by the user (usually by study or study site, depending on the context). The model is not fully Bayesian in that there is no way of specifying multilevel groupings.

<figure><img src="image_6.png" /><figcaption>**Figure 5:** *PfRecur* is designed to classify a recurrence as either a recrudescence or a reinfection using a posterior summary of *M*<sub>*C*</sub>, the number of recrudescent clones in the *r*th isolate.</figcaption></figure>
Box 1: *PfRecur* model assumptions

* **Unlinked loci**: valid if loci lie on different chromosomes (because chromosomes assort independently in meiosis) or inter-locus distances are large (because then loci are more likely separated by one or more recombination break points in meiosis).

* **Clones within isolates are independent as are non-recrudescent clones between isolates**: violated by the presence of sibling parasites within isolates [30]; additionally violated if the average population-level relatedness is high (although baseline isolates 1,...,n partially encode average relatedness [31]) and/or if the population is structured (e.g., by geographic barriers in study site, or by household effects among study participants).

* **Uniform distribution of allelic states for each isolate** (i.e., any configuration of alleles compatible with the observed MOI and the set of genotypes observed at a given locus is modelled to be equally likely): enforced in the absence of quantitative data on relative allelic abundance.

* **Reinfection and recrudescence are not mutually-exclusive**: the recurrent isolate *r* comprises a mixture of recrudescent clones drawn from the paired baseline isolate *b*, and newly-inoculated clones drawn from the contemporaneous population at large, which is approximated by the remaining baseline isolates 1,...,n.

* **Non-parametric genotyping error**: we accommodate non-parametric genotyping error in baseline isolates *relative* to the recurrent isolate, specified by the probability that each allele called in a baseline isolate matches an allele called in the recurrent isolate *r*.

* **No undetected clones in the baseline isolates 1,...,n**: we assume that the consequences of ungenotyped clones (if any) ‘average out’ over isolates 1,...,n from which baseline population allele frequencies are derived.

* **Ungenotyped clones in the paired baseline isolate b**: the allelic states of undetected clones in the paired baseline isolate are imputed using allele frequencies derived over the unpaired baseline isolates 1,...,n, with a truncated binomial model for the number of clones genotyped per locus.

* **Ungenotyped clones in the recurrent isolate r**: observed alleles in the recurrent isolate are allocated over successfully genotyped clones only, with a truncated multinomial model for the number of clones genotyped per locus.
Overview of the statistical model

PfRecur classifies infections based on an analytically tractable Bayesian model constructed around the set of assumptions detailed in Box 1. Under the model, each isolate is treated as a set of genetically-distinct clones; the cardinality of this set is referred to hereafter as the multiplicity of infection (MOI). For each isolate, we observe a set of alleles $G_\ell$ at loci $\ell \in \{1, \dots, L\}$, with locus-wise cardinality $M_\ell = |G_\ell|$. We take the MOI to be the maximum cardinality over all loci: $M = \max_{1 \le \ell \le L} M_\ell$. Given bulk genotypic data with no quantitative information on the intra-isolate abundance of each allele<fnref n="1" />, we assume a uniform distribution over the sets of allelic configurations that are compatible with $G_\ell$. We assume independence between clones within isolates, between non-recrudescent clones across isolates, and between loci.

The indices $r$, $b$ and $i$ are used to distinguish isolates that are handled differently under the model. Index $r$ corresponds to the recurrent parasite isolate for which classification is being performed; index $b$ corresponds to the paired baseline isolate (i.e., it is from the same study participant as isolate $r$). Index $i = 1, \dots, n$ iterates over $n$ baseline isolates that are not paired with isolate $r$: they are from different study participants.

The recurrent isolate $r$ (of MOI $M^{(r)}$) is modelled as a mixture of $M_C$ recrudescent clones and $M_I = M^{(r)} - M_C$ newly-inoculated clones. For each recurrent isolate $r$, the target of inference under the PfRecur statistical model is the posterior distribution over $M_C$, with state space $\{0, \dots, M^{(r)}\}$:

$$ \mathbb{P}(M_C = m | G^{(r)}) = \frac{\mathbb{P}(M_C = m) \prod_{\ell=1}^{L} \mathbb{P}(G_{\ell}^{(r)} | M_C = m, M_I = M^{(r)} - m)}{\sum_{m=0}^{M^{(r)}} \mathbb{P}(M_C = m) \prod_{\ell=1}^{L} \mathbb{P}(G_{\ell}^{(r)} | M_C = m, M_I = M^{(r)} - m)}. \quad (1) $$

Figure 6 depicts the model structure. The likelihood of the locus-wise recurrent data, $G_\ell^{(r)}$, is conditional on the random variables $M_C$ and $M_I = M^{(r)} - M_C$. It is also conditional on the locus-wise baseline data, $G_\ell^{(b)}$ and $\{G_\ell^{(i)}\}_{i=1..n}^{i \ne b}$, and on the recurrent and baseline MOIs, which are derived from the complete recurrent and baseline data. Although these variables are data derived, they are not treated as random variables under the model. The likelihood also features other user-defined inputs that are not treated as random variables. These include a genotyping detection rate, $\omega$, and a locus-wise genotyping error matrix, $\delta_\ell$. The prior on $M_C$ is conditional on $M^{(r)}$ and on a prior parameter $\beta$. There is no prior on $M_I$ in Equation (1) since $M_I$ is deterministic given $M^{(r)}$ and a realisation of $M_C$. Thanks to an analytically tractable likelihood (next section), Equation (1) is analytically tractable. As such, inference under the PfRecur framework does not require a numerical sampler or an optimisation algorithm.

<footnote n="1">¹Measures of allelic abundance are not readily available for length polymorphic markers that are genotyped using electrophoresis and currently used to perform molecular correction. However, amplicon sequencing, which is recognised by the WHO as a potential future standard for molecular correction, generates read count data on multi-allelic markers (microhaplotypes), and read count data can be used to estimate the relative abundance of within-isolate alleles.</footnote>
<figure><img src="image_7.png" /><figcaption>**Figure 6:** A graphical representation of the Bayesian statistical model within the *PfRecur* framework. Random variables are circled. Data driven quantities are shown in teal; user-specified values are shown in violet.</figcaption></figure>

## Model likelihood

We derive multiple expressions in Appendix A (see Equations (24), (26), (27) and (28)), which together can be used to evaluate the likelihood:

$$ \mathbb{P}(G_{\ell}^{(r)} | M_c = m, M_I = M^{(r)} - m). \qquad (2) $$

The various steps used to conceptualize the likelihood are sketched out below. Throughout, newly-inoculated clones are drawn from population *I* (the contemporaneous population at large), which is approximated under the model by the baseline isolates *i* = 1, ..., *n*. The recrudescent clones are drawn from population *C*, comprised of clones in the paired baseline isolate *b*. In an intermediate step, the probability of $G_{\ell}^{(r)}$ is modelled conditional on allele frequencies, whereby the probability that a clone drawn from population *S* ∈ {*I*, *C*} harbours allele *α* at locus *ℓ* is equated to the population allele frequency $\theta_{\ell}^{(S)}$(*α*). In later steps, allele counts supersede allele frequencies.

## Modelling the number of clones genotyped per locus

In the allelic data of a given isolate, some loci have cardinalities lower than the MOI. The lower cardinality could occur because multiple clones within the isolate share identical alleles at this locus, because the MOI is supuriously elevated by genotyping errors, or because some clones are not genotyped at this locus. We model the number of genotyped clones per locus in the recurrent isolate *r* at loci with cardinalities strictly less than the observed MOI (i.e. $|G_{\ell}^{(r)}| < M^{(r)}$), thereby allowing some clones in the recurrent isolate to evade detection at some loci (multiple clones sharing identical isolates is also addressed — see next section). A spuriously elevated MOI is not accounted for.
We construct a model of the number of clones genotyped per locus in the recurrent isolate $r$, i.e., the number of clones that contribute to the observation $G_{\ell}^{(r)}$, as follows. Each clone in isolate $r$, irrespective of whether it is drawn from population $I$ or $C$, is detected with probability $\omega$ at locus $\ell$. Then, given $M_C = m$ and $M_I = M^{(r)} - m$, the number of detected clones $Q_{\ell}^{(r,S)}$ derived from populations $S \in \{I, C\}$ at locus $\ell$ follow the truncated multinomial distribution

$$ \mathbb{P}(Q_{\ell}^{(r,C)} = q_C, Q_{\ell}^{(r,I)} = q_I, | M_C = m, M_I = M^{(r)} - m) = \frac{\binom{m}{q_C} \binom{M^{(r)}-m}{q_I} \omega^{q_I+q_C} (1-\omega)^{M^{(r)}-q_I-q_C}}{\sum_{j=M_{\ell}^{(r)}}^{M^{(r)}} \binom{M^{(r)}}{j} \omega^j (1-\omega)^{M^{(r)}-j}} \quad \text{for } |G_{\ell}^{(r)}| \le q_I + q_C \le M^{(r)}. \quad (3) $$

## Allocating observed alleles to genotyped clones

We average analytically over compatible allelic configurations in isolate $r$ using the inclusion-exclusion principle. Suppose there are $Q_{\ell}^{(r,S)} = q_S$ genotyped clones drawn from populations $S \in \{I, C\}$ at locus $\ell$ in isolate $r$. Then the probability of observing the set of alleles $G_{\ell}^{(r)}$ at locus $\ell$ takes the form

$$ \mathbb{P}(G_{\ell}^{(r)} | \theta_{\ell}^{(C)}, \theta_{\ell}^{(I)}, Q_{\ell}^{(r,C)} = q_C, Q_{\ell}^{(r,I)} = q_I) = \sum_{A \in \mathcal{P}(G_{\ell}^{(r)})} (-1)^{|A|} \left( \sum_{\alpha \in G_{\ell}^{(r)} \setminus A} \theta_{\ell}^{(C)}(\alpha) \right)^{q_C} \left( \sum_{\alpha \in G_{\ell}^{(r)} \setminus A} \theta_{\ell}^{(I)}(\alpha) \right)^{q_I} \quad (4) $$

where $\mathcal{P}(G_{\ell}^{(r)})$ denotes the power set of $G_{\ell}^{(r)}$; that is, the set of all subsets of $G_{\ell}^{(r)}$ encompassing the empty set. Convolving Equation (4) over the locus-wise truncated multinomial model of genotyped clones (3) yields the probability of the observation $G_{\ell}^{(r)}$ at locus $\ell$ for the recurrent isolate $r$, formulated with respect to the population allele frequencies $\theta_{\ell}^{(S)}$, $S \in \{C, I\}$.

## Modelling allele frequencies

The derivation of population allele frequencies under our model is two-fold. Denote by $H_{\ell}$ the set of possible alleles at locus $\ell$ (equipped with an arbitrary ordering) and by $\theta_{\ell}^{(S)}$ the vector of population allele frequencies over $H_{\ell}$. We begin by deriving allele frequencies $\theta_{\ell}^{(S)}$ conditional on the vector of per-isolate allele counts $C_{\ell}$ for each baseline isolate $b$ or $1, \dots, n$ over $H_{\ell}$; and later address the distribution of per-isolate allele counts $C_{\ell}$, which are not directly observed using bulk genotypic data.

For population $I$, we derive allele frequencies under a Bayesian multinomial-Dirichlet model of baseline isolates $i = 1, \dots, n$ (which excludes the paired baseline isolate $b$). Allele frequencies are formulated over clones in the baseline isolates $i = 1, \dots, n$, whereby each isolate is effectively weighted by its MOI (i.e., high MOI isolates are more informative). Under the assumption of clone-wise independence (both within and across isolates), we obtain the multinomial likelihood

$$ \sum_{i=1}^{n} C_{\ell}^{(i)} | \theta_{\ell}^{(I)} \sim \text{Multinomial} \left( \sum_{i=1}^{n} M^{(i)}, \theta_{\ell}^{(I)} \right). $$
Taking a uniform prior for $\theta_{\ell}^{(I)}$ over the $|H_{\ell}| - 1$ simplex yields the posterior

$$ \theta_{\ell}^{(I)} \mid \sum_{i=1}^{n} \mathbf{C}_{\ell}^{(i)} \sim \text{Dirichlet} \left( \sum_{i=1}^{n} M^{(i)}, \sum_{i=1}^{n} \mathbf{C}_{\ell}^{(i)} + \mathbf{1} \right). \qquad (5) $$

For population C, allele frequencies are largely informed by the paired baseline isolate b, but not exclusively because the alleles of ungenotyped clones in isolate b are imputed using population I allele frequencies, which are based on isolates $i = 1, \dots, n$ (Equation (5)). This imputation explains the omission of isolate b in the formulation of allele frequencies for population I. Akin to the model of genotyped clones for isolate r (Equation (3)), the number of genotyped clones $Q_{\ell}^{(b)}$ in isolate b that contribute to the observation $G_{\ell}^{(b)}$ is modelled to follow a truncated binomial distribution, with support $|G_{\ell}^{(b)}|, \dots, M^{(b)}$ and success probability $\omega$. Given $Q_{\ell}^{(b)} = q$, we denote by $\mathbf{C}_{\ell}^{(b,q)}$ the vector of allele counts over the q clones in isolate b that are genotyped at locus $\ell$. We model $\theta_{\ell}^{(C)}$ as a deterministic function of $Q_{\ell}^{(b)} = q$, the per-isolate allele counts $\mathbf{C}_{\ell}^{(b,q)}$ of the genotyped clones in isolate b, and $\theta_{\ell}^{(I)}$,

$$ \theta_{\ell}^{(C)} = \frac{\mathbf{C}_{\ell}^{(b,q)}}{M^{(b)}} + \left(1 - \frac{q}{M^{(b)}}\right) \cdot \theta_{\ell}^{(I)}. \qquad (6) $$

## Modelling genotyping errors

We account for genotyping error when modelling the per-isolate baseline allele counts $\mathbf{C}_{\ell}^{(i)}$ and $\mathbf{C}_{\ell}^{(b,q)}$. Genotyping errors are modelled using a user-specified right-stochastic error matrix $\delta_{\ell}$ for each locus $\ell$, with the interpretation that $\delta_{\ell}(\alpha, \alpha')$ yields the probability that an allele called as $\alpha$ in a baseline isolate b or $i = 1, \dots, n$ matches an allele called as $\alpha'$ in a recurrent isolate r at locus $\ell$. The PfRecur framework is marker-agnostic, in that the non-parametric error model $\delta_{\ell}$ can be adapted to different marker types. In the present study, we consider a normalised geometric model adapted to length polymorphic microsatellite markers (akin to [8]), parametrised by a genotyping error probability $\varepsilon$ where, for allele lengths $i, j$

$$ \delta_{\ell}(i, j) = \begin{cases} \frac{\varepsilon^{|i-j|+1}}{\sum_{k \neq i} \varepsilon^{|i-k|}} & \text{if } j \neq i \\ 1 - \varepsilon & \text{if } j = i \end{cases}, \qquad (7) $$

where the sum in the denominator is taken over the allelic lengths in the set $H_{\ell}$.

## Deriving moments of allele counts

The locus-wise probability (4) of observed genotypes for recurrent isolate r is formulated as a multinomial expression of the population allele frequencies $\theta_{\ell}^{(S)}$. Convolving Equation (4) over the modelled distributions of $\theta_{\ell}^{(I)}$ (5) and $\theta_{\ell}^{(C)}$ (6), which are conditioned on per-isolate allele counts, yields a multinomial expression of baseline allele counts (Equation (19)).

Because individual clones and the alleles they carry are not directly observable in multiclonal isolates genotyped using standard methods (e.g., not single-cell genotyping), allele counts for multiclonal isolates must be derived under an appropriate model [20–24]. By Jensen's inequality, the expected per-isolate
allele counts $\mathbb{E}[C_{\ell}^{(z)}]$ (which are straightforward to compute) cannot be directly substituted for $C_{\ell}^{(z)}$ in Equation (19). Convolving the locus-wise probability (19) over compatible allelic configurations in baseline isolates $b$ and $i = 1, \dots, n$ requires calculating moments of the per-isolate allele counts

$$ \mathbb{E} \left[ \left( \sum_{\alpha \in A} C_{\ell}^{(b,q)}(\alpha) \right)^k \right] = \sum_{k=1}^{M^{(i)}} k^m \cdot \mathbb{P} \left( \sum_{\alpha \in A} C_{\ell}^{(b,q)}(\alpha) = k \right) $$

and the sample allele counts (i.e., per-isolate allele counts summed over one or more baseline isolates)

$$ \mathbb{E} \left[ \left( \sum_{i=1}^{n} \sum_{\alpha \in A} C_{\ell}^{(i)}(\alpha) \right)^k \right] = \sum_{k=1}^{M^{(1)}+\dots+M^{(n)}} k^m \cdot \mathbb{P} \left( \sum_{i=1}^{n} \sum_{\alpha \in A} C_{\ell}^{(i)}(\alpha) = k \right) $$

aggregated over allelic subsets $A \subseteq H_{\ell}$ (Appendix A.2).

In Appendix A.3, we derive these moments from first principles under our framework. In brief, we use a simple combinatorial argument based on ordered partitions to derive probability mass functions and consequently cumulants of the per-isolate allele counts. To adjust per-isolate allele counts for genotyping error, we adopt a Poisson binomial model: the number of distinct alleles in a given set $A \subset H_{\ell}$ that are harboured by clones in a baseline isolate is modelled as a sum of $|G_{\ell}|$ independent, but not identically distributed, Bernoulli random variables with respective success probabilities $\sum_{\alpha \in A} \delta_{\ell}(g, \alpha)$, $g \in G_{\ell}^{(q)}$.

To compute moments of the sample allele counts (under the assumed independence of clones within and between baseline isolates), we first sum cumulants of the per-isolate allele counts to recover cumulants of the sample allele counts. We then exploit complete exponential Bell polynomials to map cumulants of the sample allele counts to moments of the sample allele counts, as required. We adopt this construction because the order of these moments (up to $M^{(r)}$) is likely, in practical settings, to be smaller than the number of baseline isolates $n$ over which allele counts are aggregated.

## Prior of the statistical model

Since mixtures of newly-inoculated and recrudescent clones are comparatively unlikely in low transmission transmission settings, we take a symmetric prior, which weights pure reinfections and pure recrudescences more heavily than intermediate mixtures. We implement this through a symmetric beta binomial distribution

$$ M_C \sim \text{BetaBinomial}(M^{(r)}, \beta, \beta) $$

with $0 < \beta \le 1$. In the case $\beta = 1$, we recover a uniform distribution over the breakdown of newly-inoculated vs recrudescent clones. In the limit $\beta \to 0$, the prior probability of all intermediate mixtures approaches zero, whereby reinfection and recrudescence constitute mutually exclusive categories.
## Posterior metrics

We generate two posterior metrics for each recurrence: the posterior probability of at least one recrudescent clone

$$M1 := 1 - \mathbb{P}(M_C = 0 | G^{(r)}, G^{(b)}, G^{(1)}, \dots, G^{(n)}) \quad (8)$$

and the posterior expected proportion of recrudescent clones

$$M2 := \frac{1}{M^{(r)}} \sum_{m=0}^{M^{(r)}} m \cdot \mathbb{P}(M_C = m | G^{(r)}, G^{(b)}, G^{(1)}, \dots, G^{(n)}). \quad (9)$$

## Application of posterior metrics

For conceptual consistency with the CDC model [8], we perform classification using metric M1 of *PfRecur*: a recurrent isolate *r* is classified as a recrudescence if $M1 > 0.5$, and as a reinfection otherwise. Down-stream efficacy estimates are computed using metric M1, rather than dichotomised classifications, in line with recommendations from the CDC [8, 17]. Metric M2 serves as a supplementary descriptor for recurrences comprising a mixture of newly-inoculated and recrudescent clones, and is used to validate the model against simulated data for which these mixtures are known (see below).

## R software

*PfRecur* is implemented as an R package (available at https://github.com/somyamehra/PfRecur). This largely relies on base R functionality [35], with additional dependencies on *copula::Stirling2* and *copula::Stirling1* [36] (to evaluate Stirling numbers of the second and unsigned first kind respectively); *PDQutils::cumulant2moment* and *PDQutils::moment2cumulant* [37] (to map between cumulants and moments respectively); *poisbinom::dpoisbinom* [38] (to evaluate the density function of the Poisson binomial distribution); and *VGAM::dbetabinom.ab* [39] (to evaluate the density function of the beta binomial prior). The package accommodates isolates of MOI 9 and fewer (to avoid numerical instability when evaluating the posterior).

As input, the *PfRecur* package requires categorical (presence/absence) genotypic data in a list of binary matrices, where each matrix corresponds to a marker; named matrix columns correspond to alleles; named matrix rows correspond to isolates; and matrix elements are set to 1 if the corresponding allele has been detected in the relevant isolate, and 0 otherwise. Additional user-specified parameters include the per-clone marker-wise probability of detection $\omega$, and a list of marker-wise row-stochastic genotyping error matrices $\delta_{\ell}$.

Given a recurrent isolate *r*, paired baseline isolate *b*, and baseline isolates 1, ..., *n*, the function *PfRecur::evaluate_posterior* returns the discrete posterior distribution for $M_C$ over the state space {0, ..., $M^{(r)}$} in addition to metrics M1 (8) and M2 (9), in the form of a named list.
Simulation study

To validate *PfRecur*, we simulate recurrent isolates as mixtures of newly-inoculated and recrudescent clones (Appendix B). In brief, we consider isolates with MOIs up to 9 (mean baseline MOI ≈ 3), genotyped at 7 unlinked multi-allelic markers (each with between 10 and 30 distinct alleles). We permit siblings within isolates, violating the assumed independence of clones under *PfRecur*. Each simulated clone is detected at each marker with probability ω = 0.9. Genotyping error is applied to the set of alleles harboured by detected clones in each baseline isolate with probability ε = 0.05, in accordance with the length-dependent normalised geometric model (7). We simulate 40 baseline datasets, each comprising 25 isolates. For a baseline isolate of MOI *M*⁽ᵇ⁾, we simulate paired recurrences with MOI *M*⁽ᵏ⁾ = 1, ..., 9 and *m* = 0, ..., min{*M*⁽ᵇ⁾, *M*⁽ᵏ⁾} recrudescent clones. We apply our probabilistic classifier *PfRecur* to recover the posterior distribution for the number of newly-inoculated vs recrudescent clones within each simulated recurrent isolate under a uniform prior (β = 1); the underlying parameters ω and ε are assumed to be known. The results in the main text, pertaining to metrics M1 and M2, are aggregated across 28773 simulated recurrences.

Reanalysis of Dimbu *et al*. [17]

Using *PfRecur* we perform a re-analysis of [17], with baseline isolates stratified by study site (whereby
allele frequencies for newly-inoculated clones are modelled to be site-specific). By default, we set $\beta = 0.25$
for the prior; $\omega = 0.9$ for the per-clone marker-wise probability of detection in each baseline/recurrent
pair; and $\varepsilon = 0.05$ for the genotyping error probability, under the normalised geometric model (7). We
additionally perform a sensitivity analysis for $\omega \in [0.75, 1]$ and $\varepsilon \in [0, 0.25]$. Classification is performed
with both the entire 7 microsatellite marker set, and also omitting the *TA109* marker.

We compare posterior metric M1 of *PfRecur* against the “gold-standard” CDC model [8], which was
originally used to analyse [17]. In the present study, we have re-run code provided openly by Plucinski
and colleagues [40] (with 100,000 iterations for the Gibbs sampler); posterior probabilities based on all 7
microsatellite markers may differ from those reported in [17] due to the stochastic nature of the MCMC
algorithm.

False positive recrudescence rates

To estimate false positive rates for calling recrudescence, we generate 500 artificial ‘non-recrudescence’
datasets from [17] by generating random derangements of baseline study participants labels within each
study site, whereby permuted baseline/recurrent pairs cannot be derived from the same individual and
therefore cannot represent recrudescences. The generation of permuted datasets, rather than permuted
pairs, is necessitated by the construction of the CDC model [8]. We perform classification for these
permuted datasets using both *PfRecur* (metric M1) and the CDC model [8, 40] (with 10,000 iterations
for the Gibbs sampler due to computational time constraints). For each permuted dataset, we compute
the false positive recrudescence rate by averaging the posterior probability of recrudescence (under the
CDC model) or metric M1 (under *PfRecur*) over recurrent isolates. In addition, we consider a match-counting approach, treating the presence of one or more shared alleles at least 4 or 5 markers as evidence of recrudescence. For the panel of 7 neutral microsatellites, no WHO-endorsed guidelines are available but [16, 17, 25] support the "≥ 4/7 rule"<fnref n="2" />. To avoid ambiguity, we restrict match-counting classification to baseline/recurrent pairs with at least one allele call at each of the 7 markers.

## Code and data availability

The *PfRecur* framework has been implemented in an eponymous R package, openly available in a GitHub repository: https://github.com/somyamehra/PfRecur.

This study uses open access data that has previously been published by Dimbu *et al.* [17]. Parasite densities and clinical metadata have been retrieved from Supplemental Table S4 of Dimbu *et al.* [17], while genotypic data have been retrieved from an accompanying GitHub repository [40].

We have re-run code provided openly by Dr Mateusz Plucinski and colleagues (with very minor input/output modifications), which implements the model detailed in Plucinski *et al.* [8] and is available in a GitHub repository [40]: https://github.com/MateuszPlucinski/AngolaTES2021.

For completeness, all code and data relevant to this study (including the data of Dimbu *et al.* [17] and the implementation of Plucinski *et al.* [8]) have been collated in a GitHub repository: https://github.com/somyamehra/PfTreatmentFailure.

## References

1. WHO. *World malaria report 2023* (World Health Organization, 2023).

2. Ashley, E. A. *et al.* Spread of artemisinin resistance in Plasmodium falciparum malaria. *New England Journal of Medicine* **371**, 411–423 (2014).

3. Imwong, M. *et al.* The spread of artemisinin-resistant Plasmodium falciparum in the Greater Mekong subregion: a molecular epidemiology observational study. *The Lancet Infectious Diseases* **17**, 491–497 (2017).

4. Balikagala, B. *et al.* Evidence of artemisinin-resistant malaria in Africa. *New England Journal of Medicine* **385**, 1163–1171 (2021).

<footnote n="2">²Current WHO guidelines, tailored to a 3 marker panel, stipulate a strict 3/3 match-counting rule for primary analysis [9]. We do not consider the strict 7/7 match-counting rule for the 7 neutral microsatellite panel used in [17], given the microsatellite marker *TA109* appears to be problematic. There are also several recurrences in [17] (namely, ZL21-292, ZQ21-103 and ZL21-245) where mismatch at a single marker (with a length difference plausibly attributable to either genotyping or human error) supports the use of a relaxed match-counting rule.</footnote>
5. Uwimana, A. et al. Emergence and clonal expansion of in vitro artemisinin-resistant Plasmodium falciparum kelch13 R561H mutant parasites in Rwanda. *Nature Medicine* **26**, 1602–1608 (2020).

6. WHO. *Methods for Surveillance of Antimalarial Drug Efficacy* (World Health Organization, 2009).

7. Nsanzabana, C., Djalle, D., Guérin, P. J., Ménard, D. & González, I. J. Tools for surveillance of anti-malarial drug resistance: an assessment of the current landscape. *Malaria Journal* **17**, 1–16 (2018).

8. Plucinski, M. M., Morton, L., Bushman, M., Dimbu, P. R. & Udhayakumar, V. Robust algorithm for systematic classification of malaria late treatment failures as recrudescence or reinfection using microsatellite genotyping. *Antimicrobial Agents and Chemotherapy* **59**, 6096–6100 (2015).

9. WHO. *Informal consultation on methodology to distinguish reinfection from recrudescence in high malaria transmission areas* (World Health Organization, 2021).

10. Schnoz, A., Beuret, C., Concu, M., Hosch, S., Rutaihwa, L. K., Golumbeanu, M. & Nsanzabana, C. Genotyping methods to distinguish Plasmodium falciparum recrudescence from new infection for the assessment of antimalarial drug efficacy: an observational, single-centre, comparison study. *en. Lancet Microbe*, 100914 (Oct. 2024).

11. Greenhouse, B., Dokomajilar, C., Hubbard, A., Rosenthal, P. J. & Dorsey, G. Impact of transmission intensity on the accuracy of genotyping to distinguish recrudescence from new infection in antimalarial clinical trials. *Antimicrobial Agents and Chemotherapy* **51**, 3096–3103 (2007).

12. Juliano, J. J. et al. Misclassification of drug failure in Plasmodium falciparum clinical trials in southeast Asia. *The Journal of Infectious Diseases* **200**, 624–628 (2009).

13. Messerli, C., Hofmann, N. E., Beck, H.-P. & Felger, I. Critical evaluation of molecular monitoring in malaria drug efficacy trials and pitfalls of length-polymorphic markers. *Antimicrobial Agents and Chemotherapy* **61**, 10–1128 (2017).

14. Felger, I., Snounou, G., Hastings, I., Moehrle, J. J. & Beck, H.-P. PCR correction strategies for malaria drug trials: updates and clarifications. *The Lancet Infectious Diseases* **20**, e20–e25 (2020).

15. Jones, S. et al. Improving methods for analyzing antimalarial drug efficacy trials: molecular correction based on length-polymorphic markers msp-1, msp-2, and glurp. *Antimicrobial Agents and Chemotherapy* **63**, 10–1128 (2019).
16. Jones, S., Plucinski, M., Kay, K., Hodel, E. M. & Hastings, I. M. A computer modelling approach to evaluate the accuracy of microsatellite markers for classification of recurrent infections during routine monitoring of antimalarial drug efficacy. *Antimicrobial Agents and Chemotherapy* **64**, 10–1128 (2020).

17. Dimbu, P. R. et al. Therapeutic response to four artemisinin-based combination therapies in Angola, 2021. *Antimicrobial Agents and Chemotherapy*, e01525-23 (2021).

18. Slater, M., Kiggundu, M., Dokomajilar, C., Kamya, M. R., Bakyaita, N., Talisuna, A., Rosenthal, P. J. & Dorsey, G. Distinguishing recrudescences from new infections in antimalarial clinical trials: major impact of interpretation of genotyping results on estimates of drug efficacy. *The American Journal of Tropical Medicine and Hygiene* **73**, 256–262 (2005).

19. Taylor, A. R., Foo, Y. S. & White, M. T. Plasmodium vivax relapse, reinfection and recrudescence estimation using genetic data. *medRxiv*, 2022–11 (2022).

20. Taylor, A. R. et al. Estimation of malaria haplotype and genotype frequencies: a statistical approach to overcome the challenge associated with multiclonal infections. *Malaria Journal* **13**, 1–11 (2014).

21. Chang, H.-H. et al. THE REAL McCOIL: A method for the concurrent estimation of the complexity of infection and SNP allele frequency for malaria parasites. *PLoS Computational Biology* **13**, e1005348 (2017).

22. Ju, N., Liu, J. & He, Q. SNP-Slice resolves mixed infections: Simultaneously unveiling strain haplotypes and linking them to hosts. *Bioinformatics* **40** (2024).

23. Murphy, M. & Greenhouse, B. MOIRE: a software package for the estimation of allele frequencies and effective multiplicity of infection from polyallelic data. *Bioinformatics* **40**, btae619 (2024).

24. Paschalidis, A., Watson, O. J., Aydemir, O., Verity, R. & Bailey, J. A. coiaf: directly estimating complexity of infection with allele frequencies. *PLOS Computational Biology* **19**, e1010247 (2023).

25. Plucinski, M. M. & Barratt, J. L. Nonparametric Binary Classification to Distinguish Closely Related versus Unrelated Plasmodium falciparum Parasites. *The American Journal of Tropical Medicine and Hygiene* **104**, 1830 (2021).

26. Mulligan, J.-A., Mandike, R., Palmer, N., Williams, H., Abdulla, S., Bloland, P. & Mills, A. The costs of changing national policy: lessons from malaria treatment policy guidelines in Tanzania. *Tropical Medicine & International Health* **11**, 452–461 (2006).
27. Njau, J. D. et al. The costs of introducing artemisinin-based combination therapy: evidence from district-wide implementation in rural Tanzania. *Malaria Journal* **7**, 1–14 (2008).

28. WHO. *Strategy to respond to antimalarial drug resistance in Africa* (World Health Organization, 2022).

29. White, N. J. & Chotivanich, K. Artemisinin-resistant malaria. en. *Clin. Microbiol. Rev.*, e0010924 (Oct. 2024).

30. Taylor, A. R. et al. Resolving the cause of recurrent Plasmodium vivax malaria probabilistically. *Nature Communications* **10**, 5595 (2019).

31. Mehra, S., Neafsey, D. E., White, M. & Taylor, A. R. Systematic bias in malaria parasite relatedness estimation. *bioRxiv*, 2024-04 (2024).

32. Jonckheere, A. R. A distribution-free k-sample test against ordered alternatives. *Biometrika* **41**, 133–145 (1954).

33. Pohlert, T. & Pohlert, M. *Package ‘pmcmr’. R package version 1* (2018).

34. Nei, M. Analysis of gene diversity in subdivided populations. *Proceedings of the national academy of sciences* **70**, 3321–3323 (1973).

35. R Core Team. *R: A Language and Environment for Statistical Computing* R Foundation for Statistical Computing (Vienna, Austria, 2021). https://www.R-project.org/.

36. Ivan Kojadinovic & Jun Yan. Modeling Multivariate Distributions with ContinuousMARGINS Using the copula R Package. *Journal of Statistical Software* **34**, 1–20. https://www.jstatsoft.org/v34/i09/ (2010).

37. Pav, S. E. *PDQutils: PDQ Functions via Gram Charlier, Edgeworth, and Cornish Fisher Approximations* R package version 0.1.6 (2017). https://github.com/shabbychef/PDQutils.

38. Olivella, S. & Shiraito, Y. *poisbinom: A Faster Implementation of the Poisson-Binomial Distribution* R package version 1.0.1 (2017). https://CRAN.R-project.org/package=poisbinom.

39. Yee, T. W. *Vector Generalized Linear and Additive Models: With an Implementation in R* (Springer, New York, USA, 2015).

40. Plucinski, M. *AngolaTES2021* https://github.com/MateuszPlucinski/AngolaTES2021. 2024.

41. Taylor, A. R., Vickers, E. N. & Greenhouse, B. Review of MrsFreqPhase methods: methods designed to estimate statistically malaria parasite multiplicity of infection, relatedness, frequency and phase. *Malaria Journal* **23**, 308 (2024).
42. Gerlovina, I., Gerlovin, B., Rodríguez-Barraquer, I. & Greenhouse, B. Dcifer: an IBD-based method to calculate genetic distance between polyclonal infections. *Genetics* **222**, iyac126 (2022).

43. Hong, Y. On computing the distribution function for the Poisson binomial distribution. *Computational Statistics & Data Analysis* **59**, 41–51 (2013).

44. Taylor, A. R., Jacob, P. E., Neafsey, D. E. & Buckee, C. O. Estimating relatedness between malaria parasites. *Genetics* **212**, 1337–1351 (2019).

45. Terry M. Therneau & Patricia M. Grambsch. *Modeling Survival Data: Extending the Cox Model* ISBN: 0-387-98784-3 (Springer, New York, 2000).

## Acknowledgements

We thank Dr Mateusz Plucinski from Centers for Disease Control and Prevention, Atlanta for thoughtful and constructive criticism of the manuscript.

We thank the authors of Dimbu *et al.* [17] for making their code and data openly available.

JAW is a Sir Henry Dale Fellow funded by the Wellcome Trust (223253/Z/21/Z). NJW is a Principal Research Fellow funded by the Wellcome Trust (093956/Z/10/C). A CC BY or equivalent licence is applied to the author accepted manuscript arising from this submission, in accordance with the grant's open access conditions. ART is a Marie Skłodowska-Curie Fellow (project number 101110393) funded by the European Union. Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union or granting authority. Neither the European Union nor the granting authority can be held responsible for them.

## Author contributions

Methodology: S.M. and J.A.W. Formal analysis, visualisation: S.M. Writing (original draft): S.M., A.R.T., N.J.W and J.A.W. Writing (review and editing): S.M., A.R.T., M.I., N.J.W. and J.A.W. Conceptualisation: N.J.W. and M.I. Supervision: N.J.W. and J.A.W.

## Competing interests

The authors have no competing interests to declare.
# A Derivation of model likelihood

We model an isolate $r$ from a recurrent infection as a mixture of clones drawn from two populations: population $C$ comprising clones in the paired baseline isolate $b$ (this represents reCrudescence); and population $I$, reflecting the contemporaneous at large but approximated by the baseline population for which baseline isolates 1, ..., $n$ (excluding the baseline isolate $b$, which is paired with $r$) constitute a finite sample (this represents reInfection).

Under our model, the probability that a given clone in isolate $r$ exhibits allele $\alpha$ is equated to the frequency of $\alpha$ in the population from which it is derived. For population $I$, allele frequencies are derived across clones in isolates 1, ..., $n$ under a multinomial-Dirichlet model (with a uniform prior over the standard simplex); for population $C$, we instead consider allele frequencies across clones in isolate $b$. We make the simplifying assumptions of marker-wise independence (i.e., we assume that loci are unlinked); neutral markers; and clone-wise independence, both within isolates and between non-recrudescent clones across isolates (i.e., we ignore relatedness, which amplifies allelic sharing between closely-related clones [30]).

The principal difficulty arises because per-isolate genotypic data are generated in bulk. We observe the presence or absence of detection of an allele within an isolate, but not the allocation of alleles to individual clones, and thus not the per-isolate allele counts. Here, we take a uniform distribution over all allelic configurations that are compatible with the set of observed alleles. Further complications include imperfect detectability of clones and inaccurate genotyping. Accordingly, we account for undetected clones in the paired baseline $b$ and recurrent $r$ isolates, under a marker-wise truncated binomial model. The locus-wise probability of observations for isolate $r$ is formulated over successfully genotyped clones only. In contrast, allelic states for undetected clones in isolate $b$ are imputed based on population $I$ allele frequencies (this imputation is the rationale for excluding isolate $b$ from population $I$). We additionally allow for genotyping error in baseline isolates $b$ and 1, ..., $n$ relative to the recurrent isolate $r$, specified by the probability that an allele called in a baseline isolate matches an allele called in isolate $r$.

In the following sections, we derive an analytic locus-wise likelihood of observed genotypes for isolate $r$ that addresses these complications systematically. We first formulate the locus-wise probability of observations for the recurrent isolate $r$ with respect to allele frequencies in populations $I$ and $C$ (Section A.1), followed by per-isolate allele counts across clones in baseline isolates $b$ and 1, ..., $n$ (Section A.2). Characterising per-isolate allele counts from categorical (presence/absence) data on multiclonal isolates is a non-trivial problem that has been addressed previously [20–24], which are all reviewed in [41]. However, to evaluate the likelihood, we note that it is sufficient to compute *moments* of allele counts summed over one or more baseline isolates. This is a fundamentally easier problem, and permits an explicit solution under our simplifying assumptions (Section A.3).

## Notation

We begin by introducing notation and several additional assumptions. We consider isolates that have been genotyped at a panel of $L$ markers. At locus $l \in \{1, \dots, L\}$, for each isolate $z$, we observe a set of
alleles $G_{\ell}^{(z)}$. We denote by $M_{\ell}^{(z)} := |G_{\ell}^{(z)}|$ the observed cardinality (i.e., the number of distinct alleles) at locus $\ell$ for isolate $z$. The MOI for isolate $z$ is taken to be the maximum observed cardinality across the set of $L$ loci, that is,

$$M^{(z)} = \max_{1 \le \ell \le L} |G_{\ell}^{(z)}| \qquad (10)$$

(this assumption can be readily relaxed).

We thus assume that isolate $z$ comprises a composite of $M^{(z)}$ clones. We index these clones $j$ (arbitrary ordering), and denote by $\alpha_{\ell}^{(z,j)}$ the allele carried by the $j$th clone in the $z$th isolate at locus $\ell$. If $M^{(z)} > 1$ we do not observe $\alpha_{\ell}^{(z,j)}$. Assuming perfect detectability across clones and perfect genotyping accuracy, we can at most observe whether or not at least one clone in isolate $z$ carries allele $\alpha$ at locus $\ell$, that is,

$$1\{\alpha \in G_{\ell}^{(z)}\} = 1 - \prod_{j=1}^{M^{(z)}} 1\{\alpha_{\ell}^{(z,j)} \neq \alpha\}.$$

However, we do not observe *how many* or which (under the arbitrary ordering) of the $M^{(z)}$ clones within isolate $z$ carry allele $\alpha$ at locus $\ell$.

We model genotyping error for baseline isolates as follows. Denote by $H_{\ell} \supset G_{\ell}^{(r)} \cup G_{\ell}^{(b)} \cup \bigcup_{i=1}^{n} G_{\ell}^{(i)}$ the set of known alleles at locus $\ell$. We say that an allele identified as $\alpha \in H_{\ell}$ in a baseline isolate $b$ or $1, \dots, n$ matches allele $\alpha' \in H_{\ell}$ in isolate $r$ with probability $\delta_{\ell}(\alpha, \alpha') \ge 0$ where

$$\sum_{\alpha' \in H_{\ell}} \delta_{\ell}(\alpha, \alpha') = 1 \text{ for each } \alpha \in H_{\ell}.$$

For notational convenience, we impose an arbitrary ordering on elements in the set $H_{\ell}$; the notation $\theta_{\ell}$ is used to denote a vector of allele frequencies, while $C_{\ell}$ is used to denote a vector of per-isolate allele counts across the ordered set $H_{\ell}$. $\delta_{\ell}$ can then be formulated as a right stochastic matrix (i.e.,with rows summing to 1). A summary of notation is provided in Table 2.

## A.1 A locus-wise model predicated on allele frequencies

We begin by formulating a locus-wise model of observed alleles $G_{\ell}^{(r)}$ for isolate $r$, with respect to allele frequencies in populations $I$ and $C$. The derivation of allele frequencies is deferred to Section A.2.

To accommodate incomplete detection of clones in the recurrent isolate $r$, we allocate observed alleles to successfully genotyped clones only. We adopt a locus-wise truncated multinomial model for the number of successfully genotyped clones at each locus $\ell$. We assume that each clone in isolate $r$, whether recrudescent or newly-inoculated, is successfully genotyped at locus $\ell$ with equal probability $\omega$. Suppose the recurrent isolate $r$ comprises $M_C = m$ clones from population $C$ and $M_I = M^{(r)} - m$ clones from population $I$. Then the numbers of clones $Q_{\ell}^{(r,S)}$ in isolate $r$ that are derived from populations $S \in \{C, I\}$ and successfully genotyped at locus $\ell$ (that is, contribute to the observation $G_{\ell}^{(r)}$) take the form

$$\mathbb{P}(Q_{\ell}^{(r,C)} = q_C, Q_{\ell}^{(r,I)} = q_I, | M_C = m, M_I = M^{(r)} - m)$$
<table><caption>**Table 2:** Summary of notation</caption><thead><tr><th></th><th>Quantity</th><th>Interpretation</th></tr></thead><tbody><tr><td rowspan="4">Observed</td><td>H<sub>ℓ</sub></td><td>Set of possible alleles at locus ℓ</td></tr><tr><td>G<sub>ℓ</sub><sup>(z)</sup></td><td>Set of alleles observed for isolate z at locus ℓ</td></tr><tr><td>M<sub>ℓ</sub><sup>(z)</sup> = |G<sub>ℓ</sub><sup>(z)</sup>|</td><td>Cardinality at locus ℓ for isolate z</td></tr><tr><td>M<sup>(z)</sup></td><td>Multiplicity of infection for isolate z</td></tr><tr><td rowspan="2">User-specified</td><td>ω</td><td>Per-clone locus-wise probability of being detected at locus ℓ (applied to clones in isolates b and r only)</td></tr><tr><td>δ<sub>ℓ</sub>(α, α')</td><td>Model of genotyping error: probability that an allele identified as α in a baseline isolate b or i = 1, ..., n matches allele α' in a recurrent isolate r locus ℓ</td></tr><tr><td>Estimated</td><td>M<sub>S</sub></td><td>Number of clones in isolate r derived from population S ∈ {I, C}</td></tr><tr><td>Unobserved</td><td>α<sub>ℓ</sub><sup>(z,j)</sup></td><td>True allele harboured by clone j of isolate z at locus ℓ</td></tr><tr><td rowspan="6">Modelled</td><td>θ<sub>ℓ</sub><sup>(P)</sup>(α)</td><td>Frequency of allele α at locus ℓ in population P ∈ {I, C}</td></tr><tr><td>Q<sub>ℓ</sub><sup>(r,S)</sup></td><td>Number of clones in isolate r that are derived from population S ∈ {I, C} and have been detected at locus ℓ</td></tr><tr><td>Q<sub>ℓ</sub><sup>(b)</sup></td><td>Number of clones in isolate b that have been detected at locus ℓ</td></tr><tr><td>C<sub>ℓ</sub><sup>(i)</sup>(α)</td><td>Count of allele α across clones in isolate i = 1, ..., n at locus ℓ</td></tr><tr><td>C<sub>ℓ</sub><sup>(b,q)</sup>(α)</td><td>Count of allele α across the q clones in isolate b that have been genotyped (in)correctly at locus ℓ</td></tr><tr><td>X<sub>ℓ</sub><sup>(z)</sup>(A)</td><td>Adjusting for genotyping error relative to isolate r, the number of alleles observed for baseline isolate z = b or z = 1, ..., n at locus ℓ that lie in the set A ⊂ H<sub>ℓ</sub></td></tr><tr><td rowspan="3">Other</td><td>s(m, k)</td><td>Unsigned Stirling number of the first kind</td></tr><tr><td>S(m, k)</td><td>Stirling number of the second kind</td></tr><tr><td>P(X)</td><td>Power set of X, i.e., the set of all subsets of X encompassing the empty set</td></tr></tbody></table>
$$ = \frac{\binom{m}{q_C} \binom{M^{(r)}-m}{q_I} \omega^{q_I+q_C} (1-\omega)^{M^{(r)}-q_I-q_C}}{\sum_{j=M_{\ell}^{(r)}}^{M^{(r)}} \binom{M^{(r)}}{j} \omega^j (1-\omega)^{M^{(r)}-j}} \quad \text{for } M_{\ell}^{(r)} \le q_I + q_C \le M^{(r)}. \quad (11) $$

Suppose $Q_{\ell}^{(r,C)} = q_C$ clones from population C and $Q_{\ell}^{(r,I)} = q_I$ clones from population I are successfully genotyped at locus $\ell$ in isolate r. As in [23, 42], we use the inclusion-exclusion principle to average over the different ways to allocate alleles in $G_{\ell}^{(r)}$ to $q_C+q_I$ clones; the inclusion-exclusion principle circumvents the need to allocate alleles to clones explicitly. Let $\mathcal{P}(G_{\ell}^{(r)})$ denote the power set of $G_{\ell}^{(r)}$, that is, the set of all subsets of $G_{\ell}^{(r)}$ (encompassing the empty set). Then the conditional probability of the observation $G_{\ell}^{(r)}$ at locus $\ell$ is

$$ \mathbb{P}\left(G_{\ell}^{(r)} \mid \theta_{\ell}^{(C)}, \theta_{\ell}^{(I)}, Q_{\ell}^{(r,C)} = q_C, Q_{\ell}^{(r,I)} = q_I\right) = \sum_{A \in \mathcal{P}(G_{\ell}^{(r)})} (-1)^{|A|} \left(\sum_{\alpha \in G_{\ell}^{(r)} \setminus A} \theta_{\ell}^{(C)}(\alpha)\right)^{q_C} \left(\sum_{\alpha \in G_{\ell}^{(r)} \setminus A} \theta_{\ell}^{(I)}(\alpha)\right)^{q_I}, \quad (12) $$

where $\theta_{\ell}^{(S)}(\alpha)$ is the frequency of allele $\alpha$ at locus $\ell$ in population $S \in \{I, C\}$.

Convolving Equation (12) the distribution of successfully genotyped clones (11) then yields

$$ \mathbb{P}\left(G_{\ell}^{(r)} \mid \theta_{\ell}^{(C)}, \theta_{\ell}^{(I)}, M_C = m, M_I = M^{(r)} - m\right) = \sum_{q_C=0}^{M^{(r)}} \sum_{q_I=M_{\ell}^{(r)}-q_C}^{M^{(r)}-q_C} \frac{\binom{m}{q_C} \binom{M^{(r)}-m}{q_I} \omega^{q_I+q_C} (1-\omega)^{M^{(r)}-q_I-q_C}}{\sum_{j=M_{\ell}^{(r)}}^{M^{(r)}} \binom{M^{(r)}}{j} \omega^j (1-\omega)^{M^{(r)}-j}} \left\{ \sum_{A \in \mathcal{P}(G_{\ell}^{(r)})} (-1)^{|A|} \left( \sum_{\alpha \in G_{\ell}^{(r)} \setminus A} \theta_{\ell}^{(C)}(\alpha) \right)^{q_C} \left( \sum_{\alpha \in G_{\ell}^{(r)} \setminus A} \theta_{\ell}^{(I)}(\alpha) \right)^{q_I} \right\}. \quad (13) $$

**Remark:** Rather than modelling (un)genotyped clones in recurrent isolate r directly, we could have taken a uniform distribution over all allelic configurations of $M^{(r)}$ clones such that each allele in the observed set $G_{\ell}^{(r)}$ is carried by at least one clone at locus $\ell$; this would implicitly allow for the incomplete detection of clones in isolate r. Under this construction, the conditional probability of the observation $G_{\ell}^{(r)}$ at locus $\ell$ would take the form

$$ \mathbb{P}_{\text{other}}\left(G_{\ell}^{(r)} \mid \theta_{\ell}^{(C)}, \theta_{\ell}^{(I)}, M_C = m, M_I = M^{(r)} - m\right) = \sum_{A \in \mathcal{P}(G_{\ell}^{(r)})} (-1)^{|A|} \left(\sum_{\alpha \in H_{\ell} \setminus A} \theta_{\ell}^{(C)}(\alpha)\right)^m \left(\sum_{\alpha \in H_{\ell} \setminus A} \theta_{\ell}^{(I)}(\alpha)\right)^{M^{(r)}-m} \quad (14) $$

with differences to Equation (12) highlighted in red. We have not adopted this construction because it is liable to generate false positive signals of recrudescence when the recurrent isolate r has MOI>1, and exhibits alleles matching the paired baseline isolate b at a subset of loci.
## A.2 Modelling allele frequencies

Here, we address the derivation of allele frequencies for both populations I and C. By convolving Equation (13) over the modelled distribution of population allele frequencies, we show that the locus-wise probability of observations $G_{\ell}^{(r)}$ for isolate r can be formulated as a multinomial expression of sample allele counts (i.e, allele counts summed over one or more baseline isolates b and $i = 1, \dots, n$). The consequences of genotyping error for baseline isolates relative to the recurrent isolate r, specified by the probability that alleles called in baseline isolates b and $i = 1, \dots, n$ match alleles called in the recurrent isolate r, will be addressed in Section A.3 in the sub-model of per-isolate allele counts.

### A.2.1 Allele frequencies for population I

Treating isolates 1, ..., n as a finite sample from a larger population I, we derive allele frequencies for population I under a Bayesian multinomial-Dirichlet model. We weight each isolate by its MOI, whereby each clone in isolates $i = 1, \dots, n$ is weighted equally.

For notational convenience, we introduce an arbitrary ordering on $H_{\ell}$, the set if possible alleles at locus $\ell$. For the $k^{th}$ allele $\alpha_k$ in the set $H_{\ell}$, we denote by $\theta_{\ell}^{(I)}(k)$ the baseline population frequency, and by

$$C_{\ell}^{(i)}(k) = \sum_{j=1}^{M^{(i)}} 1\{\alpha_{\ell}^{(i,j)} = \alpha_k\} \quad (15)$$

the per-isolate allele count $\in \{0, 1, \dots, M^{(i)}\}$ across clones in isolate 1, ..., n.

Under the assumed independence of clones within and between the baseline isolates 1, ..., n, we model sample allele counts aggregated over isolates 1, ..., n to follow a multinomial distribution

$$\sum_{i=1}^{n} C_{\ell}^{(i)} \mid \theta_{\ell}^{(I)} \sim \text{Multinomial} \left( \sum_{i=1}^{n} M^{(i)}, \theta_{\ell}^{(I)} \right)$$

governed by the population allele frequencies $\theta_{\ell}^{(I)}$.

Taking uniform prior for $\theta_{\ell}^{(I)}$ over the $|H_{\ell}| - 1$ simplex, Dirichlet($|H_{\ell}|, \mathbf{1}$), the posterior distribution of allele frequencies $\theta_{\ell}^{(I)}$ given the per-isolate allele counts $C_{\ell}^{(i)}$ takes the form

$$\theta_{\ell}^{(I)} \mid \sum_{i=1}^{n} C_{\ell}^{(i)} \sim \text{Dirichlet} \left( \sum_{i=1}^{n} M^{(i)}, \sum_{i=1}^{n} C_{\ell}^{(i)} + \mathbf{1} \right). \quad (16)$$

### A.2.2 Allele frequencies for population C

Allele frequencies for population C are largely governed by clones in the paired baseline isolate b. However, we accommodate undetected clones in isolate b, with allelic states for ungenotyped clones imputed using population I allele frequencies $\theta_{\ell}^{(I)}$. This imputation is the rationale for withholding the paired isolate b when deriving population I allele frequencies.
Akin to the model of (un)genotyped clones for isolate r (Equation (11)), we model each clone in isolate b to be detected at locus $\ell$ with probability $\omega$. Then the number of clones $Q_{\ell}^{(b)}$ in isolate b that are successfully genotyped at locus $\ell$ (that is, contribute to the observation $G_{\ell}^{(b)}$) is modelled to follow a truncated binomial distribution, with the lower bound of the support given by the observed cardinality $M_{\ell}^{(b)}$ at locus $\ell$:

$$ \mathbb{P}(Q_{\ell} = q | M_{\ell}^{(b)}, M^{(b)}, \omega) = \frac{\binom{M_{\ell}^{(b)}}{q} \omega^q (1-\omega)^{M_{\ell}^{(b)}-q}}{\sum_{j=M_{\ell}^{(b)}}^{M_{\ell}^{(b)}} \binom{M_{\ell}^{(b)}}{j} \omega^j (1-\omega)^{M_{\ell}^{(b)}-j}} \quad \text{for } q = M_{\ell}^{(b)}, \dots, M^{(b)}. \quad (17) $$

Given $Q_{\ell}^{(b)} = q$, let $C_{\ell}^{(b,q)}$ be the vector of allele counts across the q clones in isolate b that have been successfully genotyped at locus $\ell$. Then, conditional on $Q_{\ell}^{(b)} = q$, $C_{\ell}^{(b,q)}$ and $\theta_{\ell}^{(I)}$, we model

$$ \theta_{\ell}^{(C)} \mid \{Q_{\ell}^{(b)} = q, C_{\ell}^{(b,q)}, \theta_{\ell}^{(I)}\} = \frac{C_{\ell}^{(b,q)}}{M^{(b)}} + \left(1 - \frac{q}{M^{(b)}}\right) \cdot \theta_{\ell}^{(I)}. \quad (18) $$

**Remark:** Assuming perfect detection of clones in the paired baseline isolate b (i.e., setting $\omega = 1$) can have strong implications when paired recurrences share alleles at all but one locus, with a length difference that cannot plausibly be explained by genotyping error, but likely reflects human error; see Figure 7 for two empirical examples from [17]. Accommodating ungenotyped clones in isolate b can prevent a single locus from having a disproportionately large effect on overall classifications.

### A.2.3 A locus-wise model predicated on allele counts

Using the derived distributions of allele frequencies for population I, $\theta_{\ell}^{(I)}$ (Equation (16)), and population C, $\theta_{\ell}^{(C)}$ (Equation (18)), we can formulate the locus-wise probability of observations $G_{\ell}^{(r)}$ for isolate r in terms of sample allele counts.

Substituting Equation (18) into Equation (13) and convolving over the truncated binomial distribution $Q_{\ell}^{(b)}$ (17) governing the number of clones detected at locus $\ell$ in isolate b, we find that:

$$ \begin{aligned} & \mathbb{P}(G_{\ell}^{(r)} \mid \hat{\theta}_{\ell}^{(I)}, C_{\ell}^{(b,q)}, M_C = m, M_I = M^{(r)} - m) \\ & = \frac{1}{J_{\ell}^{(b)} J_{\ell}^{(r)}} \sum_{q_B=M_{\ell}^{(b)}}^{M^{(b)}} \sum_{q_C=0}^{M^{(r)}} \sum_{q_I=M_{\ell}^{(r)}-q_C}^{M^{(r)}-q_C} \binom{M^{(b)}}{q_B} \binom{m}{q_C} \binom{M^{(r)}-m}{q_I} \omega^{q_I+q_C+q_B} (1-\omega)^{M^{(r)}+M^{(b)}-q_I-q_C-q_B} \\ & \qquad \sum_{k=0}^{q_C} \binom{q_C}{k} \left(\frac{q_B}{M^{(b)}}\right)^k \left(1-\frac{q_B}{M^{(b)}}\right)^{q_C-k} \sum_{A \in \mathcal{P}(G_{\ell}^{(r)})} (-1)^{|A|} \left(\sum_{\alpha \in G_{\ell}^{(r)} \setminus A} \frac{C_{\ell}^{(b,q_B)}(\alpha)}{q_B}\right)^k \left(\sum_{\alpha \in G_{\ell}^{(r)} \setminus A} \theta_{\ell}^{(I)}(\alpha)\right)^{q_I+q_C-k} \end{aligned} \quad (19) $$

where we set

$$ J_{\ell}^{(s)} = \sum_{j=M_{\ell}^{(s)}}^{M^{(s)}} \binom{M^{(s)}}{j} \omega^j (1-\omega)^{M^{(s)}-j}. \quad (20) $$
Convolving Equation (19) over the derived distribution of $\theta_{\ell}^{(I)}$ (16) (with integration performed over the $|H_{\ell}| - 1$ simplex) yields the locus-wise probability of observations $G_{\ell}^{(r)}$ for isolate $r$ formulated with respect to the per-isolate allele counts $C_{\ell}^{(b,q)}$ across successfully genotyped clones in isolate $b$, and the sample allele counts $\sum_{i=1}^{n} C_{\ell}^{(i)}$:

$$ \begin{aligned} & \mathbb{P}\left(G_{\ell}^{(r)} \middle| \sum_{i=1}^{n} C_{\ell}^{(i)}, C_{\ell}^{(b,q)}, M_C = m, M_I = M^{(r)} - m\right) \\ &= \int_{\boldsymbol{\theta}_{\ell}} \mathbb{P}\left(G_{\ell}^{(r)} \middle| \boldsymbol{\theta}_{\ell}^{(I)}, C_{\ell}^{(b,q)}, M_C = m, M_I = M^{(r)} - m\right) \cdot f\left(\boldsymbol{\theta}_{\ell} \middle| \sum_{i=1}^{n} C_{\ell}^{(i)}\right) d\boldsymbol{\theta}_{\ell} \\ &= \frac{1}{J_{\ell}^{(b)} J_{\ell}^{(r)}} \sum_{q_B=M_{\ell}^{(b)}}^{M^{(b)}} \sum_{q_C=0}^{M^{(r)}} \sum_{q_I=M_{\ell}^{(r)}-q_C}^{M^{(r)}-q_C} \binom{M^{(b)}}{q_B} \binom{M^{(r)}-m}{q_I} \omega^{q_I+q_C+q_B} (1-\omega)^{M^{(r)}+M^{(b)}-q_I-q_C-q_B} \\ & \qquad \sum_{k=0}^{q_C} \binom{q_C}{k} \left(\frac{q_B}{M^{(b)}}\right)^k \left(1-\frac{q_B}{M^{(b)}}\right)^{q_C-k} \frac{\Gamma\left(\sum_{i=1}^{n} M^{(i)} + |H_{\ell}|\right)}{\Gamma\left(\sum_{i=1}^{n} M^{(i)} + |H_{\ell}| + q_I + q_C\right)} \\ & \qquad \left\{ \sum_{A \in \mathcal{P}(G_{\ell}^{(r)})} (-1)^{|A|} \left( \sum_{\alpha \in G_{\ell}^{(r)} \setminus A} \frac{C_{\ell}^{(b,q_B)}(\alpha)}{q_B} \right)^k \left( |H_{\ell}| - |A| + \sum_{\alpha \in G_{\ell}^{(r)} \setminus A} \sum_{i=1}^{n} C_{\ell}^{(i)}(\alpha) \right)^{(q_I+q_C-k)} \right\} \end{aligned} \quad (21) $$

where $(\cdot)^{(\cdot)}$ denotes the rising factorial. Here, we have used the aggregation property of the Dirichlet distribution.

To write Equation (21) in terms of simple powers of $C_{\ell}^{(b,q)}$ and $\sum_{i=1}^{n} C_{\ell}^{(i)}$, we apply the binomial expansion to yield

$$ \left( |H_{\ell}| - |A| + \sum_{\alpha \in G_{\ell}^{(r)} \setminus A} \sum_{i=1}^{n} C_{\ell}^{(i)}(\alpha) \right)^{(q_I+q_C-k)} \\ = \sum_{j=0}^{q_I+q_C-k} \beta_{(q_I+q_C-k, j)} (|H_{\ell}| - |A|) \cdot \left( \sum_{\alpha \in G_{\ell}^{(r)} \setminus A} \sum_{i=1}^{n} C_{\ell}^{(i)}(\alpha) \right)^j \quad (22) $$

where the coefficients

$$ \beta_{(m,j)}(x) = \sum_{k=j}^{m} s(m, k) \binom{k}{j} x^{k-j} \quad (23) $$

are formulated with respect to unsigned Stirling numbers of the first kind $s(m, k)$. Substituting Equation (22) into Equation (21) yields a multinomial expression of the per-isolate allele counts $C^{(b,q)}$ and the sample allele counts $\sum_{i=1}^{n} C_{\ell}^{(i)}$, as desired.

## A.3 Moments of allele counts

Using bulk genotypic data, allele counts are not directly observed for multiclonal isolates. To recover the locus-wise likelihood of observations $G_{\ell}^{(r)}$ for isolate $r$ given it comprises $M_C = m$ recrudescent and $M_I = M^{(r)} - m$ newly-inoculated clones, we must convolve Equation (21) over the underlying isolate-wise distribution of allele counts. Because the locus-wise probability of observations $G_{\ell}^{(r)}$ for isolate
r can be written as a multinomial expression of the per-isolate allele counts $C^{(b,q)}$ and sample allele counts $\sum_{i=1}^{n} C_{\ell}^{(i)}$ (Equations (21) and (22)), evaluation of the locus-wise likelihood distils down to the computation of moments of these allele counts. This is a tractable combinatorial problem which can be solved explicitly under our simplifying assumptions.

Substituting Equation (22) into Equation (21) and convolving over the distribution of per-isolate allele counts for baseline isolates b and i = 1, ..., n, we obtain

$$ 
\begin{aligned}
& \mathcal{L}(G_{\ell}^{(r)} | M_C = m, M_I = M^{(r)} - m) \\
&= \sum_{C_{\ell}^{(b,q)}} \sum_{C_{\ell}^{(i)}} \mathbb{P}\left(G_{\ell}^{(r)} \mid \sum_{i=1}^{n} C_{\ell}^{(i)}, C_{\ell}^{(b,q)}, M_C = m, M_I = M^{(r)} - m\right) \mathbb{P}\left(\sum_{i=1}^{n} C_{\ell}^{(i)}\right) \mathbb{P}(C_{\ell}^{(b,q)}) \\
&= \frac{1}{J_{\ell}^{(b)} J_{\ell}^{(r)}} \sum_{q_B=M_{\ell}^{(b)}}^{M^{(b)}} \sum_{q_C=0}^{M^{(r)}} \sum_{q_I=M_{\ell}^{(r)}-q_C}^{M^{(r)}-q_C} \binom{M^{(b)}}{q_B} \binom{M^{(r)}-m}{q_I} \omega^{q_I+q_C+q_B} (1-\omega)^{M^{(r)}+M^{(b)}-q_I-q_C-q_B} \\
& \qquad \times \sum_{k=0}^{q_C} \binom{q_C}{k} \left(\frac{q_B}{M^{(b)}}\right)^k \left(1-\frac{q_B}{M^{(b)}}\right)^{q_C-k} \frac{\Gamma\left(\sum_{i=1}^{n} M^{(i)} + |H_{\ell}|\right)}{\Gamma\left(\sum_{i=1}^{n} M^{(i)} + |H_{\ell}| + q_I + q_C\right)} \\
& \qquad \times \left\{ \sum_{A \in \mathcal{P}(G_{\ell}^{(r)})} (-1)^{|A|} \mathbb{E}\left[\left(\sum_{\alpha \in G_{\ell}^{(r)} \setminus A} \frac{C_{\ell}^{(b,q_B)}(\alpha)}{q_B}\right)^k\right] \right. \\
& \qquad \qquad \left. \sum_{j=0}^{q_I+q_C-k} \beta_{(q_I+q_C-k, j)} (|H_{\ell}| - |A|) \mathbb{E}\left[\left(\sum_{\alpha \in G_{\ell}^{(r)} \setminus A} \sum_{i=1}^{n} C_{\ell}^{(i)}(\alpha)\right)^j\right] \right\}
\end{aligned} 
\quad (24) 
$$

where we have evoked the assumed independence of (genotyped) clones across baseline isolates b and 1, ..., n.

Evaluating the locus-wise likelihood of observations $G_{\ell}^{(r)}$ for isolate r (24) necessitates the calculation of moments of the collective counts of allelic subsets $A'$ summed over one or more baseline isolates. We can derive these moments explicitly under three key assumptions:

1. We consider the distribution of per-isolate allele counts $C_{\ell}^{(i)}$ for each isolate $i = 1, ..., n$, such that all of the $M^{(i)}$ clones in isolate $i$ contribute to the observation $G_{\ell}^{(i)}$ at locus $\ell$, whereby we assume implicitly that the effects of ungenotyped clones (if any) ‘average out’ across isolates $i = 1, ..., n$ when formulating allele frequencies for population I; note that $C_{\ell}^{(b,q)}$ corresponds explicitly to the $Q_{\ell} = q$ clones within isolate b that have been successfully genotyped at locus $\ell$ and therefore contribute to the observation $G_{\ell}^{(b)}$.

2. We model clones to be independent within and between baseline isolates.

3. We take a uniform distribution over compatible allelic configurations within each isolate, that is, any allocation of the alleles $G_{\ell}^{(i)}$ ($G_{\ell}^{(b)}$) over $M^{(i)}$ ($q = M_{\ell}^{(b)}, ..., M^{(b)}$) clones, such that at least one clone harbours each allele in the set $G_{\ell}^{(i)}$ ($G_{\ell}^{(b)}$), is modelled to be equally likely.

As foreshadowed, we adjust allele counts in baseline isolates for genotyping error relative to the recurrent
isolate $r$, in accordance with the non-parametric model $\delta_{\ell}(\alpha, \alpha')$ which yields the probability that an allele identified as $\alpha \in H_{\ell}$ in a baseline isolate matches allele $\alpha' \in H_{\ell}$ in a recurrent isolate $r$.

### A.3.1 Per-isolate allele counts

We begin by characterising per-isolate allele counts $C_{\ell}^{(b,q)}$ for the $q$ clones in isolate $b$ that are successfully genotyped at locus $\ell$. Denote by

$$X_{\ell}^{(b)}(A') = \left| \left\{ \bigcup_{j=1}^{q} \alpha_{\ell}^{(b,j)} \right\} \cap A' \right|$$

the number of distinct alleles, harboured by the $q$ clones in isolate $b$ that are successfully genotyped at locus $\ell$, that lie in a specified set $A' \subset H_{\ell}$. Accommodating genotyping error relative to the recurrent isolate $r$, each observed allele $g \in G_{\ell}^{(b)}$ matches an allele within the set $A' \subset G_{\ell}^{(r)}$ with probability

$$\delta(g, A') = \sum_{\alpha \in A'} \delta_{\ell}(g, \alpha).$$

We model genotyping error independently for each allele in the observed set $G_{\ell}^{(b)}$. Accordingly, we model $X_{\ell}^{(b)}(A')$ as a sum of $M_{\ell}^{(b)}$ independent but not identically distributed Bernoulli random variables (i.e., a Poisson binomial distribution) with respective success probabilities $\delta(g, A')$, $g \in G_{\ell}^{(b)}$.

Define

$$C_{\ell}^{(b,q)}(A') = \sum_{\alpha \in A'} C_{\ell}^{(b,q)}(\alpha) = \sum_{j=1}^{q} 1\{\alpha_{\ell}^{(b,j)} \in A'\} \quad (25)$$

to be the number of clones within isolate $b$ that are genotyped (correctly or incorrectly) at locus $\ell$ and harbour an allele in the set $A'$ at locus $\ell$. We can recover the distribution of $C_{\ell}^{(b,q)}(A')$ using a simple combinatorial argument as follows.

The number of allelic configurations $N_{\ell}^{(b,q)}$ of $q$ clones compatible with the observation $G_{\ell}^{(b)}$ at locus $\ell$ is equivalent to the number of ordered partitions of $q$ objects into $M_{\ell}^{(b)}$ non-empty subsets. We can write this in the form

$$N_{\ell}^{(b,q)} = S(q, M_{\ell}^{(b)}) \cdot M_{\ell}^{(b)!}$$

where $S(n, k)$ denotes the Stirling number of the second kind.

Given $X_{\ell}^{(b)}(A') = h$, allelic configurations of $q$ clones such that $C_{\ell}^{(q,b)}(A') = k \ge h$ clones harbour alleles in the set $A'$ can be generated by:

* Separating the $q$ clones into two groups of size $k$ and $q-k$ respectively, which can be performed in $\binom{q}{k}$ ways;

* For the group of size $k$ (corresponding to clones with alleles in the set $A'$), constructing ordered partitions into $h$ non-empty subsets, which can be performed in $S(k, h) \cdot h!$ ways;
* For the group of size $(q-k)$ (corresponding to clones with alleles outside the set $A'$), constructing ordered partitions into $(M_{\ell}^{(b)} - h)$ non-empty subsets, which can be performed in $S(q-k, M_{\ell}^{(b)} - h) \cdot (M_{\ell}^{(b)} - h)!$ ways.

Therefore, adjusting for genotyping error relative to the recurrent isolate $r$, the number of allelic configurations of $q$ clones, such that the unique set of observed alleles at locus $\ell$ is $G_{\ell}^{(b)}$ and precisely $X_{\ell}^{(b)}(A') = k$ clones harbour alleles within the set $A'$, takes the form

$$N_{\ell}^{(b,q)}(A', k) = \sum_{h=\max\{0, M_{\ell}^{(b)} - q + k\}}^{\min\{k, M_{\ell}^{(b)}\}} \binom{q}{k} \cdot S(k, h)h! \cdot S(q-k, M_{\ell}^{(b)} - h)(M_{\ell}^{(b)} - h)! \cdot \mathbb{P}(X_{\ell}^{(b)}(A') = h).$$

Efficient computation of the probability mass function $\mathbb{P}(X_{\ell}^{(b)}(A') = h)$ of the Poisson binomial distribution can be performed using the DFT method of Hong [43], with an implementation available in the R package `poisbinom` [38]. The consequences of genotyping error, relative to the recurrent isolate $r$, are embedded implicitly within the term $\mathbb{P}(X_{\ell}^{(b)}(A') = h)$.

We thus obtain the probability mass function for the number of clones $C_{\ell}^{(b,q)}(A')$ (of the $q$ clones that are genotyped at locus $\ell$ in isolate $b$) with alleles lying in the set $A'$ at locus $\ell$:

$$\mathbb{P}(C_{\ell}^{(b,q)}(A') = k) = \sum_{h=\max\{0, M_{\ell}^{(b)} - q + k\}}^{\min\{k, M_{\ell}^{(b)}\}} \frac{\binom{q}{k}}{\binom{M_{\ell}^{(b)}-h}{h}} \cdot \frac{S(k, h)S(q-k, M_{\ell}^{(b)} - h)}{S(q, M_{\ell}^{(b)})} \cdot \mathbb{P}(X_{\ell}^{(b)}(A') = h) \quad (26)$$

for $k \in \{0, 1, \dots, q\}$. We can use the PMF to (26) to compute moments of $C^{(b,q)}(A')$ for relevant allelic subsets $A'$ to be substituted into Equation (24).

### A.3.2 Sample allele counts

The likelihood (24) additionally requires computation of the moments of sample counts of allelic sets $A'$:

$$C_{\ell}^{(I)}(A') := \sum_{i=1}^{n} \sum_{\alpha \in A'} C_{\ell}^{(i)}(\alpha).$$

Here, we exploit the fact that the order of these moments (up to $M^{(r)}$) is likely to be smaller than the size of the baseline sample $n$. With minor modifications in notation, we can use the PMF (26) for per-isolate allele counts to directly compute the $y^{th}$ cumulant $\kappa_y(C_{\ell}^{(i)}(A'))$ for the number of clones $C_{\ell}^{(i)}(A') = \sum_{\alpha \in A'} C_{\ell}^{(i)}(\alpha)$ within each isolate $i = 1, \dots, n$ that harbour alleles within the set $A'$. Due to the independence of clones across baseline isolates, it follows that

$$\kappa_y(C_{\ell}^{(I)}(A')) = \sum_{i=1}^{n} \kappa_y(C_{\ell}^{(i)}(A')) \quad (27)$$

To recover the $m^{th}$ moment of $C_{\ell}^{(I)}(A')$, we exploit the mapping between moments and cumulants:

$$\mathbb{E}[C_{\ell}^{(I)}(A')^m] = B_m(\kappa_1(C_{\ell}^{(I)}(A')), \dots, \kappa_m(C_{\ell}^{(I)}(A'))). \quad (28)$$

where $B_m(\cdot)$ denotes the $m^{th}$ complete exponential Bell polynomial. Equation (28) can then be substituted into the likelihood (24) for relevant allelic subsets $A'$.
## B Validation against simulated data

To validate our classification model *PfRecur*, we simulate recurrent isolates as mixtures of newly-inoculated and recrudescent clones. A summary of simulation parameters is provided in Table 3.

<table><caption>**Table 3:** Summary of simulation parameters</caption><thead><tr><td>Parameter</td><td>Interpretation</td><td>Value</td></tr></thead><tbody><tr><td>L</td><td>Number of independent (unlinked) and neutral markers</td><td>7</td></tr><tr><td>h</td><td>Vector of marker cardinalities (number of possible alleles |H<sub>l</sub>| at each marker l = 1, ..., L)</td><td>(30, 30, 20, 20, 20, 10)</td></tr><tr><td>ε</td><td>Probability of genotyping error per allele call in each baseline isolate</td><td>0.05</td></tr><tr><td>M<sub>max</sub></td><td>Maximum MOI (number of not necessarily distinct clones) within each isolate</td><td>9</td></tr><tr><td>μ</td><td>Mean of the geometric distribution (with non-zero state space) governing the isolate MOI (truncated at M<sub>max</sub>)</td><td>3</td></tr><tr><td>p<sub>sib</sub></td><td>Parameter governing within-isolate sibling structure (represented by non-crossing partitions); zero corresponds to the absence of sibling structure, i.e., independent clones</td><td>0.1</td></tr><tr><td>ω</td><td>Per-clone marker-wise probability of detection</td><td>0.9</td></tr><tr><td>n</td><td>Number of baseline isolates</td><td>25</td></tr></tbody></table>

We generate baseline samples comprising *n* independent isolates, accompanied by paired recurrent isolates. For each isolate, we simulate genotypic data across *L* markers, indexed *l* = 1, ..., *L*.

Each marker *l* is uniquely specified by a cardinality $h_l$. The set of possible alleles at marker *l* is denoted $H_l = \{1, ..., h_l\}$, with the integer label interpreted as a length (in arbitrary units).

We accommodate structured genotyping error at each marker *l*. We adopt a normalised geometric model, based on the absolute length difference between an observed and underlying allele [8], with success parameter $(1-\epsilon)$. We normalise this model such that each allele is erroneously genotyped with probability $\epsilon$, that is,

$$ \delta_l(i, j) = \begin{cases} \frac{\epsilon^{|i-j|+1}}{\sum_{k \neq i} \epsilon^{|i-k|}} & \text{if } j \neq i \\ 1 - \epsilon & \text{if } j = i. \end{cases} \qquad (29) $$

We then generate a random permutation $\sigma$ of the set $H_l$. The population frequency of the $\sigma(i)^{th}$ allele at marker *l* is taken to be

$$ \theta_l(\sigma(i)) = \frac{2i}{h_l(h_l + 1)}, $$
yielding the effective cardinality [44]

$$ \text{EffectiveCardinality}_{\ell} = \left( \sum_{i=1}^{h_{\ell}} \theta_{\ell}(i)^2 \right)^{-1} = \frac{3h_{\ell}}{4} + \frac{3}{8} \left( 1 - \frac{1}{1 + 2h_{\ell}} \right). $$

This is motivated by empirical data, which appear to yield an approximately-linear trend in ordered (expected) sample allele frequencies.

The MOI of each isolate (baseline or recurrent) is sampled from the geometric distribution with mean $\mu$ and non-zero state space, truncated at $M_{\max}$, that is,

$$ P(M = m) = \frac{\mu^{M_{\max} - m} (\mu - 1)^{m-1}}{\mu^{M_{\max}} - (\mu - 1)^{M_{\max}}} \quad \text{for } m = 1, \dots, M_{\max}. \qquad (30) $$

We accommodate the presence of siblings within simulated isolates. To generate a baseline set of clones for isolate $i$ with MOI $M^{(i)}$, we first simulate a non-crossing partition of the set $\{1, \dots, M^{(i)}\}$, governed by the partition parameter $p_{\text{sib}}$: given

$$ b_j \stackrel{i.i.d.}{\sim} \text{Bernoulli}(p_{\text{sib}}), \quad j = 1, \dots, (M^{(i)} - 1) $$

clones $j$ and $(j+1)$ are assigned to the same block if and only if $b_j = 0$. For each block of siblings $B$, we then generate a pair of independent parents. The allelic state $\rho_{\ell}^{(q)}(B)$ of parent $q \in \{1, 2\}$ of block $B$ at locus $\ell$ is sampled from a categorical distribution governed by the population allele frequencies

$$ \rho_{\ell}^{(q)}(B) \sim \text{Categorical}(\theta_{\ell}). $$

We then recover each clone in block $B$ as a mosaic of these independent parents, retaining the assumption of unlinked markers: for each clone $j \in B$, the allelic state at marker $\ell$

$$ \alpha_{\ell}^{(i,j \in B)} \sim \text{Uniform}[\rho_{\ell}^{(1)}(B), \rho_{\ell}^{(2)}(B)]. $$

For each clone $j = 1, \dots, M^{(i)}$, we simulate a Bernoulli random variable

$$ \omega^{(i,j)} \stackrel{i.i.d.}{\sim} \text{Bernoulli}(\omega) $$

governing detection due to genotyping at locus $\ell$. We then aggregate the set of detected alleles present at locus $\ell$

$$ G_{\ell}^{(\text{detected})} = \{ \alpha^{(i,j)} | \omega^{(i,j)} = 1 \}. $$

For each detected allele $g \in G_{\ell}^{(\text{detected})}$, we sample an observed allele from the categorical distribution

$$ \gamma(g) \sim \text{Categorical}(\delta_{\ell}(g, 1), \dots, \delta_{\ell}(g, h_{\ell})); $$

classification and population allele frequency estimation is performed based on the observed set of genotypes $\gamma(G_{\ell}^{(\text{detected})})$ for the baseline isolate.

For the $i^{\text{th}}$ baseline isolate of MOI $M^{(i)}$, we generate a set of paired recurrences as follows. For each $M_{\text{recurrent}}^{(i)} = 1, \dots, M_{\max}$ and $m = 0, \dots, \min(M^{(i)}, M_{\text{recurrent}}^{(i)})$, we simulate a mixture of $m$ recrudescent
and ($M_{\text{recurrent}}^{(i)} - m$) newly-inoculated clones. The $m$ recrudescent clones are sampled uniformly at random without replacement from the paired baseline isolate; the ($M_{\text{recurrent}}^{(i)} - m$) newly-inoculated clones are simulated under the same framework as a baseline isolate. We account for the incomplete detection of clones in the recurrent isolate (with the per-clone marker-wise probability of detection $\omega$), but do not account for genotyping error for the recurrent isolate.

Paired simulated recurrences are classified using our novel model *PfRecur*, with a uniform prior over the number of newly-inoculated vs recrudescent clones (i.e., with $\beta = 1$). Classification is based on the maximum observed cardinality across simulated markers, rather than the true MOI of each isolate. We perform classification in an idealised setting with the user-specified parameters $\omega = 0.9$ and $\varepsilon = 0.05$ set to the values under which data were simulated.

Results in the main text are aggregated across 40 simulation runs (i.e., 40 simulated sets of baseline samples), amounting to 28773 simulated recurrent isolates.
# C Empirical data: a re-analysis of Dimbu et al. [17]

## C.1 Within-host diversity

To ascertain the contribution of each microsatellite marker to the MOI of a given isolate, we derive a within-host diversity metric. For an isolate with MOI $M^{(i)}$ and cardinality $M_{\ell}^{(i)}$ at marker $\ell$, we define the within-host diversity of marker $\ell$ to be the expected value of the complement of Nei's gene identity metric (or equivalently, the probability of differing alleles at marker $\ell$ when a pair of clones is sampled with replacement from the isolate), taking a uniform distribution over all possible allelic configurations of $M^{(i)}$ clones that harbour precisely $M_{\ell}^{(i)}$ distinct alleles at locus $\ell$. Following the reasoning in Appendix A.3.1, this takes the form

$$ \xi_{\ell}^{(i)} = 1 - \sum_{k=1}^{M^{(i)} - M_{\ell}^{(i)}} \binom{M^{(i)}}{k} \frac{S(k, 1)S(M^{(i)} - k, M_{\ell}^{(i)} - 1)}{S(M^{(i)}, M_{\ell}^{(i)})} \left(\frac{k}{M^{(i)}}\right)^2 \quad (31) $$

where $S(\cdot, \cdot)$ denotes a Stirling number of the second kind.

## C.2 Sensitivity of PfRecur to user-specified parameters

Here, we consider the sensitivity of our classification model PfRecur to the user-specified parameters $\varepsilon$ (governing the genotyping error probability under the normalised geometric model (7), applied to each allele called in a baseline isolate relative to the recurrent isolate of interest) and $\omega$ (governing the marker-wise probability of detection for each clone in the paired baseline and recurrent isolates of interest). A complete sensitivity analysis is available in an accompanying GitHub repository (https://github.com/somyamehra/PfTreatmentFailure); below, we highlight isolates for which classification is sensitive to user-specified parameters.

### C.2.1 Classification of recurrences in Dimbu et al. [17]

The imperfect detection of clones in the paired baseline and recurrent isolate proves a key consideration for a subset of recurrences in Dimbu et al. [17] using PfRecur. The posterior probability of at least one recrudescent clone (metric M1) is inflated as ungenotyped clones are increasingly accommodated (i.e., the per-clone marker-wise probability of detection $\omega$ is decreased, Figure 7A). However, accommodating undetected clones in the paired baseline isolate protects against plausible human error, such as the incorrect entry of a single allele in a putative recrudescence (Figure 7B). Inferred mixtures of newly-inoculated vs recrudescent clones, quantified through the posterior proportion of recrudescent clones (metric M2), are also sensitive to per-clone marker-wise probability of detection $\omega$ in several cases (Figure 7C). A dependence on the genotyping error probability $\varepsilon$, under the normalised geometric model (7), is also apparent for a subset of isolates (Figure 8); however, our results are largely insensitive to genotyping error probability within the plausible range $0.01 \le \varepsilon \le 0.1$. By default, we set $\varepsilon = 0.05$ and $\omega = 0.9$.
(A) Classification sensitive to per-clone marker-wise probability of detection

(A1)

<table><thead><tr><th></th><th>Baseline<br>(LL21-061D0)<br>Parasitemia: 5.62e+04/µL<br>Temperature: 37.1</th><th>Recurrent<br>(LL21-061D21)<br>Parasitemia: 3.67e+03/µL<br>Temperature: 36.6</th></tr></thead><tbody><tr><td>M313</td><td>242, 248</td><td>242</td></tr><tr><td>M383</td><td>141</td><td>137</td></tr><tr><td>TA1</td><td>165, 168</td><td>174</td></tr><tr><td>POLYA</td><td>153</td><td>147</td></tr><tr><td>PFPK2</td><td>168, 171, 186</td><td>186</td></tr><tr><td>M2490</td><td>81</td><td>81</td></tr><tr><td>TA109</td><td>160, 172, 175</td><td>172</td></tr></tbody></table>

All markers

Genotyping error rate (normalised geometric model): 0.05

<figure><img src="image_8.png" /></figure>

Excluding TA109

Genotyping error rate (normalised geometric model): 0.05

(A2)

<table><thead><tr><th></th><th>Baseline<br>(ZL21-233D0)<br>Parasitemia: 2.06e+03/µL<br>Temperature: 37.5</th><th>Recurrent<br>(ZL21-233D28)<br>Parasitemia: 4.15e+03/µL<br>Temperature: 36.4</th></tr></thead><tbody><tr><td>M313</td><td>234</td><td>236</td></tr><tr><td>M383</td><td>147</td><td>147</td></tr><tr><td>TA1</td><td>177</td><td>168</td></tr><tr><td>POLYA</td><td>141, 144</td><td>144, 162</td></tr><tr><td>PFPK2</td><td>186</td><td>186</td></tr><tr><td>M2490</td><td>72</td><td>81</td></tr><tr><td>TA109</td><td>160, 175</td><td>163</td></tr></tbody></table>

All markers

Genotyping error rate (normalised geometric model): 0.05

<figure><img src="image_9.png" /></figure>

Excluding TA109

Genotyping error rate (normalised geometric model): 0.05

<table><thead><tr><th></th><th>Baseline<br>(LL21-010D0)<br>Parasitemia: 4.79e+04/µL<br>Temperature: 35.4</th><th>Recurrent<br>(LL21-010D14)<br>Parasitemia: 4.75e+01/µL<br>Temperature: 36.7</th></tr></thead><tbody><tr><td>M313</td><td>256</td><td>246, 284</td></tr><tr><td>M383</td><td>123</td><td>123</td></tr><tr><td>TA1</td><td>171</td><td>171</td></tr><tr><td>POLYA</td><td>159</td><td>111, 159</td></tr><tr><td>PFPK2</td><td>183</td><td>183</td></tr><tr><td>M2490</td><td>75</td><td>72</td></tr><tr><td>TA109</td><td>187</td><td>124, 187, 199</td></tr></tbody></table>

All markers

Genotyping error rate (normalised geometric model): 0.05

<figure><img src="image_10.png" /></figure>

Excluding TA109

Genotyping error rate (normalised geometric model): 0.05

<table><thead><tr><th></th><th>Baseline<br>(ZQ21-085D0)<br>Parasitemia: 5.01e+04/µL<br>Temperature: 38.9</th><th>Recurrent<br>(ZQ21-085D28)<br>Parasitemia: 2.88e+04/µL<br>Temperature: 36.5</th></tr></thead><tbody><tr><td>M313</td><td>232, 258</td><td>208, 234</td></tr><tr><td>M383</td><td>125, 137</td><td>125, 139</td></tr><tr><td>TA1</td><td>168</td><td>162, 168</td></tr><tr><td>POLYA</td><td>150, 153</td><td>132, 138, 150</td></tr><tr><td>PFPK2</td><td>171, 183</td><td>171</td></tr><tr><td>M2490</td><td>81</td><td>81</td></tr><tr><td>TA109</td><td>160, 172, 184</td><td>175, 181</td></tr></tbody></table>

All markers

Genotyping error rate (normalised geometric model): 0.05

<figure><img src="image_11.png" /></figure>

Excluding TA109

Genotyping error rate (normalised geometric model): 0.05

<table><thead><tr><th></th><th>Baseline<br>(LL21-077D0)<br>Parasitemia: 2.31e+04/µL<br>Temperature: 39</th><th>Recurrent<br>(LL21-077D14)<br>Parasitemia: 3.90e+01/µL<br>Temperature: 36.1</th></tr></thead><tbody><tr><td>M313</td><td>224, 236</td><td>236</td></tr><tr><td>M383</td><td>123, 143</td><td>143</td></tr><tr><td>TA1</td><td>165</td><td></td></tr><tr><td>POLYA</td><td>153, 183</td><td>153</td></tr><tr><td>PFPK2</td><td>165, 186</td><td>186</td></tr><tr><td>M2490</td><td>81</td><td>78</td></tr><tr><td>TA109</td><td>163, 172</td><td>109, 124, 163, 184, 199</td></tr></tbody></table>

All markers

Genotyping error rate (normalised geometric model): 0.05

<figure><img src="image_12.png" /></figure>

Excluding TA109

Genotyping error rate (normalised geometric model): 0.05

<table><thead><tr><th></th><th>Baseline<br>(BD21-053D0)<br>Parasitemia: 9.79e+04/µL<br>Temperature: 39.6</th><th>Recurrent<br>(BD21-053D35)<br>Parasitemia: 1.76e+03/µL<br>Temperature: 36.1</th></tr></thead><tbody><tr><td>M313</td><td>210</td><td>218</td></tr><tr><td>M383</td><td>123, 129</td><td>129</td></tr><tr><td>TA1</td><td>165, 168, 177, 183</td><td></td></tr><tr><td>POLYA</td><td>153, 165</td><td>129, 138</td></tr><tr><td>PFPK2</td><td>162, 168, 180</td><td>168</td></tr><tr><td>M2490</td><td>75, 78, 81</td><td></td></tr><tr><td>TA109</td><td>160</td><td>124, 184, 196, 199</td></tr></tbody></table>

All markers

Genotyping error rate (normalised geometric model): 0.05

<figure><img src="image_13.png" /></figure>

Excluding TA109

Genotyping error rate (normalised geometric model): 0.05

<caption>**Figure 7:** Classifications under *PfRecur* sensitive to the per-clone marker-wise probability of detection $\omega$ in the paired baseline/recurrent isolates.</caption>
## (B) Accommodating imperfect detection protects against possible human error

(B1)

<table><thead><tr><th></th><th>Baseline<br>(ZQ21-103D0)<br>Parasitemia: 1.59e+05/µL<br>Temperature: 38.8</th><th>Recurrent<br>(ZQ21-103D18)<br>Parasitemia: 1.89e+05/µL<br>Temperature: 36</th></tr></thead><tbody><tr><td>M313</td><td>226, 238, 256</td><td>238</td></tr><tr><td>M383</td><td>133, 139</td><td>139</td></tr><tr><td>TA1</td><td>159, 165</td><td>195</td></tr><tr><td>POLYA</td><td>150, 156</td><td>150</td></tr><tr><td>PFPK2</td><td>168, 171, 186, 189, 192</td><td>171</td></tr><tr><td>M2490</td><td>72, 78, 81</td><td>78</td></tr><tr><td>TA109</td><td>160, 163, 172, 175, 184</td><td>184</td></tr></tbody></table>

All markers

Genotyping error rate (normalised geometric model): 0.05

Excluding TA109

Genotyping error rate (normalised geometric model): 0.05

<figure><img src="image_14.png" /></figure>

(B2)

<table><thead><tr><th></th><th>Baseline<br>(ZL21-245D0)<br>Parasitemia: 1.58e+05/µL<br>Temperature: 37.9</th><th>Recurrent<br>(ZL21-245D21)<br>Parasitemia: 2.20e+04/µL<br>Temperature: 35.6</th></tr></thead><tbody><tr><td>M313</td><td>228, 234</td><td>234</td></tr><tr><td>M383</td><td>125, 131</td><td>151</td></tr><tr><td>TA1</td><td>162, 165, 183</td><td>165</td></tr><tr><td>POLYA</td><td>156, 159, 177</td><td>159</td></tr><tr><td>PFPK2</td><td>159, 165</td><td>159</td></tr><tr><td>M2490</td><td>78, 81</td><td>81</td></tr><tr><td>TA109</td><td>163, 175, 184</td><td>163</td></tr></tbody></table>

All markers

Genotyping error rate (normalised geometric model): 0.05

Excluding TA109

Genotyping error rate (normalised geometric model): 0.05

<figure><img src="image_15.png" /></figure>

## (C) Resolution of mixtures sensitive to per-clone marker-wise probability of detection

(C1)

<table><thead><tr><th></th><th>Baseline<br>(ZL21-203D0)<br>Parasitemia: 1.32e+05/µL<br>Temperature: 37.7</th><th>Recurrent<br>(ZL21-203D14)<br>Parasitemia: 3.73e+02/µL<br>Temperature: 35.8</th></tr></thead><tbody><tr><td>M313</td><td>240, 246</td><td>240</td></tr><tr><td>M383</td><td>173</td><td>159, 173</td></tr><tr><td>TA1</td><td>168, 171</td><td>168</td></tr><tr><td>POLYA</td><td>159, 180</td><td>180</td></tr><tr><td>PFPK2</td><td>165, 171</td><td>171</td></tr><tr><td>M2490</td><td>81</td><td>81</td></tr><tr><td>TA109</td><td>172, 178</td><td>172, 178, 199</td></tr></tbody></table>

All markers

Genotyping error rate (normalised geometric model): 0.05

Excluding TA109

Genotyping error rate (normalised geometric model): 0.05

<figure><img src="image_16.png" /></figure>

(C2)

<table><thead><tr><th></th><th>Baseline<br>(ZL21-269D0)<br>Parasitemia: 5.41e+04/µL<br>Temperature: 38.9</th><th>Recurrent<br>(ZL21-269D21)<br>Parasitemia: 1.94e+04/µL<br>Temperature: 37.2</th></tr></thead><tbody><tr><td>M313</td><td>250</td><td>210, 250</td></tr><tr><td>M383</td><td>135</td><td>135</td></tr><tr><td>TA1</td><td>168</td><td>168</td></tr><tr><td>POLYA</td><td>159</td><td>150, 159</td></tr><tr><td>PFPK2</td><td>171, 183</td><td>171, 183</td></tr><tr><td>M2490</td><td>78, 81</td><td>78, 81</td></tr><tr><td>TA109</td><td>172</td><td>172</td></tr></tbody></table>

All markers

Genotyping error rate (normalised geometric model): 0.05

Excluding TA109

Genotyping error rate (normalised geometric model): 0.05

<figure><img src="image_17.png" /></figure>

(C3)

<table><thead><tr><th></th><th>Baseline<br>(ZL21-304D0)<br>Parasitemia: 1.51e+04/µL<br>Temperature: 37</th><th>Recurrent<br>(ZL21-304D14)<br>Parasitemia: 3.67e+02/µL<br>Temperature: 36.8</th></tr></thead><tbody><tr><td>M313</td><td>218, 232</td><td>218</td></tr><tr><td>M383</td><td>141, 151</td><td>139, 151</td></tr><tr><td>TA1</td><td>138, 168</td><td>138, 168</td></tr><tr><td>POLYA</td><td>153, 156, 162</td><td>156, 162</td></tr><tr><td>PFPK2</td><td>159, 162, 177, 186</td><td>159, 177</td></tr><tr><td>M2490</td><td>78, 81, 84</td><td>78, 81</td></tr><tr><td>TA109</td><td>163, 175</td><td>175, 199</td></tr></tbody></table>

All markers

Genotyping error rate (normalised geometric model): 0.05

Excluding TA109

Genotyping error rate (normalised geometric model): 0.05

<figure><img src="image_18.png" /></figure>

(C4)

<table><thead><tr><th></th><th>Baseline<br>(LL21-054D0)<br>Parasitemia: 2.04e+04/µL<br>Temperature: 37.9</th><th>Recurrent<br>(LL21-054D14)<br>Parasitemia: 6.30e+01/µL<br>Temperature: 38</th></tr></thead><tbody><tr><td>M313</td><td>264</td><td>264</td></tr><tr><td>M383</td><td>131</td><td>131</td></tr><tr><td>TA1</td><td>171</td><td>171</td></tr><tr><td>POLYA</td><td>159</td><td>159</td></tr><tr><td>PFPK2</td><td>183, 186</td><td>186</td></tr><tr><td>M2490</td><td>84</td><td></td></tr><tr><td>TA109</td><td>196</td><td>124, 199, 223</td></tr></tbody></table>

All markers

Genotyping error rate (normalised geometric model): 0.05

Excluding TA109

Genotyping error rate (normalised geometric model): 0.05

<figure><img src="image_19.png" /><figcaption>**Figure 7:** Classifications under *PfRecur* sensitive to the per-clone marker-wise probability of detection $\omega$ in the paired baseline/recurrent isolates.</figcaption></figure>
Classification sensitive to genotyping error

<table><caption>(A1)</caption><thead><tr><th></th><th>Baseline<br>(ZL21-2600)<br>Parasitemia: 2.29e+03/µL<br>Temperature: 37.2</th><th>Recurrent<br>(ZL21-260D14)<br>Parasitemia: 2.82e+02/µL<br>Temperature: 36</th></tr></thead><tbody><tr><td>M313</td><td>218, 226, 238</td><td>226, 238</td></tr><tr><td>M383</td><td>123, 137</td><td>123</td></tr><tr><td>TA1</td><td>159, 162, 165, 168, 171</td><td>141, 174</td></tr><tr><td>POLYA</td><td>150, 153, 171</td><td>153</td></tr><tr><td>PFPK2</td><td>165, 168, 189</td><td>165, 171</td></tr><tr><td>M2490</td><td>81, 84</td><td>81</td></tr><tr><td>TA109</td><td>175, 184, 199</td><td>175, 196, 199</td></tr></tbody></table>

<table><caption>(A2)</caption><thead><tr><th></th><th>Baseline<br>(ZQ21-085D0)<br>Parasitemia: 5.01e+04/µL<br>Temperature: 38.9</th><th>Recurrent<br>(ZQ21-085D28)<br>Parasitemia: 2.88e+04/µL<br>Temperature: 36.5</th></tr></thead><tbody><tr><td>M313</td><td>232, 258</td><td>208, 234</td></tr><tr><td>M383</td><td>125, 137</td><td>125, 139</td></tr><tr><td>TA1</td><td>168</td><td>162, 168</td></tr><tr><td>POLYA</td><td>150, 153</td><td>132, 138, 150</td></tr><tr><td>PFPK2</td><td>171, 183</td><td>171</td></tr><tr><td>M2490</td><td>81</td><td>81</td></tr><tr><td>TA109</td><td>160, 172, 184</td><td>175, 181</td></tr></tbody></table>

<table><caption>(A3)</caption><thead><tr><th></th><th>Baseline<br>(ZL21-233D0)<br>Parasitemia: 2.06e+03/µL<br>Temperature: 37.5</th><th>Recurrent<br>(ZL21-233D28)<br>Parasitemia: 4.15e+03/µL<br>Temperature: 36.4</th></tr></thead><tbody><tr><td>M313</td><td>234</td><td>236</td></tr><tr><td>M383</td><td>147</td><td>147</td></tr><tr><td>TA1</td><td>177</td><td>168</td></tr><tr><td>POLYA</td><td>141, 144</td><td>144, 162</td></tr><tr><td>PFPK2</td><td>186</td><td>186</td></tr><tr><td>M2490</td><td>72</td><td>81</td></tr><tr><td>TA109</td><td>160, 175</td><td>163</td></tr></tbody></table>

<table><caption>(A4)</caption><thead><tr><th></th><th>Baseline<br>(LL21-061D0)<br>Parasitemia: 5.62e+04/µL<br>Temperature: 37.1</th><th>Recurrent<br>(LL21-061D21)<br>Parasitemia: 3.67e+03/µL<br>Temperature: 36.6</th></tr></thead><tbody><tr><td>M313</td><td>242, 248</td><td>242</td></tr><tr><td>M383</td><td>141</td><td>137</td></tr><tr><td>TA1</td><td>165, 168</td><td>174</td></tr><tr><td>POLYA</td><td>153</td><td>147</td></tr><tr><td>PFPK2</td><td>168, 171, 186</td><td>186</td></tr><tr><td>M2490</td><td>81</td><td>81</td></tr><tr><td>TA109</td><td>160, 172, 175</td><td>172</td></tr></tbody></table>

<table><caption>(A5)</caption><thead><tr><th></th><th>Baseline<br>(BD21-053D0)<br>Parasitemia: 9.79e+04/µL<br>Temperature: 39.6</th><th>Recurrent<br>(BD21-053D35)<br>Parasitemia: 1.76e+03/µL<br>Temperature: 36.1</th></tr></thead><tbody><tr><td>M313</td><td>210</td><td>218</td></tr><tr><td>M383</td><td>123, 129</td><td>129</td></tr><tr><td>TA1</td><td>165, 168, 177, 183</td><td>129, 138</td></tr><tr><td>POLYA</td><td>153, 165</td><td>168</td></tr><tr><td>PFPK2</td><td>162, 168, 180</td><td>168</td></tr><tr><td>M2490</td><td>75, 78, 81</td><td></td></tr><tr><td>TA109</td><td>160</td><td>124, 184, 196, 199</td></tr></tbody></table>

<figure><img src="image_20.png" /><figcaption>**Figure 8:** Classifications under *PfRecur* sensitive to the genotyping error probability $\varepsilon$ (applied to each allele called in a baseline isolate relative to the recurrent isolate of interest) under the normalised geometric model (7).</figcaption></figure>
### C.2.2 False positive recrudescence rates

Relaxing the per-clone marker-wise probability of detection $\omega$ can amplify the posterior probability of recrudescence under *PfRecur*. However, false positive recrudescence rates across 500 permuted artificial 'not-recrudescence' datasets (obtained by generating derangements of baseline study participant labels within each study site in [17]) remain consistently low for $\varepsilon = 0.05$ and $\omega \in [0.75, 1]$. Empirical complementary CDFs for metric the posterior probability of recrudescence (under the CDC model) vs metric M1 (under *PfRecur*) using which false positive rates are computed are shown in Figure 10.

<figure><img src="image_21.png" /></figure>

**Figure 9:** False positive recrudescence rates across 500 permuted artificial 'not-recrudescence' datasets under *PfRecur* (metric M1) with genotyping error probability $\varepsilon = 0.05$ and per-clone marker-wise detection probability $\omega \in [0.75, 1]$ vs the CDC model vs a match-counting approach. Points show medians; bold error bars show the interquartile range; grey error bars show 95% confidence intervals.

<figure><img src="image_22.png" /></figure>

**Figure 10:** Complementary empirical CDFs for the posterior probability of recrudescence (CDC model) vs metric M1 (*PfRecur*) aggregated over 500 permuted artificial 'not-recrudescence' datasets.

### C.3 Comparison against CDC model

We find that the CDC model [8] typically yields elevated posterior probabilities of recrudescence for multiclonal isolates with matching alleles at a subset of markers. Recurrences with differing classifications under *PfRecur* vs the CDC model [8] are summarised in Figure 11.
<figure><img src="image_23.png" /><figcaption>**Figure 11:** Differing classifications under *PfRecur* (metrics M1 and M2) vs the CDC model [8]. For *PfRecur*, we take a genotyping error probability of $\varepsilon = 0.05$ under the normalised geometric model (29) (applied to each allele called in each baseline isolate), and assume that the per-clone marker-wise probability of detection (applied to clones in the paired baseline/recurrent isolates) is $\omega = 0.9$.</figcaption></figure>
### C.3.1 Efficacy estimates

To generate efficacy estimates adjusted for classification uncertainty and right censoring, we perform the Kaplan-Meier correction used in Dimbu et al. [17]<fnref n="3" />: for each study arm, we perform 10,000 iterations of bootstrap resampling whereby we sample a binary reinfection or recrudescence state for each patient under a Bernoulli distribution with success parameter given by the posterior probability of recrudescence (or metric M1 under *PfRecur*); and then compute Kaplan-Meier estimates (including the lower 2.5% and upper 97.5% confidence interval under the *peto* method implemented in the R function *survival::survfit* [45]) for these binary infection states. Survival estimates for the last day of follow-up (28 days for AL and ASAQ; 42 days for DP and ASPY) are then averaged across these 10,000 iterations to yield an adjusted efficacy estimate and an accompanying 95% confidence interval. Efficacy estimates under *PfRecur* vs the CDC model (re-run on our system) are summarised in Table 4. Efficacy estimates generated using the CDC model in Table 4 differ from those reported in Dimbu et al. [17] due to the stochastic nature of the MCMC algorithm, and because early treatment failures (1 patient treated with AL, and 3 patients treated with ASAQ in Zaire) have not been taken into account in our re-analysis.

<table><caption>**Table 4:** Efficacy estimates without PCR correction (as reported in Table 4 of [17]) vs based on posterior probabilities of recrudescence under the CDC model vs metric M1 under *PfRecur* (with default parameters $\omega = 0.9, \varepsilon = 0.05$), adjusted for loss to follow-up under a Kaplan-Meier model [8, 17].</caption><thead><tr><th></th><th>Study site</th><th>Uncorrected</th><th>CDC model</th><th>PfRecur</th></tr></thead><tbody><tr><td rowspan="2">AL (day 28)</td><td>Zaire</td><td>77.1 (69–86)</td><td>88.9 (82.8–95.3)</td><td>90.4 (84.5–96.4)</td></tr><tr><td>Lunda Sul</td><td>87.3 (81–94)</td><td>94.4 (89.9–99.0)</td><td>96.3 (92.6–100)</td></tr><tr><td rowspan="2">ASAQ (day 28)</td><td>Zaire</td><td>82.9 (76–91)</td><td>93.9 (89.1–98.9)</td><td>94.4 (89.8–99.2)</td></tr><tr><td>Lunda Sul</td><td>100</td><td>100</td><td>100</td></tr><tr><td>SPY (day 42)</td><td>Benguela</td><td>87.2 (81–94)</td><td>99.7 (99.2–100)</td><td>100 (100–100)</td></tr><tr><td>DP (day 42)</td><td>Benguela</td><td>94.2 (90–99)</td><td>98.3 (95.8–100)</td><td>99.0 (97.1–100)</td></tr></tbody></table>

<footnote n="3">³Implemented in https://github.com/MateuszPlucinski/AngolaTES2021 [40]</footnote>
## D Comparison against CDC model structure [8]

An extended comparison of the *PfRecur* model structure against the CDC model [8] is provided below. In brief, the CDC model [8] explicitly estimates allelic configurations using a Gibbs sampler; averages the likelihood of recrudescence across pairs of clones in a baseline and recurrent isolate; and concurrently generates estimates of parametric (geometric) genotyping error probabilities and a per-clone unobserved allele “penalty”. In contrast, in *PfRecur* we average analytically over compatible allelic configurations in each isolate to obtain an analytically tractable posterior; model recurrent isolates as mixtures of newly-inoculated vs recrudescent clones; and treat both genotyping error and the per-clone marker-wise probability of detection for paired baseline/recurrent isolates as user-specified parameters.

<table><thead><tr><td>Feature</td><td>CDC model [8]</td><td>PfRecur</td></tr></thead><tbody><tr><td>Model structure</td><td>Construct joint likelihood for baseline and recurrent isolates, allowing the concurrent estimation of allele frequencies and genotyping error</td><td>Construct likelihood for observed genotypes in recurrent isolates only</td></tr><tr><td>Likelihood of reinfection</td><td>Allele-frequency based likelihood predicated on population allele frequencies estimated over all available isolates, excluding recrudescent alleles</td><td>Allele-frequency based likelihood for newly-inoculated clones predicated on population allele frequencies derived from baseline isolates only (excluding the paired baseline isolate)</td></tr><tr><td>Likelihood of recrudescence</td><td>For recrudescent alleles, replace allele frequencies with a genotyping error term (probability that the allele is identical to a baseline allele)</td><td>Allele-frequency based likelihood for recrudescent clones predicated on allele frequencies largely derived from the paired baseline isolate, but with allelic states for ungenotyped clones in the paired baseline isolate imputed based on allele frequencies for newly-inoculated clones</td></tr><tr><td>Adjustment for multiple comparisons</td><td>Average likelihood of recrudescence over pairwise comparisons of baseline/recurrent clones at each locus</td><td>Likelihood averaged over all compatible allelic configurations in the recurrent isolate (assume a uniform distribution across compatible allelic configurations); derived allele frequencies adjusted for allelic configurations in baseline isolates</td></tr></tbody></table>
<table><tr><td>Population allele frequencies</td><td>Estimated under a multinomial-Dirichlet model based on sample allele frequencies pooled across all available baseline and recurrent isolates, including unobserved alleles, but without double-counting recrudescent alleles (these are defined to be the pair of alleles at each marker in a baseline/recurrent pair that have the smallest length difference; may be repeated in a given isolate); initialise sample allele frequencies by summing over the observed set of alleles for each isolate (i.e., weighting each isolate by its marker-wise cardinality)</td><td>Allele frequencies for newly-inoculated clones derived under a multinomial-Dirichlet model based on allele counts across baseline isolates excluding the paired baseline isolate, averaging over compatible allelic configurations in each relevant baseline isolate (assuming perfect detection of clones) and adjusted for genotyping error; allele frequencies for recrudescent clones derived over clones in the paired baseline isolate, with allelic states for undetected clones in the paired baseline isolate imputed using allele frequencies for newly-inoculated clones</td></tr><tr><td>Genotyping error structure</td><td>Parametric (approximately geometric with respect to allele lengths); estimated under a beta-beta model in the Gibbs sampler, based on the minimal distance between the closest alleles at each locus in recurrences deemed to be recrudescences</td><td>Non-parametric; user-specified; formulated relative to the recurrent isolate and applied to all baseline isolates (as the probability that an allele called in the baseline isolate matches an allele called in the recurrent isolate); per-isolate allele counts (used to derive allele frequencies) are adjusted for genotyping error</td></tr><tr><td>Penalising unobserved alleles</td><td>Number of hidden alleles across all isolates defined to be the difference between MOI and the locus cardinality; implement a per-clone unobserved allele penalty q whenever a hidden allele lies outside the set of observed alleles; q estimated under a beta-binomial model in the Gibbs sampler, based on the number of hidden alleles (fixed across iterations)</td><td>Implement a marker-wise probability ω of detecting each clone in the paired baseline and recurrent isolates; ω is user-specified</td></tr><tr><td>Adjusting for unobserved alleles</td><td>Update one randomly chosen hidden allele in a baseline/recurrent pair per iteration; sample hidden alleles uniformly at random over the set of possible alleles at that marker; likelihood ratio for Gibbs sampler based on the convolution over frequency of each allele (true underlying allele) and genotyping error (probability true allele appears as sampled hidden allele), multiplied by the per-clone unobserved allele penalty q if relevant; if a switch is to be made, add a length variability term to sampled hidden allele (adjusted for allelic binning: based on the mean SD within allelic classes at the marker)</td><td>Impute allelic states for undetected clones in the paired baseline isolate based on derived population allele frequencies for newly-inoculated clones, with a marker-wise truncated binomial model for the number of clones detected at each marker; derive the locus-wise likelihood of observed genotypes for the recurrent isolate over detected clones only, with a truncated multinomial model for the number of newly-inoculated vs recrudescent clones detected at each marker</td></tr></table>