The Casino Inside a Pandemic: How Unpredictable Timing Creates Predictable Catastrophe
The Invisible Hand Shaking the Dice
Somewhere in the mathematics of disease spread, there's a hidden casino. Every time a pathogen encounters a new host, the odds of transmission don't just depend on biology—they depend on timing. On whether an infected person coughs during the three days before they feel sick, or the two weeks after they've recovered. On whether recovery comes swiftly for some and drags on agonizingly for others. On the entire messy, unpredictable rhythm of infection and recovery across a population.
And in that casino, the house doesn't always win. Sometimes—rarely, but often enough to matter—the pathogen hits the jackpot. A single case becomes a hundred. A hundred become a hundred thousand. The tail of the outbreak-size distribution, the extreme events that overwhelm hospitals and collapse economies, has remained stubbornly difficult to predict.
Until now, perhaps.
A pair of physicists at the Hebrew University of Jerusalem have cracked a problem that has eluded epidemiologists for decades: how to predict the full distribution of epidemic outcomes—including the catastrophic outliers—when infection and recovery times don't follow the neat exponential patterns that make most models tractable. Their solution, published in July 2026, doesn't just improve predictions for one specific disease or one specific network. It reveals something unexpected: that beneath the chaos of real-world epidemics, a surprising simplicity lurks.
Ami Taitelbaum and Michael Assaf show that no matter how convoluted the actual timing of infection and recovery might be—no matter how far reality strays from the textbook assumptions—epidemic outcomes across wildly different scenarios collapse onto a single, predictable curve. The implications ripple from COVID forecasting to pandemic preparedness to our fundamental understanding of how diseases spread through human contact networks.
The Science
To understand why this matters, you need to understand what epidemiologists have been struggling with.
The workhorse of epidemic modeling is the SIR model—Susceptible, Infected, Recovered. Think of it as a three-state machine: people start susceptible, they can catch the disease from an infected neighbor, and once they've recovered, they're immune. The math is elegant, the logic is sound, and for decades it has underpinned our understanding of everything from measles to malaria.
But there's a catch. In the standard formulation, the model assumes Markovian dynamics—specifically, that the time an individual spends infected follows an exponential distribution. This is a mathematical convenience, not a biological fact. In reality, infections don't respect neat exponential rhythms. COVID's generation interval clustered around five to seven days, but with enormous variance. Some people infected with tuberculosis can carry the disease for years before transmitting it; others never transmit at all. Ebola's infectious period varied wildly depending on whether a patient survived and, if so, for how long.
These departures from exponentiality—technically called non-Markovian dynamics—create enormous analytical headaches. When the time between infection and recovery doesn't follow a predictable pattern, the system loses a crucial mathematical property that makes many predictions possible. The distribution of possible outcomes, including the probability of catastrophic outbreaks, becomes intractable.
"Non-Markovian epidemics arising from nonexponential waiting-time distributions render the dynamics more intricate and create substantial modeling and analytical challenges," Taitelbaum and Assaf note dryly in their paper. That's academic understatement for a problem that has stymied the field for years.
The challenge is compounded when you add network structure. In reality, people don't interact uniformly with everyone else. They have friends, colleagues, family members—specific connections through which disease can travel. A model that treats everyone as equally likely to infect everyone else (a "well-mixed" population) misses crucial features of real epidemics. Superspreaders exist. Some nodes in a contact network have dozens of connections; others have few. The topology of these networks dramatically shapes how outbreaks unfold.
Previous research has tackled these problems piecemeal. Some studies examined non-Markovian dynamics but ignored network structure. Others looked at network effects but assumed Markovian (exponential) waiting times. A few pushed the boundaries of both. But a unified theory predicting the full distribution of outbreak sizes—including the rare, extreme events that keep epidemiologists awake at night—remained elusive.
Taitelbaum and Assaf's approach is to find a bridge between these worlds. Their key insight is elegant: any arbitrary waiting-time distribution can be reduced to a single number called the per-edge transmissibility, denoted as T. Think of this as the effective probability that infection successfully passes along a single contact between a susceptible person and an infected one. Once you have this number, the full complexity of non-exponential infection and recovery times "maps" onto a simpler Markovian process—the kind that mathematicians know how to analyze.
The transmissibility T is computed by integrating the infection waiting-time distribution against the recovery survival function—in essence, asking: what's the probability that infection occurs before recovery? If infections happen quickly relative to recovery, T is high. If recovery is swift, T is low. This single number captures everything that matters about timing for the purpose of outbreak prediction.
To test their theory, Taitelbaum and Assaf ran extensive simulations. They generated networks of thousands of nodes—some regular (everyone has roughly the same number of connections), others following the Erdős–Rényi random graph model, and still others with heavy-tailed degree distributions mimicking real social networks. They assigned infection and recovery times following gamma distributions, which can be broader or narrower than exponential depending on their shape parameter α (alpha). And they tracked what happened when a single infected individual was introduced into each network.
The simulations were computationally intensive. For each combination of network structure, waiting-time parameters, and reproduction number (R₀—the average number of secondary infections each case generates), they ran thousands of epidemic realizations to map out the full distribution of outcomes. This gave them the data needed to test whether their theoretical predictions matched reality.
What They Found
The first major result is almost beautiful in its simplicity: for networks with weak degree heterogeneity—meaning most nodes have similar numbers of connections—outbreak statistics from radically different scenarios all collapse onto the same curve.
When Taitelbaum and Assaf rescaled outbreak-size distributions by centering them at their mode (the most likely outcome) and scaling by their conditional standard deviation, the data collapsed. Whether the infection times followed a gamma distribution with shape parameter α = 0.6 (much broader than exponential) or α = 1.2 (narrower), whether the network was regular or Erdős–Rényi, whether the mean degree was 6 or 16—the normalized distributions were indistinguishable.
Universal Collapse of Outbreak Statistics
| Label | Value |
|---|---|
| 0.10 | 0.1 |
| 0.15 | 0.15 |
| 0.20 | 0.2 |
| 0.25 | 0.25 |
| 0.30 | 0.3 |
| 0.35 | 0.35 |
| 0.40 | 0.4 |
| 0.45 | 0.45 |
This universal collapse is captured by the well-mixed WKB theory (named after the Wentzel–Kramers–Brillouin approximation from quantum mechanics). The theory predicts that the probability of observing an outbreak of fractional size x* follows P(x*) ~ e^(-N·𝒮(x*)), where N is the population size and 𝒮 is a large-deviation action that depends on the effective reproduction number. For weakly heterogeneous networks, the position on this universal curve is determined by what the researchers call R₀^eff—an effective reproduction number that encodes both the network structure and the waiting-time distributions.
The bond-percolation reproductive number plays a crucial role here. In bond percolation, each edge in the network is randomly designated as either "occupied" (transmission occurs) or "empty" (no transmission), with probability equal to the transmissibility T. The final outbreak is simply the connected component of the initial infected node in this subgraph. For weakly heterogeneous networks, this mapping works almost perfectly: the outbreak-size distribution predicted by bond percolation matches simulations, and the effective reproduction number R₀^eff can be computed from the deterministic final size.
The second key finding concerns what happens when you actually compute transmissibility for realistic waiting-time distributions. Taitelbaum and Assaf focused on gamma distributions, parameterized by shape α and rate r, with mean 1/r. When they plugged these into their transmissibility formula, they found something striking: the effect is highly asymmetric.
Broader-than-exponential infection times (α_inf < 1) dramatically increase transmissibility. Imagine a disease where some people become infectious almost immediately, while others take much longer. The early transmitters create chains of infection that bypass the slower individuals, effectively raising the overall transmission probability. Narrower-than-exponential infection times (α_inf > 1) have the opposite effect—they suppress transmissibility by concentrating infectiousness in a narrower window.
Recovery times matter too, but in a more limited way. Broader-than-exponential recovery (α_rec < 1) decreases transmissibility because it gives infection more time to occur before the infector recovers. But this effect is significant only when the distribution is extremely broad; for moderate deviations from exponentiality, recovery timing barely matters.
Mean Outbreak Size vs. Infection Shape Parameter
| Label | Value |
|---|---|
| α = 0.6 | 0.45 |
| α = 0.8 | 0.36 |
| α = 1.0 | 0.28 |
| α = 1.2 | 0.23 |
When the researchers computed the full outbreak-size distribution for Erdős–Rényi networks with various infection shape parameters, they found that their WKB theory with R₀^eff matched both the non-Markovian network simulations and well-mixed Markovian simulations with the same effective reproduction number. The mapping works. A non-Markovian epidemic on a network behaves, in terms of its outcome distribution, like a Markovian epidemic in a well-mixed population with an appropriately chosen reproduction number.
But this is only half the story. When degree heterogeneity increases—when the network has a broad distribution of connection counts, as real social networks do—the well-mixed mapping starts to break down.
For networks with high heterogeneity (coefficient of variation σ_k/k̄ > 0.5), the universal collapse persists for the conditional mean and standard deviation, but the full distribution shape starts to deviate from the well-mixed WKB prediction. The researchers found that on a gamma network with σ_k = 8 (highly heterogeneous), the well-mixed WKB theory fails badly: it overestimates the conditional standard deviation and mispredicts how the mean outbreak size changes with the infection shape parameter.
The solution is to retain the network structure but use an effective Markovian dynamics with a different effective reproduction number: R₀^net = T/(1-T) · λ_c^(-1), where λ_c is the Markovian epidemic threshold. When Taitelbaum and Assaf simulated Markovian epidemics on the same heterogeneous network with this R₀^net, they recovered the non-Markovian outbreak statistics almost perfectly.
This works even for empirical networks. The researchers tested their framework on the Hamsterster network—a real social network from a now-defunct online community, with 24,000 nodes, heavy-tailed degree distribution, and strong degree-degree correlations (assortative mixing). The network has properties that break many simplifying assumptions. Yet when they extended their bond-percolation mapping to account for degree correlations, the effective Markovian dynamics on the network captured the full non-Markovian outbreak-size distribution.
The practical implication is stark: even when the universal well-mixed curve fails, the effective Markovian approach succeeds. You don't need to abandon network structure to handle non-Markovian dynamics. You just need to recompute the effective reproduction number for the specific network and waiting-time distribution.
Why This Changes Things
Let's talk about what this actually means for the world.
Consider a pandemic preparedness scenario. Public health officials need to plan for the worst-case scenarios—not just the most likely outbreak size, but the probability that a pathogen causes catastrophic harm. Currently, this requires either detailed agent-based simulations (computationally expensive, requiring vast computing resources) or simpler models that assume exponential waiting times (potentially misleading for real diseases).
Taitelbaum and Assaf's framework offers a third path. Take measured waiting-time distributions from contact tracing data—Ferretti et al. (2020) have shown how to estimate these from genomic and epidemiological data—and compute the effective transmissibility T. Use this to calculate R₀^eff or R₀^net. Then apply the well-mixed WKB theory (for weakly heterogeneous populations) or the network-adapted Markovian mapping (for highly heterogeneous ones) to predict not just the mean outbreak, but the full distribution including extreme tails.
The numbers from the paper illustrate the stakes. On an Erdős–Rényi network of 3,000 people with mean degree 10 and R₀ = 1.5 (roughly the transmissibility of SARS-CoV-2 Delta variant), with infection times narrower than exponential (α_inf = 1.2), the mean outbreak infects about 28% of the population. The standard deviation is around 7%. But here's what the distribution actually looks like: approximately 25% of outbreaks fall more than 20% below the mean, 32% fall within 20% of the mean, and 21% fall more than 20% above the mean.
That 21% figure is the terrifying part. Even in a scenario with moderate average transmissibility and only modest deviations from exponentiality, there's a one-in-five chance of an outbreak that devastates 35% or more of the population. Given the right conditions—broader infection times, higher heterogeneity—the tail risk grows even larger.
The Hamsterster network results tell a similar story. With R₀ = 3 and α_inf = 1.4, the mean outbreak is smaller (~13% of the population) but the right tail still carries about 19% of the probability mass. "In both examples, the right tail carries substantial probability mass of about 20%, indicating a significant risk of unusually large outbreaks," the researchers write.
This matters because standard epidemic planning often focuses on the mean or the threshold condition (R₀ > 1 vs. R₀ < 1). The existence of heavy tails means that focusing on these summary statistics systematically underestimates the risk of extreme events. It's the epidemiological equivalent of planning for average weather while ignoring hurricanes.
There are broader theoretical implications too. The discovery that outbreak statistics collapse onto a universal curve—for weakly heterogeneous networks—is a surprising simplification. It suggests that the details of waiting-time distributions and network structure matter far less than the effective reproduction number. This echoes similar universalities discovered in other areas of statistical physics, where complex microscopic dynamics give rise to simple macroscopic behavior.
The framework also unifies previously separate threads of research. Studies of non-Markovian SIS dynamics (where individuals can be reinfected) treated transmissibility differently from SIR dynamics (where recovery confers permanent immunity). Taitelbaum and Assaf show that for SIR, the shape-dependence of transmissibility is essential—it's not just the effective rate that matters, but how the waiting-time distribution is shaped. Their framework captures both infection and recovery non-Markovianity simultaneously, something previous approaches achieved only partially.
From a practical standpoint, the work opens new avenues for risk assessment. Real epidemic data—contact tracing studies, genomic surveillance, wastewater monitoring—can be used to estimate waiting-time distributions. These can be fed into the transmissibility calculation, yielding effective reproduction numbers that account for the full complexity of timing. The resulting predictions should be more accurate than methods that assume exponentiality, without requiring the computational overhead of full agent-based simulation.
What's Next
Several questions remain open. The WKB theory assumes large population sizes; for small populations (say, a single household or a small community), finite-size effects become important and the universal collapse may break down. The left tail of the distribution—corresponding to very small outbreaks or rapid fade-outs—is less universal than the right tail, and understanding these finite-size corrections is an important direction for future work.
The framework also assumes tree-like network structure (no short cycles), which holds well for sparse random graphs but may fail for dense networks or those with strong clustering. Real social networks often exhibit substantial clustering—your friends are also friends of each other—which could affect the distribution shape. Extending the theory to clustered networks is a natural next step.
Another avenue is temporal networks, where connections change over time. Real contact patterns aren't static; they vary by time of day, day of week, and season. The current framework treats the network as fixed; adapting it to time-varying contact patterns would increase its practical utility for real-time epidemic forecasting.
There's also the question of multi-strain dynamics. If a pathogen evolves, or if there are multiple variants circulating simultaneously, the mapping to a single transmissibility value may no longer suffice. The interplay between non-Markovian dynamics and strain competition remains unexplored.
Most ambitiously, the framework could be extended to other compartmental models beyond SIR. The SEIR model (adding an Exposed compartment between Susceptible and Infected) is widely used for diseases like COVID-19. Whether the universal collapse persists for SEIR dynamics, and whether the transmissibility mapping can be generalized, are questions for future research.
The practical path forward is clear: calibrate the framework to real pathogen data. Contact tracing studies have already produced estimates of generation intervals and infectious periods for numerous diseases. These can be used to compute effective transmissibility and reproduction numbers, which can then be fed into the outbreak-size prediction framework. The result would be actionable risk bounds—probability distributions for outbreak sizes under various intervention scenarios—that go beyond the point estimates currently used in pandemic planning.
Taitelbaum and Assaf are careful to note that their framework doesn't predict what any specific outbreak will do—it predicts the distribution of possible outcomes. But that's exactly what risk assessment needs. We're not asking "will this particular outbreak be catastrophic?" We're asking "what's the probability of a catastrophic outbreak, and how does that probability change if we intervene?" The framework gives us tools to answer those questions with unprecedented mathematical rigor.
The invisible casino of epidemic dynamics has rules. Taitelbaum and Assaf have started to decode them.
The original paper, "Extreme outbreaks in non-Markovian epidemics on complex networks" by Ami Taitelbaum and Michael Assaf (Hebrew University of Jerusalem), was published on arXiv on July 27, 2026.