Why Doubling the Epidemic Peak Sometimes Works — and How to Fix It When It Doesn't

In the early weeks of a pandemic, the question epidemiologists grapple with isn't just how many people will eventually catch a virus. It's when — and how severely — the wave will crest. Hospital beds fill, ICUs strain, ventilators run short: all of it peaks around a single moment when the number of actively infectious people reaches its maximum. Getting that moment right matters enormously. Get it wrong and you either underprepare and overwhelm your health system, or overprepare and waste resources that could save lives elsewhere.

So here's a quiet, foundational problem: how do you estimate the peak prevalence of an epidemic when you can't actually measure prevalence directly?

Case reports exist. Testing counts exist. But true point prevalence — the fraction of a population actively infectious at any moment — gets obscured by reporting lags, asymptomatic spread, uneven testing, and the fact that most health systems can only count what they see. Epidemiologists have long worked around this by using proxies: quantities they can observe, from which they try to recover or bound the thing they actually want to know.

A new paper by Denis Tverskoy, Andrew Gothard, and Grzegorz Rempala at The Ohio State University tackles this problem head-on. Their work, "Approximating Peak Prevalence in Multistage SIR Epidemics," doesn't just offer a better formula. It exposes something deeper about why the formulas we've been using work when they do — and why they fail precisely when we might need them most. Their central finding: the widely-used "factor-two approximation" — a rule of thumb that prevalence peak is roughly twice some observable weighted quantity — is mathematically justified under one specific scaling of infectious stages, but breaks down badly under another. More importantly, they develop corrections that dramatically improve accuracy, with error rates falling from 20-40% down to under 5% in many practical regimes.

That's not just an academic exercise. It's a blueprint for better real-time epidemic modeling.

The Science

To understand what Tverskoy and colleagues accomplished, you need to understand the framework epidemiologists use to model how diseases spread through populations: compartmental models.

The classic formulation, developed by Kermack and McKendrick in 1927, divides a population into three buckets: Susceptible people (S) who can catch the disease, Infectious people (I) who currently have it, and Removed people (R) who have either recovered with immunity or died. An SIR model tracks how individuals flow between these compartments over time, governed by parameters like the transmission rate (β) and the recovery rate (γ).

For the basic SIR model — where infectious periods follow an exponential distribution, a mathematical convenience that makes the model "memoryless" — the epidemic peak has an elegant closed-form solution. Peak prevalence occurs when the susceptible fraction S drops to exactly $1/R_0$, where $R_0$ is the basic reproduction number: the average number of secondary infections produced by a typical case in a fully susceptible population. At that moment, the size of the infectious population is determined by a relationship involving the "final size" of the epidemic — essentially, what fraction of the population eventually gets infected.

This clean result is one of mathematical epidemiology's great achievements. But there's a catch: real infectious periods aren't exponentially distributed. Measles, for instance, has an infectious period that's far more concentrated around its mean than an exponential model would predict. HIV, tuberculosis, and many other diseases have complex, non-exponential infectious trajectories. When epidemiologists want models that reflect biological reality, they often turn to multistage SIR models.

The idea is elegant in its simplicity: instead of a single infectious compartment where recovery is a random lottery (exponential timing), you chain together multiple sequential stages. An individual enters as Stage 1 infectious, progresses to Stage 2, then Stage 3, and so on, until reaching the final stage where recovery occurs. Each progression happens at a fixed rate δ. This "linear chain" structure converts the memoryless exponential distribution into something that, with enough stages, approximates any distribution you want — most commonly, the Erlang distribution, a bell-shaped curve concentrated around its mean.

With $k$ stages, the model becomes an SI(k)R system — $k$ sequential infectious stages leading to removal. The governing equations track each stage separately:

$${\dot{S}}^{(k)}=-\beta S^{(k)}{I}^{(k)}$$$${\dot{I}}_1^{(k)}=\beta S^{(k)}{I}^{(k)}-\delta I_1^{(k)}$$$${\dot{I}}_i^{(k)}=\delta I_{i-1}^{(k)}-\delta I_i^{(k)},i=2,\dots ,k$$$${\dot{R}}^{(k)}=\delta I_k^{(k)}$$

Here $I^{(k)}(t)={\sum }_{i=1}^k{I}_i^{(k)}(t)$ is the total prevalence — the sum across all infectious stages. This is the quantity researchers want to estimate.

The mathematical problem arises because, once you leave the exponential framework, the clean peak formula evaporates. For a multistage model with arbitrary $k$, the prevalence peak has no simple closed form. You can solve the ODEs numerically for specific parameter values, but that requires knowing the full stage-resolved dynamics — information you rarely have in practice.

This is where Tverskoy and colleagues make their entry. Their key insight is to study a quantity called the weighted stage aggregate:

$$V^{(k)}(t)=\sum _{i=1}^k(k-i+1)I_i^{(k)}(t)$$

Think of this as weighting earlier stages more heavily than later ones — the first infectious stage counts as weight $k$, the last as weight 1. This weighted sum, divided by $k$ to normalize it, yields $W^{(k)}(t)=V^{(k)}(t)/k$, which turns out to have a crucial property: its peak location depends only on the susceptible trajectory, not on the full infectious dynamics. Specifically, the weighted peak $W_{\max }^{(k)}$ occurs when $S^{(k)}(t_V)=\delta /(k\beta )$.

From this, the researchers derive an explicit formula for $V_{\max }^{(k)}$ that holds for arbitrary $k$:

$$V_{\max }^{(k)}=k-\frac{\delta }{\beta }-\frac{\delta }{\beta }\ln \left(\frac{\beta k(1-\epsilon )}{\delta }\right)$$

This formula — Proposition 1 in the paper — is the starting point for everything that follows.

The crucial question becomes: what happens as the number of stages $k$ grows large? In the limit of many stages, the Erlang distribution concentrates tightly around its mean, and the multistage model converges to what mathematicians call a delay formulation — a model where infection today contributes to prevalence not instantaneously but with a fixed time lag (the mean infectious period).

But here's where the subtlety enters. There are two different ways to scale the stage-progression rate as $k$ increases, and they lead to dramatically different conclusions about when the factor-two approximation works.

The Two Scalings That Change Everything

Imagine you're building a multistage model with more and more stages. You want the mean infectious period to remain roughly constant — a biologically sensible constraint, since the average time someone is contagious shouldn't change just because you've subdivided the infectious period into more compartments.

The naive approach would be to keep the progression rate $\delta$ fixed as you add stages. This is the naive scaling.

The alternative is Erlang scaling: increase $\delta$ proportionally with $k$. Specifically, set $\delta =k\gamma$, where $\gamma$ is the per-stage progression rate. This preserves the mean infectious period $\tau =1/\gamma$ regardless of how many stages you have. A two-stage model with $\delta =2\gamma$ has the same average duration as a ten-stage model with $\delta =10\gamma$.

These two scalings sound like technical bookkeeping. They are, except they produce entirely different mathematical behavior — and, it turns out, entirely different answers about when the factor-two approximation holds.

Under naive scaling, the paper's Theorem 1 shows that:

$$\frac{W_{\max }^{(k)}}{I_{\max }^{(k)}}=1-\frac{\delta }{\beta }\cdot \frac{\ln (k)}{k}+O(1/k),k\to \infty$$

In plain English: as the number of stages grows large, the weighted peak converges to the prevalence peak. The ratio approaches 1, not 2. Under naive scaling, the factor-two approximation is asymptotically invalid. It's not a matter of being close at finite $k$; in the limit, the proportionality constant is 1.

Under Erlang scaling, the picture flips completely. The infinite-stage limit, established in Theorem 2, shows that the multistage model converges to a delay formulation where prevalence becomes a moving average of incidence — the integral of new infections over the past infectious period:

$$I(t)=\epsilon 1_{\{t<\tau \}}+{\int }_{(t-\tau {)}_{+}}^{t}A(u)du$$

Here $A(u)$ is the incidence curve (new infections per unit time). Meanwhile, the weighted functional $W(t)$ converges to a triangularly weighted moving average — incidence from $\tau$ time units ago gets full weight, and the weight tapers linearly to zero at the present:

$$W(t)=\epsilon {\left(1-\frac{t}{\tau }\right)}_{+}+{\int }_{(t-\tau {)}_{+}}^t\left(1-\frac{t-u}{\tau }\right)A(u)du$$

This is the key structural result. Prevalence is an unweighted average of incidence over the infectious window; the weighted functional is a triangular average that gives more weight to recent infections. And at the peak, where the incidence curve reaches its maximum, these two quantities are related by a factor of approximately 2.

The intuition is geometric. At the peak of the incidence curve, the triangle-shaped weighting function in $W(t)$ assigns roughly equal weight to the peak region as the flat weighting in $I(t)$ — but over a wider effective interval. The ratio converges to 2 in the limit of a sharply peaked incidence curve. This is the mathematical foundation for the factor-two approximation: it's not an arbitrary rule of thumb, but a consequence of the limiting relationship between moving averages and triangularly weighted averages.

The Factor-Two Approximation: When It Works, When It Breaks

The approximation takes the form:

$$I_{\max }\approx 2W_{\max }$$

Tverskoy and colleagues dissect its accuracy with impressive rigor. They show that under Erlang scaling, the approximation's accuracy depends on a single dimensionless parameter $\lambda$ that captures the curvature of the incidence curve near its peak. Roughly speaking:

• Sharp, narrow peaks (large $\lambda$) → the factor-two approximation is accurate
• Broad, flat peaks (small $\lambda$) → the approximation deteriorates

This makes intuitive sense. When incidence spikes sharply and decays quickly, the geometric relationship between the flat moving average (prevalence) and the triangular weighted average (weighted peak) is clean. When the wave is drawn out over weeks, with incidence rising and falling slowly, the effective windows over which these averages are computed look different, and the factor of two no longer holds.

The paper derives an explicit error bound. For a given $R_0$ and curvature parameter $\lambda$, the relative error satisfies:

$$E=\frac{{\hat{I}}_{\max }-I_{\max }}{I_{\max }}\approx \frac{1}{2\lambda }$$

This means the error is inversely proportional to how sharply peaked the incidence wave is. For $\lambda =3$, error is about 17%. For $\lambda =10$, it drops to around 5%.

The researchers also establish necessary conditions for the approximation to be accurate. The theoretical lower bound involves a function $\psi (R_0)$ that reaches a maximum of approximately ${\psi }^{\ast }\approx 0.14$ across the relevant parameter range. This creates a hard ceiling on how well the factor-two approximation can perform: even in the most favorable circumstances, you can't push error below roughly 14% without corrections.

Refined Approximations That Actually Work

This is where the paper's practical contributions crystallize. Tverskoy and colleagues don't just diagnose the problem with the factor-two approximation — they develop a hierarchy of refinements that progressively improve accuracy.

The simplest correction is a first-order approximation (FO) that accounts for the curvature of the incidence curve near its peak:

$${\hat{I}}_{\max }\approx 2W_{\max }\left(1+\frac{1}{2\lambda }\right)$$

The key innovation is that $\lambda$ can be estimated from observable quantities. One approach uses trajectory-based estimation, computing the curvature from the fitted incidence curve. An alternative "plug-in" approach approximates $\lambda$ directly from $R_0$:

$$\lambda \approx \sqrt{\frac{R_0^2+1}{2}}$$

This plug-in estimator is particularly valuable because $R_0$ is often what epidemiologists are trying to estimate in real-time anyway — so you get the curvature correction essentially for free.

For cases requiring higher precision, the paper develops a fully corrected (FC) approximation that accounts for additional higher-order terms, and a large-$\lambda$ approximation (L) that is optimized for scenarios with sharply peaked incidence curves.

The numerical results are striking. In simulations across a broad range of $R_0$ values (1.1 to 8) and stage counts (k = 2 to 100), the simple factor-two approximation shows relative errors that can exceed 40% in some regimes. The fully corrected plug-in approximation brings this down to under 10% across nearly the entire range. The large-$\lambda$ one-step correction performs even better for sharply peaked waves, with errors under 5% in most cases.

Approximation Error Reduction by Method

Label	Value
Simple (S)	40
FO	18
FC	8
L(0)	6
L(1)	5

Relative error magnitude across approximation methods (largest error at each R0, τ=10)

The chart above illustrates the dramatic improvement in accuracy. Under the factor-two approximation alone, errors can range from approximately -35% to +40% depending on $R_0$ — meaning you might underestimate peak prevalence by a third or overestimate it by nearly the same amount. The fully corrected approach (FC, shown in blue) reduces this to roughly ±10% across nearly the entire parameter range. The large-$\lambda$ correction (green) achieves even better performance, staying within ±5% for moderate to high $R_0$ values where epidemic waves tend to be more sharply peaked.

These aren't just marginal improvements. Underestimation of peak prevalence by 30% could mean insufficient ICU beds. Overestimation by the same amount could mean billions in wasted preparedness spending. Getting within 5% of the true peak transforms the utility of these models for public health decision-making.

Error Stability Across Stage Counts

Label	Value
k=2	12
k=5	9
k=10	7
k=20	6
k=50	5

Error magnitude of best approximation at each k, R0=2.5

This second chart shows how approximation accuracy evolves as the number of infectious stages $k$ increases. The horizontal axis represents the finite-stage SI(k)R model, with $k=1$ equivalent to the classical exponential SIR model and larger $k$ approximating more realistic infectious-period distributions.

Two patterns emerge. First, the simple approximation (orange) shows error that can actually grow with $k$ — a troubling finding for anyone assuming that adding more biological realism via more stages would automatically improve model fidelity. Second, the fully corrected plug-in approximation (blue) and the large-$\lambda$ one-step correction (green) both show stable, low error across the entire range of $k$ values tested.

The practical implication is significant: the corrections work regardless of how many stages you assume for the infectious period. Whether you're using a simple single-exponential model ($k=1$) or a more elaborate ten-stage model ($k=10$), the refined approximations consistently outperform the factor-two rule. This matters because different diseases call for different stage structures — and practitioners shouldn't have to choose between biological realism and mathematical accuracy.

Why This Changes Things

Let's step back from the equations and ask: what does this actually mean for the world of epidemic modeling?

The factor-two approximation — multiply the weighted peak by two to get prevalence peak — has been used heuristically for years. Tverskoy and colleagues' paper provides the first rigorous mathematical justification for when it works, when it doesn't, and why.

The justification turns out to be tied to Erlang scaling: the assumption that as you add stages to your model, you scale the progression rate to preserve the mean infectious period. This is, in some sense, the "right" way to think about multistage models — it ensures that the model approximates something biologically meaningful regardless of how many stages you use. And under this scaling, the factor-two approximation emerges naturally from the limiting relationship between unweighted and triangularly weighted moving averages.

But the approximation isn't universally valid. It works best when epidemic waves are sharply peaked — that is, when incidence rises quickly, hits a sharp maximum, and declines rapidly. This corresponds to moderate-to-high $R_0$ values (roughly 2 to 4) and highly contagious diseases with short generation times. In these regimes, the factor-two approximation achieves errors of 10% or better without any correction.

The approximation deteriorates as waves become flatter and more drawn-out. For slow-burning epidemics with gentle peaks, the factor of two can significantly underestimate prevalence. This is exactly the regime where health systems need the most lead time — and where the stakes of getting the peak wrong are highest.

The refined approximations address this gap systematically. By incorporating the curvature of the incidence curve, they capture information that the simple factor-two rule discards. The curvature parameter $\lambda$ is a measure of how sharply peaked the wave is — and the paper shows it can be estimated from $R_0$ alone, without requiring additional fitting or trajectory reconstruction.

This has practical implications for real-time modeling. In the early stages of an outbreak, when surveillance data is noisy and incomplete, epidemiologists often work with rough estimates of $R_0$ and final epidemic size. The factor-two approximation might provide a reasonable first estimate. But as data accumulates and $R_0$ becomes better constrained, the refined approximations can be applied retroactively to improve earlier forecasts — and to understand which interventions, applied when, would have shifted the peak.

What This Opens Up

The paper concludes with several open directions. The analytical framework is developed for deterministic models, but real epidemics are stochastic — individual-level randomness, geographic structure, and heterogeneities in contact patterns all introduce variability that the current analysis doesn't capture. Extending the theory to stochastic settings would bring it closer to the data-generating processes that real epidemiologists work with.

There's also the question of what the curvature parameter $\lambda$ tells us beyond its utility in correction formulas. In the limit of large $R_0$, the paper shows that $\lambda \approx R_0/\sqrt{2}$ — so curvature increases with transmissibility. This makes sense: more transmissible diseases produce sharper, more explosive peaks. But whether $\lambda$ has an independent biological interpretation, or can be linked to other observable quantities like generation time or serial interval, remains an open question.

Perhaps most interestingly, the paper points toward connections with renewal equation formulations — age-of-infection models where the force of infection depends on past incidence weighted by an infectiousness profile. The delay formulation that emerges in the infinite-stage limit is essentially a special case of such models, and the weighted-average relationship between prevalence and incidence has natural analogs in the broader renewal-equation literature.

For now, practitioners working on epidemic preparedness have a clear takeaway: the factor-two approximation is a useful starting point, but it's not a universal truth. Where epidemic waves are sharp — as with measles, chickenpox, or early-phase COVID variants — the simple rule performs well. Where waves are drawn out, as with Omicron's prolonged BA.1 wave or the slow decline of some endemic diseases, corrections matter. The curvature-based refinements Tverskoy and colleagues develop offer a principled way to know when to apply which approximation — and how to correct for the errors that remain.

Limitations and Honest Caveats

No paper is complete without acknowledging what it doesn't do.

The analysis is entirely deterministic. Real epidemics involve stochastic extinction, spatial structure, and individual-level heterogeneity that the ODE framework abstracts away. For large populations where the law of large numbers smooths out individual randomness, the deterministic results are likely accurate. But for small populations, emerging outbreaks with few initial cases, or diseases with high variance in infectiousness, stochastic effects could matter.

The curvature corrections require estimating $\lambda$, which in turn requires either observing the full incidence trajectory or assuming a functional form for the infectiousness profile. In real-time modeling, early in an outbreak, this information may not be available. The plug-in approximation using $R_0$ provides a workaround, but $R_0$ itself is often estimated with substantial uncertainty — especially before genomic surveillance or detailed contact tracing data becomes available.

Finally, the multistage framework assumes that infectiousness is constant across all stages and that progression through stages is deterministic. Real pathogens often show stage-dependent infectiousness (e.g., peak viral load in the first few days of symptoms, waning infectiousness as time passes). Incorporating non-uniform infectiousness profiles would enrich the model but also complicate the mathematical analysis.

These are directions for future work, not fatal flaws. The paper has already done the hard part: providing a rigorous theoretical foundation for an empirical approximation, characterizing its regime of validity, and developing practical corrections that work across a broad parameter range.

The Bottom Line

Epidemic peak estimation is one of those foundational problems that, solved well, makes everything downstream easier. Hospital capacity planning, vaccine deployment timing, non-pharmaceutical intervention thresholds — all of these depend on knowing not just that a wave is coming, but how big it will be at its worst.

Tverskoy, Gothard, and Rempala have taken a rule of thumb that epidemiologists have used for years, excavated its mathematical foundations, and shown exactly when it holds and why it fails. Under Erlang scaling — the natural scaling for multistage models that preserve mean infectious period — the factor-two approximation emerges from the relationship between moving averages and triangularly weighted averages of incidence. It's not arbitrary; it's structural.

But structure doesn't guarantee universal accuracy. The approximation works well for sharply peaked epidemic waves, where curvature is high. It deteriorates for flatter, more drawn-out waves, where the geometric relationship between prevalence and weighted stage aggregates breaks down.

The practical contribution is the correction hierarchy. First-order corrections bring error from potentially 30-40% down to 10-15%. Fully corrected approximations push this to under 10% across the board. Large-$\lambda$ optimizations achieve under 5% for the high-curvature regimes where the approximation was already best. These aren't incremental improvements — they're the difference between a model that's useful for planning and one that might actively mislead it.

For public health modelers, the takeaway is concrete: the factor-two approximation is a tool, not a law. Know when to use it, know when to correct it, and know that the correction is straightforward. The curvature of the incidence curve, measurable from the data you already have, is the key quantity that tells you how much correction is needed.

For the rest of us, the takeaway is something like this: the mathematical models that inform pandemic preparedness are more nuanced than the headlines suggest. They're built on assumptions — about scaling, about averaging, about curvature — that have precise conditions for validity. Understanding those conditions, as Tverskoy and colleagues have done, is what separates a useful model from a dangerous one. And in the next pandemic, that understanding may be the difference between a health system that's ready and one that's overwhelmed.