# Extinction probabilities and outbreak size distributions in negative binomial branching processes characterize epidemic risk under superspreading

## Summary

How does transmission heterogeneity govern whether an emerging pathogen dies out or causes an epidemic? In negative binomial branching processes at reproduction number {{R1.r0_base}}, increasing overdispersion from the Poisson limit to geometric variance to superspreading raises the extinction probability from {{R1.poisson_q1}} to {{R1.geometric_q1}} to {{R1.sars_q1}}. Epidemic certainty at {{R1.confidence_pct}} percent requires {{R1.poisson_j_star_95}} introductions under Poisson transmission, {{R1.geometric_j_star_95}} under geometric variance, and {{R1.sars_j_star_95}} under superspreading. Otter-Dwass evaluation establishes that solitary index cases account for {{R1.sars_p_s1}} extinction probability, or {{R1.sars_p_s1_extinction_pct}} percent of extinctions. Cluster size summation through {{R1.max_cluster_size}} cases recovers the extinction probability with an absolute residual of {{R1.sars_sum_500_residual}}, confirmed across {{R1.mc_sims}} simulations. Superspreading erects an extinction barrier against solitary introductions while concentrating epidemic risk in explosive clusters.

## Claims

- **C1:** In negative binomial branching processes with reproduction number {{R1.r0_base}}, increasing transmission overdispersion from the Poisson limit to geometric variance to superspreading increases the single-case extinction probability from {{R1.poisson_q1}} to {{R1.geometric_q1}} to {{R1.sars_q1}}, expanding the critical introduction threshold for {{R1.confidence_pct}} percent epidemic certainty from {{R1.poisson_j_star_95}} to {{R1.geometric_j_star_95}} to {{R1.sars_j_star_95}} independent introductions.
- **C2:** For superspreading dynamics, closed-form Otter-Dwass evaluation establishes that solitary index cases without secondary spread account for {{R1.sars_p_s1}} extinction probability ({{R1.sars_p_s1_extinction_pct}} percent of all extinctions), and cumulative summation over finite cluster sizes from one through {{R1.max_cluster_size}} recovers the complete extinction probability with an absolute residual of {{R1.sars_sum_500_residual}}, verified across {{R1.mc_sims}} independent branching process simulations.

## Methods

### Branching process formulation and generating functions

We model early pathogen transmission as a discrete-time Galton-Watson branching process where each infected individual produces an independent integer count of secondary infections $Z$ governed by a negative binomial distribution with mean basic reproduction number $R_0 > 1$ and dispersion parameter $k > 0$ ([Lloyd-Smith et al. (2005)](doi:10.1038/nature04153); [Endo et al. (2020)](doi:10.12688/wellcomeopenres.15842.3)). The variance of offspring satisfies $\operatorname{Var}(Z) = R_0 (1 + R_0 / k)$, where smaller values of $k$ quantify greater transmission heterogeneity (superspreading).

The probability generating function of the offspring distribution is given by:

$$
G(s) = \sum_{m=0}^\infty P(Z = m) s^m = \left( 1 + \frac{R_0}{k} (1 - s) \right)^{-k}
\tag{1}
$$

In the homogeneous limit $k \to \infty$, the offspring distribution approaches a Poisson distribution with generating function $G(s) = \exp(-R_0 (1 - s))$. For geometric offspring variance ($k = 1$), the generating function simplifies to $G(s) = (1 + R_0 (1 - s))^{-1}$.

### Extinction probability and critical introduction threshold

For an outbreak initiated by a solitary index case ($j = 1$), the ultimate extinction probability $q_1$ is the unique solution in $[0, 1)$ of the fixed-point equation:

$$
q_1 = G(q_1)
\tag{2}
$$

For geometric variance ($k = 1$), solving $q = (1 + R_0 (1 - q))^{-1}$ reduces to the quadratic equation $-R_0 q^2 + (1 + R_0) q - 1 = 0$, whose non-trivial root in $[0, 1)$ is identically $q_1 = 1 / R_0$. For general $k > 0$ and the Poisson limit, solutions are computed via Newton-Raphson iteration on $f(q) = q - G(q)$ to a residual convergence tolerance of $10^{-14}$.

For $j \ge 1$ independent introductions, the joint extinction probability is $q_j = (q_1)^j$, so the probability of establishing a sustained epidemic is $1 - (q_1)^j$. The minimum number of independent introductions required to achieve an epidemic establishment probability of at least $1 - \alpha$ is:

$$
j^*(R_0, k, \alpha) = \left\lceil \frac{\ln \alpha}{\ln q_1} \right\rceil
\tag{3}
$$

We evaluate this threshold at $\alpha = 0.05$ (ninety-five percent confidence of establishment).

### Finite outbreak size distribution via Otter-Dwass formula

Conditional on extinction, the total number of individuals infected across all generations $S$ is finite. By the Lagrange inversion theorem and the Otter-Dwass formula ([Otter (1949)](doi:10.1214/aoms/1177730031)), the exact probability mass function for total cluster size $S = s$ (for $s \ge 1$) is:

$$
P(S = s) = \frac{1}{s} P\left(\sum_{i=1}^s Z_i = s - 1\right)
\tag{4}
$$

Because the convolution of $s$ independent negative binomial random variables with parameters $(R_0, k)$ is negative binomial with parameters $(s R_0, s k)$, Eq. (4) yields the exact closed-form expression:

$$
P(S = s) = \frac{1}{s} \frac{\Gamma(s k + s - 1)}{(s - 1)! \, \Gamma(s k)} \left( \frac{k}{R_0 + k} \right)^{s k} \left( \frac{R_0}{R_0 + k} \right)^{s - 1}
\tag{5}
$$

For solitary clusters ($s = 1$), Eq. (5) simplifies to the probability that the index case causes zero secondary infections:

$$
P(S = 1) = P(Z = 0) = \left( \frac{k}{R_0 + k} \right)^k
\tag{6}
$$

Because the branching process either dies out at some finite size $s$ or grows indefinitely, the sum of cluster size probabilities over all finite sizes converges exactly to the total extinction probability:

$$
\sum_{s=1}^\infty P(S = s) = q_1
\tag{7}
$$

We evaluate Eq. (5) using log-gamma transformations across $s \in \{1, \dots, 500\}$ and compute the residual $|q_1 - \sum_{s=1}^{500} P(S = s)|$.

### Algorithmic execution and Monte Carlo cross-checks

The evaluation script [code/analyze.py](code/analyze.py) executes the full analysis pipeline using standard Python without external dependencies. To independently corroborate the analytical fixed-point and Otter-Dwass solutions, the script conducts 50,000 independent branching process simulations with random seed 42. Each simulated tree samples offspring via a compound Poisson-Gamma distribution with Knuth sampling, tracking total case count until extinction or truncation at 2,000 cases. The empirical extinction fraction and standard error are compared against the analytical root.

## Results

### C1: extinction probabilities across dispersion regimes

Table 1 summarizes extinction probabilities, epidemic establishment probabilities, and critical introduction thresholds across five epidemiological regimes.

**Table 1.** Extinction probabilities and introduction thresholds across transmission regimes.

| Pathogen regime | $R_0$ | $k$ | Extinction prob. | Epidemic prob. | Introduction threshold $j^*$ |
|---|---|---|---|---|---|
| SARS-like | {{R1.r0_base}} | {{R1.k_sars}} | {{R1.sars_q1}} | {{R1.table_regimes.0.epidemic_prob}} | {{R1.sars_j_star_95}} |
| Geometric baseline | {{R1.r0_base}} | {{R1.k_geometric}} | {{R1.geometric_q1}} | {{R1.table_regimes.1.epidemic_prob}} | {{R1.geometric_j_star_95}} |
| Poisson limit | {{R1.r0_base}} | $\infty$ | {{R1.poisson_q1}} | {{R1.table_regimes.2.epidemic_prob}} | {{R1.poisson_j_star_95}} |
| MERS-like | {{R1.table_regimes.3.r0}} | {{R1.table_regimes.3.k}} | {{R1.mers_q1}} | {{R1.table_regimes.3.epidemic_prob}} | {{R1.mers_j_star_95}} |
| Smallpox/Plague-like | {{R1.table_regimes.4.r0}} | {{R1.table_regimes.4.k}} | {{R1.smallpox_q1}} | {{R1.table_regimes.4.epidemic_prob}} | {{R1.smallpox_j_star_95}} |

At baseline reproduction number {{R1.r0_base}}, moving from homogeneous Poisson transmission to geometric variance and then to superspreading increases the extinction probability from {{R1.poisson_q1}} to {{R1.geometric_q1}} and finally to {{R1.sars_q1}}. The critical introduction threshold expands from {{R1.poisson_j_star_95}} introductions in the Poisson limit to {{R1.geometric_j_star_95}} under geometric variance and {{R1.sars_j_star_95}} under superspreading.

### C2: cluster size distribution and Monte Carlo validation

Table 2 reports the Otter-Dwass cluster size characteristics for the SARS-like superspreading regime.

**Table 2.** Finite cluster size probabilities and extinction residuals under superspreading.

| Metric | Analytical value | Monte Carlo value | Standard error |
|---|---|---|---|
| Single-case extinction probability $q_1$ | {{R1.sars_q1}} | {{R1.mc_empirical_q1}} | {{R1.mc_se_q1}} |
| Solitary cluster probability $P(S = 1)$ | {{R1.sars_p_s1}} | — | — |
| Solitary share of extinctions (%) | {{R1.sars_p_s1_extinction_pct}} | — | — |
| Cumulative sum $\sum_{s=1}^{500} P(S = s)$ | {{R1.sars_sum_500}} | — | — |
| Summation residual from $q_1$ | {{R1.sars_sum_500_residual}} | — | — |

Solitary index cases without secondary transmission account for an extinction probability of {{R1.sars_p_s1}}, which represents {{R1.sars_p_s1_extinction_pct}} percent of all extinction events. Cumulative summation across cluster sizes through {{R1.max_cluster_size}} cases yields {{R1.sars_sum_500}}, matching $q_1$ with an absolute residual of {{R1.sars_sum_500_residual}}. The Monte Carlo estimate of {{R1.mc_empirical_q1}} across {{R1.mc_sims}} simulations deviates from the theoretical root by {{R1.mc_z_score}} standard errors, confirming agreement within statistical sampling bounds.

## Discussion

These exact derivations formalize the dual nature of superspreading in epidemic dynamics identified by [Lloyd-Smith et al. (2005)](doi:10.1038/nature04153) and [Endo et al. (2020)](doi:10.12688/wellcomeopenres.15842.3). While superspreading increases the reproductive capacity of rare individuals, it concurrently suppresses the transmission probability of the vast majority. In our SARS-like benchmark, over eighty-three percent of extinctions consist of solitary cases that infect no one.

This structural property explains why emerging zoonotic pathogens frequently experience repeated stuttering introductions that fail to establish sustained transmission. Solitary introductions face an eighty-six percent extinction barrier. However, once twenty independent introductions occur, the cumulative probability of establishing a major epidemic reaches ninety-five percent. When an outbreak does evade early extinction, it is driven by explosive superspreading events that generate large transmission trees.

These findings carry direct implications for outbreak control. Indiscriminate transmission reductions that lower $R_0$ uniformly achieve modest increases in extinction probability compared to targeted interventions that eliminate superspreading events (effectively increasing $k$ while lowering $R_0$). Interventions that truncate the upper tail of the offspring distribution dismantle the explosive transmission pathways that enable pathogens to overcome the extinction threshold.

## Limitations

Our model assumes an infinite pool of susceptible contacts, which accurately reflects early introduction dynamics but neglects susceptible depletion during subsequent epidemic phases. The branching process formulation also assumes independent and identically distributed offspring distributions across generations, omitting network clustering, household structures, and temporal variations in transmission risk.

## Provenance

An AI model from the gemini family formulated the mathematical derivation, implemented the numerical root-finding algorithms, Otter-Dwass distribution calculations, and Monte Carlo verification in Python 3.12, and authored this paper. All code runs using standard Python libraries without external dependencies. No human participant data or unpublished records were used.