Main

For sexually transmitted and vector-borne diseases, host contact rates have long served as surrogates for individual infectiousness3,12,13,14, leading to the assertion of a general ‘20/80 rule’ (whereby 20% of cases cause 80% of transmission13) and to the influential concept of high-risk ‘core groups’3,12,13. For directly transmitted infections, however, the overall infectiousness of each case—that is, the number of other individuals infected during the infectious lifetime of a single individual—arises from a complex mixture of host, pathogen and environmental factors (see Supplementary Notes). Consequently, the degree of infectiousness is distributed continuously in any population4,7,11,15,16 and, crucially, distinct risk groups often cannot be defined a priori2,11. This impedes the conventional approach to adding heterogeneity to epidemic models, in which populations are divided into homogeneous subgroups2,3,4,17. Research on continuous individual variation in infectiousness for directly transmitted infections has been largely restricted to within-household transmission18,19, or to variation in infectious period20,21 or social network22. Some recent studies have used contact tracing data to investigate specific questions in light of observed variation8,16, but a broad understanding of the role of individual variation in outbreak dynamics is lacking.

As a theoretical basis, we introduce the ‘individual reproductive number’, ν, as a random variable representing the expected number of secondary cases caused by a particular infected individual. Values for ν are drawn from a continuous probability distribution with population mean R0 that encodes all variation in infectious histories of individuals, including properties of the host and pathogen and environmental circumstances. In this framework, superspreading events (SSEs) are not exceptional events9, but important realizations from the right-hand tail of a distribution of ν (refs 7, 15). Stochastic effects in transmission are modelled using a Poisson process4, so that the number of secondary infections caused by each case, Z, is described by an ‘offspring distribution’ Pr(Z = k) where Z∼Poisson(ν).

By considering three possible distributions of ν, we generate three candidate models for the offspring distribution: (1) in generation-based models neglecting individual variation, ν = R0 for all cases, yielding Z∼Poisson(R0); (2) in differential-equation models with homogeneous transmission and constant recovery rates, ν is exponentially distributed, yielding Z∼geometric(R0); (3) in a more general formulation, we let ν be gamma-distributed with mean R0 and dispersion parameter k, yielding Z∼negative binomial(R0,k) (ref. 23). The negative binomial model includes the conventional Poisson (k → ∞) and geometric (k = 1) models as special cases. It has variance R0(1 + R0/k), so smaller values of k indicate greater heterogeneity.

We gathered empirical offspring distributions from detailed contact tracing or surveillance data sets, and challenged the candidate models using model selection techniques24 (see Supplementary Notes). For SARS outbreaks in Singapore and Beijing, the negative binomial model is unequivocally favoured (Fig. 1a and Supplementary Table 1). Conventional models assuming homogeneity cannot reproduce the observed transmission patterns. For the Singapore outbreak, the maximum-likelihood estimate k̂ is 0.16 (90% confidence interval 0.11–0.64), indicating an underlying distribution of ν that is highly overdispersed (Fig. 1a, inset). According to this analysis, the great majority of SARS cases in Singapore were barely infectious (73% had ν < 1) but a small proportion were highly infectious (6% had ν > 8). This is consistent with field reports from SARS-afflicted regions5,6 but contrasts with published SARS models9,10,25,26.

Figure 1: Evidence for variation in individual reproductive number ν.
figure 1

a, Transmission data from the SARS outbreak in Singapore in 2003 (ref. 5). Bars show observed frequency of Z, the number of individuals infected by each case. Lines show maximum-likelihood fits for Z∼Poisson (squares), Z∼geometric (triangles), and Z∼negative binomial (circles). Inset, probability density function (solid) and cumulative distribution function (dashed) for gamma-distributed ν (corresponding to Z∼negative binomial) estimated from Singapore SARS data. b, Expected proportion of all transmission due to a given proportion of infectious cases, where cases are ranked by infectiousness. For a homogeneous population (all ν = R0), this relation is linear. For five directly transmitted infections (based on k̂ values in Supplementary Table 1), the line is concave owing to variation in ν. c, Proportion of transmission expected from the most infectious 20% of cases, for 10 outbreak or surveillance data sets (triangles). Dashed lines show proportions expected under the 20/80 rule (top) and homogeneity (bottom). Superscript ‘v’ indicates a partially vaccinated population. d, Reported superspreading events (SSEs; diamonds) relative to estimated reproductive number R (squares) for twelve directly transmitted infections. Lines show 5–95 percentile range of Z∼Poisson(R), and crosses show the 99th-percentile proposed as threshold for SSEs. Stars represent SSEs caused by more than one source case. ‘Other’ diseases are: 1, Streptococcus group A; 2, Lassa fever; 3, Mycoplasma pneumonia; 4, pneumonic plague; 5, tuberculosis. R is not shown for ‘other’ diseases, and is off-scale for monkeypox. See Supplementary Notes for details.

Full size image

Comparing results for eight directly transmitted infections reveals the differing degree of individual variation among diseases and outbreak settings (Fig. 1b, c and Supplementary Tables 1, 2). The Poisson offspring distribution is almost always strongly rejected. The geometric model has considerable support for several data sets, which indicates significant individual variation in transmission rates because real infectious periods are less dispersed than the exponential distribution20,21. The negative binomial model is selected decisively for several data sets, and enables comparative study of diseases through the dispersion parameter. Like SARS, measles in highly vaccinated populations shows high variation in two surveillance data sets, with narrow confidence intervals excluding the conventional models (note that heterogeneous vaccination coverage is an important environmental factor contributing to this pattern). Monkeypox and smallpox viruses show intermediate variation, consistent across multiple data sets, and pneumonic plague transmission is slightly less variable. Data limitations prevent definitive conclusions for other diseases. Comparing our findings to the 20/80 rule proposed for sexually transmitted and vector-borne diseases13, no general rule emerges but the core principle of heterogeneous transmission is certainly supported (Fig. 1c).

Numerous reports of superspreading events provide further evidence for variation in ν. We reviewed 37 published accounts of SSEs for 11 directly transmitted infections (Fig. 1d; see Supplementary Notes). Unrecognized or misdiagnosed illness is the most common cause of these SSEs, followed by alternative modes of spread (especially airborne), high contact rates, and co-infections that aid transmission. High pathogen load or shedding rates are occasionally implicated, but are rarely measured. A consistent and general definition of SSEs is currently lacking—for SARS, an SSE has been arbitrarily defined as Z ≥ 8 (ref. 6), Z ≥ 10 (ref. 5), Z > 10 (ref. 26) or ‘many more than the average number’9, and different thresholds are surely needed for measles (R0∼11–18; ref. 3) or monkeypox (R0 < 1).

We propose this general protocol for defining a superspreading event: (1) estimate the effective reproductive number, R, for the disease and population in question; (2) construct a Poisson distribution with mean R, representing the expected range of Z due to stochasticity without individual variation; (3) define an SSE as any infected individual who infects more than Z(n) others, where Z(n) is the nth percentile of the Poisson(R) distribution. A 99th-percentile SSE is then any case causing more infections than would occur in 99% of infectious histories in a homogeneous population (Fig. 1d). This approach complements a priori identification of potential superspreaders when that is feasible, as for sexually transmitted diseases (where promiscuity drives risk)3,12. In addition, the definition enables prediction of the frequency of SSEs once R0 and k have been estimated (Supplementary Fig. 1)—an outstanding challenge in emerging disease epidemiology8,9.

To assess the effect of individual variation on disease outbreaks, we analyse a branching process model with negative binomial offspring distribution, corresponding to gamma-distributed ν (Fig. 2a; see Supplementary Notes). Of primary interest is the probability of stochastic extinction, q, after the introduction of a single infected individual (Fig. 2b). For R0 < 1, all invasions die out, as in standard models. For R0 > 1, increased variation strongly favours extinction8. For example, if R0 = 3 then q = 0.06 under the assumption of homogeneous ν (k → ∞), or q = 0.33 if k = 1, but if k = 0.16 (as estimated for SARS) then q = 0.76. Extinction risk rises owing to a higher proportion of non-transmitting cases when ν is overdispersed (Figs 1a, 2a and Supplementary Fig. 2a). This effect thwarts invasion by diseases that are very potent spreaders on average: for arbitrarily high R0, q → 1 as k → 0 (Supplementary Fig. 2b). The expected number of cases before extinction is hardly affected by k (Supplementary Fig. 2c), because low-k outbreaks that fail probably lacked SSEs and thus resemble homogeneous outbreaks with lower R0. Accordingly, when individual variation is large, extinction occurs rapidly or not at all (Supplementary Fig. 2d).

Figure 2: Outbreak dynamics with different degrees of individual variation in infectiousness.
figure 2

a, The individual reproductive number ν is drawn from a gamma distribution with mean R0 and dispersion parameter k. Probability density functions are shown for six gamma distributions with R0 = 1.5 (‘k = Inf’ indicates k → ∞). b, Probability of stochastic extinction of an outbreak, q, versus population-average reproductive number, R0, following introduction of a single infected individual. The value of k increases from top to bottom (values and colours as in a). c, Growth of simulated outbreaks with R0 = 1.5 and one initial case, conditional on non-extinction. Boxes show median and interquartile range (IQR) of the first disease generation with 100 cases; whiskers show most extreme values within 1.5 × IQR of the boxes, and crosses show outliers. Percentages show the proportion of 10,000 simulated outbreaks that reached the 100-case threshold (roughly 1 - q).

Full size image

For outbreaks avoiding stochastic extinction, epidemic growth rates strongly depend on variation in ν (Fig. 2c and Supplementary Fig. 2e, f). Diseases with high individual variation show infrequent but explosive epidemics after introduction of a single case. This pattern recalls SARS in 2003, for which many settings experienced no epidemic despite unprotected exposure to SARS cases27,28, whereas a few cities suffered explosive outbreaks8,9,10,15,26. Our results, using k̂ = 0.16 for SARS, explain this simply by the presence or absence of high-ν individuals in the early generations of each outbreak6. In contrast, conventional models (with k = 1 or k → ∞) cannot simultaneously generate frequent failed invasions and rapid growth rates without additional, subjective model structure.

Disease control interventions could increase or decrease individual variation in infectiousness. Infected individuals might reduce their number of non-essential contacts, or governments might impose quarantine or isolation on particular individuals. Here we explore several idealized cases theoretically, for an outbreak with offspring distribution Z∼negative binomial(R0,k) before control (see Supplementary Notes). Consider the effect of control effort c, where c = 0 reflects no control and c = 1 reflects complete blockage of transmission. Under population-wide control, the infectiousness of every individual in the population is reduced by a factor c (that is, for all individuals). Under random, individual-specific control, a proportion c of infected individuals (chosen at random) is traced and isolated completely such that they cause zero infections (that is for a proportion c of infected individuals, and for the rest). Individual-specific control raises the degree of heterogeneity in the outbreak as measured by the variance-to-mean ratio of Z, whereas population-wide control reduces heterogeneity. Both approaches yield effective reproductive number R = (1 - c)R0, so the threshold control effort for guaranteed disease extinction is c ≥ 1 - 1/R0 as in conventional models. For intermediate values of c, however, the individual-specific approach always works better (Fig. 3a and Supplementary Fig. 3a, b), consistent with our finding that higher variation favours disease extinction (Fig. 2b). Branching process theory confirms that qind > qpop whenever c∈(0,1 - 1/R0) (see Supplementary Notes).

Figure 3: Implications for control measures.
figure 3

a, Increase in extinction probability (qind - qpop) under individual-specific control compared to population-wide control, for diseases with R0 = 3 and different degrees of individual variation, k, subject to control effort c. With population-wide control, the infectiousness of all individuals is reduced by a factor c. With individual-specific control, a proportion c of infectious individuals (selected at random) have their infectiousness reduced to zero. The outbreak is assumed to begin with one case, with control present from the outset. b, Estimates of R̂ and k̂ from outbreak data sets before and after control measures were initiated (joined by solid lines; Supplementary Table 2), and post-control values of kc estimated from theoretical models of control as described in the Supplementary Notes. c, Effect of random versus targeted control measures. The probability of outbreak containment (defined as never reaching the 100-case threshold) for four diseases with R0 = 3 and k = 0.1 (blue), k = 0.5 (green), k = 1 (black) or k → ∞ (purple). Control policies are population-wide (solid lines), random individual-specific (dotted lines), or targeted individual-specific (dashed lines, where half of all control effort is focused on the most infectious 20% of cases). For k → ∞, all individuals are identical, so targeting has no effect and dotted and dashed lines overlay one another. d, The factor by which targeting increases the effect of control on preventing a major outbreak, relative to random individual-specific control (see Supplementary Notes), when 20%, 40% or 60% of the total population is controlled. Results in c and d are the mean of 10,000 simulations, with control beginning in the second generation of cases.

Full size image

To assess the realism of these idealized control scenarios, we analysed contact tracing data from four outbreaks before and after imposition of control measures. Control always lowered the estimated dispersion parameter (that is k̂c < k̂) as predicted by the individual-specific model (Fig. 3b), although small sample sizes often led to overlapping confidence intervals (Supplementary Table 2). This increased skew in transmission arose chiefly from undiagnosed or misdiagnosed individuals, who continued to infect others (and even cause SSEs), whereas controlled individuals infected very few. To further examine our control theories, we calculated and for each data set; k̂c was always closer to , although twice it fell between the two predictions, indicating a possible combination of control mechanisms (Fig. 3b). Real-world control thus seems to increase individual variation, favouring extinction but risking ongoing SSEs. Larger data sets are needed to establish this pattern definitively.

If highly infectious individuals can be identified predictively (see Supplementary Notes) then the efficiency of control could be greatly increased (Fig. 3c, d). Focusing half of all control effort on the most infectious 20% of cases is up to threefold more effective than random control (Fig. 3d). When k = 0.1 or 0.5, outbreak containment is assured for targeted control levels at roughly half the threshold level of c = 1 - 1/R0 for random control. Gains in efficiency increase with more intense targeting of high-ν cases, but saturate as overall coverage c increases (Supplementary Fig. 3c, d). Again, branching process theory generalizes these findings: for a given proportion c of individuals controlled, greater targeting of higher-ν individuals leads to lower effective reproductive number R and higher extinction probability q (see Supplementary Notes).

The data sets analysed here were collected from published literature, and may be subject to selection bias for successful invasions and SSEs rather than typical disease behaviour. Surveillance data sets are less vulnerable to this bias, but may under-report isolated cases. We urge that detailed transmission tracing data be collected and made public whenever possible, even if unexceptional. At a minimum, we propose a new measure for inclusion in outbreak reports: the proportion of cases not transmitting (p0), which, together with R0 is sufficient to estimate the degree of variation in ν (Supplementary Fig. 4). As more data become available, trends may emerge in the degree of variation present, for example, for different modes of spread or levels of virulence. Richer data sets may also enable testing of the branching process assumption that case outcomes are independent and identically distributed, by detecting possible correlations in ν values within transmission lineages or systematic changes as outbreaks progress.

Our results have broad implications for emerging disease epidemiology, and open challenges for further work. Explosive epidemics demand rapid action by authorities and can strain health infrastructures. High extinction probabilities indicate that disease introductions or host species jumps may be more frequent than currently suspected. Cluster-size surveillance for pathogen adaptation29 or dwindling population immunity30 should be tuned to observed levels of variation. Realization of targeted control measures requires a better understanding of factors determining individual infectiousness. This work must be integrated with established theory of sexually transmitted diseases and social networks, where high-risk groups exert nonlinear influence on R0 because contact rates affect infectiousness and susceptibility equally3,4,12,13,22. All diseases probably show intermediate degrees of covariation between infectiousness and susceptibility, a topic demanding empirical and theoretical study17. The central role of R0 in epidemic analysis is unassailable, but our findings show that emerging disease outbreaks cannot be fully understood if individual variation in infectiousness is neglected. Examination of other population processes dependent on small numbers of individuals may yield similar insights.

Methods

Analysis of disease data

For data sets including the full distribution of Z, we estimated R̂0 and k̂ using maximum-likelihood methods. The candidate models were compared using Akaike's information criterion (AICc) modified for small sample size. Confidence intervals for k̂ were estimated by bias-corrected non-parametric bootstrapping and corroborated by four other methods. For data sets including only estimates of R̂0 and the proportion of cases not transmitting (p̂0), we estimated k̂ by solving p̂0 = (1 + R̂0/k)-k numerically, and evaluated the candidate models using confidence intervals calculated by two methods. Expected proportions of transmission due to particular groups of infectious individuals (Fig. 1b, c) were calculated using the gamma distribution of ν with estimated values of R̂0 and k̂. See Supplementary Notes for details, and for descriptions of data sets.

Branching process analysis

Analysis of branching process models centres on the probability generating function (pgf) of the offspring distribution, , defined for |s| ≤ 1. When R0 > 1, the long-term probability of disease extinction after introduction of a single infected individual is the unique solution of q = g(q) on the interval (0,1). For a negative binomial offspring distribution Z ≈ NegB(R0,k), the pgf is . Under population-wide control, and therefore , and the variance-to-mean ratio is 1 + (1 - c)R0/k. Under random individual-specific control, the exact pgf is with variance-to-mean ratio 1 + R0/k + cR0. This scenario can be approximated by , where is the solution to and decreases monotonically as c increases. Further details, descriptions of outbreak simulations and formal analysis of control measures are found in the Supplementary Notes.