Malaria and similar infectious diseases often exhibit a complex landscape of pathogen genetic diversity within individual infections. In malaria, the presence of multiple genetically distinct pathogen variants within a single infection is typically referred to as complexity of infection (COI) or multiplicity of infection (MOI) (Read and Taylor, 2001; Chang et al., 2017; Schneider, 2018). MOI is commonly reported in the context of malaria molecular surveillance, which has become increasingly popular in the last years as molecular assays become more affordable in endemic settings and due to WHO recommendations (World Health Organization, 2022; Sinha et al., 2023).

Epidemiologically, MOI is an important parameter as it scales with transmission intensities, however not necessarily in a linear way (Pacheco et al., 2020; Sinha et al., 2023). The effect of interactions of different pathogen variants within infections (intra-host dynamics) has been considered important in several theoretical models (e.g., Hastings and Watkins, 2005; Gurarie and McKenzie, 2006). However, so far, empirical evidence on the impact of MOI on the clinical pathogenesis of malaria remains inconclusive (Pacheco et al., 2016). Nevertheless, the distribution of MOI influences the evolutionary dynamics of malaria (Schneider, 2021; Schneider and Salas, 2022) and mediates evolutionary–genetic patterns and pathogen genetic diversity (e.g., patterns of genetic hitchhiking and linkage disequilibrium) as it affects the effective rate of recombination as illustrated in Figure 1 (Schneider and Kim, 2010; Alizon et al., 2013). Moreover, there is an important link between the frequency distribution of pathogen variants, their occurrence within infections (i.e., prevalence), and MOI. Specifically, given the frequency of a certain pathogen variant, its prevalence increases with MOI. This is particularly relevant in the context of anti-malarial drug resistance and seasonal malaria with varying transmission intensities (Geiger et al., 2014; Schneider, 2021; Schneider et al., 2022).

Importantly, MOI or COI is not unambiguously defined. Several formulations of MOI exist, with discrepancies between verbal and formal definitions (which typically underlay statistical models) as discussed in detail in Schneider et al. (2022). Here we follow the suggested definition in Schneider et al. (2022), which is used in most theoretical frameworks. Particularly, MOI is not defined as the number of distinct pathogen variants within an infection but as the number of “super-infections” during one disease episode. More precisely, MOI is the number of independent infectious events, assuming that exactly one pathogen variant (lineage) is transmitted per event (Hill and Babiker, 1995; Schneider and Escalante, 2014). This implies that MOI is an unobservable quantity because a host can be infected multiple times with the same pathogen variant ( Figure 1 ), and these infectious events cannot be reconstructed from molecular assays ( Figure 2 ). Notably, as pointed out in Schneider et al. (2022), if the distribution of MOI in the population is known or estimated, the distribution of different pathogen variants within infections can be derived (but not vice versa). The definition of MOI used here only approximately accounts for “co-infections”, i.e., the co-transmission of several pathogenic variants during an infective episode (Schneider, 2021). Focusing on super-infections (and only approximately covering co-infections) has the pragmatic advantage that no explicit model of vector–host transmission has to be specified (Schneider, 2021) considering that a co-infection is particularly relevant if one aims to determine the genetic relatedness of pathogen variants within an infection. To obtain enough resolution to study genetic relatedness, these approaches require high-quality genomic data. Although such approaches have become increasingly popular in the last years (cf. Nkhoma et al., 2012; Wong et al., 2018; Zhu et al., 2019; Nkhoma et al., 2020; Dia and Cheeseman, 2021; Neafsey et al., 2021), these approaches are inappropriate if only a handful of genetic markers are available.

Heuristic methods to estimate the distribution of MOI are typically biased; methods based on a solid statistical framework are preferable (Schneider, 2021). Several such methods (which are essentially based on the same statistical framework) have been proposed to estimate the MOI and lineage frequencies for various assumptions concerning the genetic architecture of the underlying molecular data. All these methods make certain assumptions on the distribution of MOI in the population. The maximum likelihood-based method of Hill and Babiker (1995) assumes that MOI follows either a (conditional) Poisson or negative binomial distribution and is based on one or two genetic markers. In the case of MOI following a positive Poisson distribution, this method was refined by applying a bias correction in the case of a single molecular marker by Hashemi and Schneider (2021) and to an arbitrary number of biallelic molecular markers (Tsoungui Obama and Schneider, 2022). Moreover, the approach of Li et al. (2007) assumes that MOI follows a (conditional) Poisson distribution. The Bayesian method in the program COIL (Galinsky et al., 2015) and its generalization THE REAL McCOIL (Chang et al., 2017) does not explicitly specify a specific distribution for MOI per se, but the implementation imposes a uniform distribution, which constrains the resulting posterior distribution of MOI.

Under the assumption that infective mosquito bites are rare and independent in a population with homogeneous exposure, MOI follows a Poisson distribution. This renders this parametric choice an important null model. However, if mosquito biting rates are heterogeneous in the population, the distribution of MOI will more likely follow a mixture of Poisson or a negative binomial distribution (Schneider et al., 2022). In fact, empirical evidence indicates that mosquito biting patterns are heterogeneous, with certain individuals experiencing more bites than others. This is influenced by confounding factors such as environmental conditions or individual attractiveness to mosquitoes (Noor et al., 2014; Guelbeogo et al., 2018). Consequently, the number of mosquito bites per person tends to be over-dispersed compared to the Poisson distribution (Irvine et al., 2018). (Note, however, that an over-dispersed mosquito biting rate does not imply that the MOI distribution is over-dispersed, as it is concerned only with the number of infective bites.) In any case, significant deviations from the Poisson assumptions suggest that the negative binomial distribution might be more suitable (Lloyd-Smith, 2007). However, maximum likelihood estimation of the negative binomial distribution in general is problematic. As mentioned in Adamidis (1999), if the empirical variance is not larger than the mean of count data, the maximum likelihood estimates of the parameters of the negative binomial distribution are degenerate. A way to overcome the problem that estimation of both parameters can compromise the stability and interpretability of a negative binomial model (Bandara et al., 2019) is to estimate only one parameter, while fixing the other, as suggested in Piegorsch (1990), Saha and Paul (2005), and Lloyd-Smith (2007). In practice, this means that there needs to be prior information on one parameter. However, this is typically not feasible in the context of MOI because it would require additional information such as mosquito biting rates or host exposure to mosquitoes. Such information is typically beyond the scope of malaria molecular surveillance.

If there is a prior belief that MOI does not follow a Poisson distribution, rather than assuming that MOI falls into a different class of parametric distributions, such as the negative binomial distribution, no particular class of distributions has to be imposed. Such a non-parametric approach offers a valid alternative if the MOI is completely unknown because of its flexibility. In fact, a non-parametric approach is the most flexible approach in this context.

Here we introduce a non-parametric statistical model to estimate the distribution of MOI and pathogen lineage (allele) frequencies from a single molecular marker by maximum likelihood. Non-parametric refers to the fact that the MOI distribution is not assumed to fall into a class of parametric distributions. The statistical model is first introduced in “Materials and methods”. Because the resulting likelihood function is too complex to have a closed-form solution, we derive the expectation maximization algorithm (EM algorithm) to derive the maximum likelihood estimate (MLE) numerically (Couvreur, 1997; Ng et al., 2012). The EM algorithm provides a numerically stable iteration to derive the MLE. By numerical simulations, we further investigate the performance of the non-parametric estimator in terms of bias and variance and contrast it to MOI estimates based on the assumption of an underlying Poisson distribution. The method proposed here is further applied to three data sets from Cameroon, Kenya, and Venezuela as an illustration. The method is implemented as an R script available in the Supplementary Material and at https://github.com/Maths-against-Malaria/Non-parametric-MOI-estimation .

In the following mathematical notation, we use oblique letters, e.g., m, p, x , to indicate vectors, and italic fonts, e.g., m, p, x, to refer to integers or scalars.

We consider a pathogen population with lineages A ₁ , …, A_n , detected at a single marker locus. Each lineage (or allele) A_i has relative frequency p_i in the pathogen population, jointly denoted by the vector p = (p ₁, …, p_n ). We assume that at each infective event, the mosquito vector transmits exactly one lineage to the host. This corresponds to randomly sampling one lineage from the pathogen population. Hence, co-infections (cf. Figure 1 in Schneider et al., 2022), i.e., the co-transmission of several distinct lineages during one infective event, are ignored. However, hosts can be infected multiple times by different mosquitoes (super-infections) during one disease episode. It is assumed that super-infections occur during relatively short time periods, e.g., a few days, so that all infecting variants reach detectable concentrations. Infective events, in which the variants do not reach detectable frequencies, do not count as super-infection as these (i) are irrelevant for the clinical pathogenesis of the disease and (ii) are undetectable by molecular assays. Following Schneider et al. (2022), we refer to the number of infective events during one disease episode as multiplicity of infection (MOI). Importantly, the hosts might be infected multiple times with the same lineage. Formally, if m_i represents the number of times an individual was infected with lineage A_i , then m_i = 0 if the host was not infected with lineage A_i . Summing over all m_i yields MOI m, i.e., MOI m is defined by:

where m = (m ₁, …, m_n ). The m lineages infecting a host are randomly sampled (with replacement) from the pathogen population. Therefore, within an infection, the configuration of pathogen lineages m follows a multinomial distribution with parameters m = | m | and p , i.e., m ∼ Multi (m,p ). Hence, given that a host has MOI m, the probability of configuration m , i.e., of being infected m_i times with lineage A_i (i = 1,…, n), is:

where ( m m ) : = m ! m 1 ! … m n ! in Equation 1 is the multinomial coefficient and p m : = p 1 m 1 … p n m n .

The configuration of infecting pathogen lineages ( m ) and even MOI (m = | m |) is unobservable. Specifically, from a blood sample of an infected person, the observation is limited to the absence/presence of the infecting lineages ( Figure 2 ). Formally, we represent an observation as the vector x = (x_i ), i = 1,…, n, such that the entries x_i are 0 or 1 (formally, x_i ∈ {0,1}), where 0 denotes the absence and 1 denotes the presence of the lineages in the infection ( Figure 2 ). Thus, x_i is the sign of m_i , i.e., x_i = 0 if m_i = 0 and x_i = 1 if m_i ≥ 1 (formally, x_i := sign m_i , or in compact notation x = sign m ).

Because we consider only disease-positive samples, the observations x correspond to vectors of length n with entries 0 and 1 excluding the vector that contains only zeros, 0, which corresponds to a disease-negative sample. In mathematical notation, x are elements of the set 𝒪 := {0,1} ⁿ \ {0}.

The underlying assumption is that molecular/genetic methods are not quantifying the concentration of lineages but rather detect their presence. Here it is ignored that lineages remain undetected. While this can be included in a statistical model (see, e.g., Hashemi and Schneider, 2024), here it is ignored. For more discussion on undetected or erroneously detected variants, see Schneider et al. (2022).

To express the probability of x , further notation is needed. We call an observation y ∈ 𝒪 a sub-observation of x (denoted y ≼ x ); all lineages observed in y are also observed in x , i.e., if y_i ≤ x_i for i = 1,…, n. The set of all sub-observations of x is denoted by:

We define κ_m := P[MOI = m] as the probability that a host is infected exactly m times (MOI = m) and collectively denote the MOI distribution by κ : = (κ ₁ ,κ ₂,…). Furthermore, the probability generating function (PGF) of the MOI distribution evaluated at a point z is denoted by G(z) (see “Probability distribution of observations” in the Appendix ); the probability of observation x is derived to be as:

where we jointly denote the model parameters by θ = ( κ , p ) and P [ x | θ ] for P [ x ] to emphasize the dependency on the model parameters whenever necessary. Clearly, the probability in (Equation 2) depends on the model parameters p and κ (through the PGF).

Note that while p_i is the frequency of lineage A_i , as shown in Schneider et al. (2022), its prevalence, i.e., the probability that this lineage occurs in an infection, is given by

i.e., the PGF of the MOI distribution links frequency and prevalence.

The model parameters θ , i.e., the distribution of MOI and the lineage frequency distribution, can be estimated from the probabilistic model (Equation 2). We proceed with maximum likelihood (ML) estimation.

First, it is shown how the maximum likelihood estimate (MLE) is derived. This is followed by results on the estimator’s performance in terms of bias and variance.

Estimating the multiplicity of infection (MOI) became popular in malaria molecular surveillance. Ad hoc methods are generally biased and have undesirable statistical properties. Methods based on probabilistic models are often based on similar assumptions. A popular assumption in many models is that MOI follows a Poisson distribution. This corresponds to the assumption that infectious events are rare and independent. (Remember, MOI is defined here as the number of super-infections within one disease episode.) Because there is empirical evidence that mosquito biting rates are over-dispersed, the Poisson assumption is challenged.

Here we explored the most flexible alternative, namely, the situation where no parametric assumption about the distribution of MOI is made (except that a maximum MOI exists). Note that the statistical model itself makes a number of simplifying assumptions. First, the duration of how long individuals are infected is not considered. Rather, it is assumed that super-infections occur within a short period of time, such that all variants which were successfully transmitted to a host, reach sufficient frequencies in the infection to be detectable by PCR. Thus, super-infections which do not contribute to the clinical pathogenesis of an infection are ignored. Furthermore, it is assumed that molecular assays have perfect sensitivity and specificity for all relevant lineages in an infection. Notably, the statistical model can be extended to account for missing information (imperfect sensitivity) as in Hashemi and Schneider (2024) (for the Poisson model); however, if all samples contain molecular information (i.e., no sample with missing data), the extended model reduces to the present one. More alternatives can be found in Okell et al. (2017). Extending the model to include imperfect specificity, i.e., including erroneously detected variants, is more challenging. A model that yields appropriate estimates is the one given by Plucinski et al. (2015). However, it is designed to distinguish recrudescence of reinfections and hence requires paired samples (i.e., two or more sample points for at least some patients). Furthermore, it is assumed that the molecular assays used can only provide absence/presence data but cannot quantify the relative abundance of lineages within an infection.

Given this framework, the resulting non-parametric model is more complicated than the corresponding Poisson model (cf. Schneider and Escalante, 2014), which falls into the class of exponential families (Hashemi and Schneider, 2024). This implies the usual desirable properties of maximum likelihood estimators for the Poisson model (existence and uniqueness of the MLE, efficiency, and consistency, cf. Hashemi and Schneider, 2024). Unfortunately, the non-parametric model is no longer within an exponential family and there is no proof for the same desirable theoretical properties. However, our numerical investigations suggest that they hold. In the case of the Poisson model, the MLE has an intuitive interpretation. Namely, the MLE is the parameter choice for which the empirical prevalences coincide with the expected prevalences. Moreover, the empirical prevalences form a sufficient statistic for the parameter estimation. For the non-parametric model, this interpretation is lost and the prevalences no longer form a sufficient statistic. In fact, the model used more information, namely not just how often lineages occur in infections but also in which configuration they are observed. As for the Poisson model, there exists no closed solution for the MLE of the non-parametric model, and it has to be derived numerically. For this purpose, the EM algorithm was employed here. It is a numerically stable and efficient algorithm to derive the MLE.

Notably, the non-parametric model also has certain advantages compared to the Poisson model. If the data lies on the boundary of the closed convex hull of the admissible sample space (when rewritten as a natural exponential family), the MLE of the Poisson model is degenerate (or does not exist; Hashemi and Schneider, 2023). This corresponds to the cases where only a single lineage is detected in every sample (i.e., no evidence of super-infections), or if one lineage is observed in all samples. In such a situation, an MLE is still numerically found for the non-parametric model.

The performance of the non-parametric estimator is comparably good to that of the Poisson model even if the true MOI distribution is well approximated by a Poisson distribution. Despite the Poisson model being the true model for low true MOI, the non-parametric model yields less biased estimates of this parameter. This is also true for high average MOI, at least for a sufficiently large sample size. Moreover, if the true MOI distribution is highly over-dispersed, the non-parametric estimator is preferential. Irrespective of the differences, the results obtained by both models for the empirical examples from Cameroon, Kenya, and Venezuela are very similar. This and the fact that the differences between the estimators are small suggests that there is no need to drop the Poisson assumption, except for extreme parameters or in cases, in which there is clear evidence that the Poisson assumption is unjustified. Therefore, as a recommendation, it seems appropriate to use the non-parametric model in areas of low transmission, particularly, if transmission is highly heterogeneous. In such situations, the Poisson model is highly sensitive to outliers (Schneider, 2018). Areas of low transmission are becoming increasingly important with the ongoing malaria eradication efforts to decrease the global malaria burden by 90% and eradicate the disease in at least 35 countries by 2023 (World Health Organization, 2021). However, in the empirical example from Venezuela, at a time of low transmission, both methods yielded comparable results. In the example, this was due to the relatively large sample size (N = 97). In practice, molecular surveillance data might not always appropriately reflect the heterogeneity in transmission but is rather conducted in homogeneous populations. This renders the Poisson model to be sufficiently accurate.

The results here suggest that the non-parametric approach to estimating MOI is comparably good as the Poisson model. A particularly convenient property of the Poisson distribution is that it is characterized by a single parameter. This is no longer true for the non-parametric model, where the complete probability mass function of MOI needs to be estimated. Although there are many other parametric alternatives to the Poisson model, every other meaningful distribution is characterized by more than a single parameter. Moreover, from the structure outlined in Tsoungui Obama and Schneider (2022), such alternatives also do not retain the desirable properties of the Poisson model.

Given the data, i.e., the absence and presence of alleles in samples, the non-parametric model, is the most flexible assumption for MOI, since it does not require a certain distribution. In principle, the model should be able to approximate any parametric distribution, such as the conditional Poisson or conditional negative binomial distributions. Given that the non-parametric model does not clearly outperform the Poisson model implies that the latter already utilizes the information of data effectively, and any other model assumptions regarding over-dispersion (e.g., assuming a negative binomial distribution) will not outperform the Poisson model either. However, the statistical model can be extended to include further information in addition to the original data—for instance, one could extend the non-parametric model to account for patient-specific risks (e.g., one could group patients into different risk strata based on patient queries determining the risk of exposure, such as the use of bed nets, indoor residual spraying, window screens, etc., or recent malaria cases within the household). In the case of the Poisson model, it could be extended to estimate different MOI parameters for each stratum (risk group). Notably, such information has to be collected in addition to molecular data. In any case, our results suggest that, unless transmission is very heterogeneous, it is not necessary to extend the non-parametric approach to the case in which multiple genetic markers are studies at the same time as it was done for the Poisson model (e.g. Tsoungui Obama and Schneider, 2022; Obama and Schneider, 2023).

In summary, our results suggest that the non-parametric model to estimate MOI and lineage frequencies from single molecular markers (e.g., microsatellites, SNPs, and micro-haplotypes) introduced here is a valid approach. An implementation of the method as an R script is available in the Supplementary Materials and at GitHub https://github.com/Maths-against-Malaria/Non-parametric-MOI-estimation . Clearly, it is possible to extend the non-parametric model to more complex genetic architectures than a single molecular marker following the methodology in Tsoungui Obama and Schneider (2022).

The original contributions presented in the study are included in the article/ Supplementary Material . Further inquiries can be directed to the corresponding author.

LK: Conceptualization, Data curation, Formal analysis, Investigation, Methodology, Software, Visualization, Writing – original draft, Writing – review & editing. KS: Conceptualization, Formal analysis, Funding acquisition, Investigation, Methodology, Project administration, Resources, Software, Supervision, Validation, Visualization, Writing – original draft, Writing – review & editing.