
==== Front
Proc Natl Acad Sci U S A
Proc Natl Acad Sci U S A
PNAS
Proceedings of the National Academy of Sciences of the United States of America
0027-8424
1091-6490
National Academy of Sciences

38277438
202313708
10.1073/pnas.2313708120
research-articleResearch Articleapp-mathApplied Mathematicspop-bioPopulation Biology404
430
Physical Sciences
Applied Mathematics
Biological Sciences
Population Biology
The probability of epidemic burnout in the stochastic SIR model with vital dynamics
Parsons Todd L. a https://orcid.org/0000-0002-2599-8415

Bolker Benjamin M. b c
Dushoff Jonathan b d https://orcid.org/0000-0003-0506-4794

Earn David J. D. earn@math.mcmaster.ca
c d 1 https://orcid.org/0000-0002-7562-1341

aLaboratoire de Probabilités, Statistique et Modélisation, Sorbonne Université, CNRS UMR 8001, Paris 75005, France
bDepartment of Biology, McMaster University, Hamilton, Ontario L8S 4K1, Canada
cDepartment of Mathematics & Statistics, McMaster University, Hamilton, Ontario L8S 4K1, Canada
dMichael G. DeGroote Institute for Infectious Disease Research, McMaster University, Hamilton, Ontario L8S 4K1, Canada
1To whom correspondence may be addressed. Email: earn@math.mcmaster.ca.
Edited by Alan Hastings, University of California, Davis, CA; received August 9, 2023; accepted November 17, 2023

26 1 2024
30 1 2024
26 7 2024
121 5 e231370812009 8 2023
17 11 2023
Copyright © 2024 the Author(s). Published by PNAS.
2024
https://creativecommons.org/licenses/by-nc-nd/4.0/ This article is distributed under Creative Commons Attribution-NonCommercial-NoDerivatives License 4.0 (CC BY-NC-ND).

Significance

If a new pathogen causes a large epidemic, then it might “burn out” before causing a second epidemic. The burnout probability can be estimated from large numbers of computationally intensive simulations, but an easily computable formula for the burnout probability has never been found. Using a conceptually simple approach, we derive such a formula for the standard SIR epidemic model with vital dynamics (host births and deaths). With this formula, we show that the burnout probability is always smaller for diseases with longer infectious periods, but is bimodal with respect to transmissibility (the basic reproduction number). Our analysis shows that the persistence of typical human infectious diseases cannot be explained by births of new susceptibles, clarifying an important epidemiological puzzle.

We present an approach to computing the probability of epidemic “burnout,” i.e., the probability that a newly emergent pathogen will go extinct after a major epidemic. Our analysis is based on the standard stochastic formulation of the Susceptible-Infectious-Removed (SIR) epidemic model including host demography (births and deaths) and corresponds to the standard SIR ordinary differential equations (ODEs) in the infinite population limit. Exploiting a boundary layer approximation to the ODEs and a birth-death process approximation to the stochastic dynamics within the boundary layer, we derive convenient, fully analytical approximations for the burnout probability. We demonstrate—by comparing with computationally demanding individual-based stochastic simulations and with semi-analytical approximations derived previously—that our fully analytical approximations are highly accurate for biologically plausible parameters. We show that the probability of burnout always decreases with increased mean infectious period. However, for typical biological parameters, there is a relevant local minimum in the probability of persistence as a function of the basic reproduction number R0. For the shortest infectious periods, persistence is least likely if R0≈2.57; for longer infectious periods, the minimum point decreases to R0≈2. For typical acute immunizing infections in human populations of realistic size, our analysis of the SIR model shows that burnout is almost certain in a well-mixed population, implying that susceptible recruitment through births is insufficient on its own to explain disease persistence.

epidemics
stochastic processes
SIR model
extinction
Gouvernement du Canada | Natural Sciences and Engineering Research Council of Canada (NSERC) 501100000038 RGPIN-2021-04068 Benjamin M. BolkerJonathan DushoffDavid J D Earn Gouvernement du Canada | Natural Sciences and Engineering Research Council of Canada (NSERC) 501100000038 RGPIN-2016-06493 Benjamin M. BolkerJonathan DushoffDavid J D Earn Gouvernement du Canada | Natural Sciences and Engineering Research Council of Canada (NSERC) 501100000038 RGPIN-2016-05488 Benjamin M. BolkerJonathan DushoffDavid J D Earn
==== Body
pmcIt is well known that solutions of the standard ordinary differential equations (ODEs) describing a Susceptible-Infectious-Removed (SIR) epidemic with host births and deaths (aka “vital dynamics” or “demography”) eventually converge on a globally asymptomatically stable equilibrium (1). Approach to the endemic equilibrium (EE) typically occurs via damped oscillations, motivating the use of the SIR model with demography as a basis for models of observed recurrent epidemics of childhood infections such as measles (2–6). For many biologically reasonable parameter values and population sizes, however, the troughs of these oscillations pass through infectious-host densities corresponding to a small fraction of an individual—the so-called “atto-fox problem” (7)—calling into question the appropriateness of the deterministic SIR model.

Here, we estimate the probability that a pathogen disappears at the end of a major epidemic in a stochastic individual-based SIR model, in a population of finite size. In the large population limit, the densities of each type (S, I, R) are asymptotically deterministic and governed by the standard SIR ODEs (8). We will refer to pathogen extinction soon after introduction as fizzle, whereas if the pathogen escapes fizzle, we will refer to extinction at the end of a major epidemic as epidemic burnout,* following the terminology of Dushoff (9). We will say that the pathogen persists if it has a subsequent epidemic wave, although it is worth mentioning that we always expect eventual extinction in a stochastic model with a finite population (10); the time to extinction of a pathogen that has survived to a state near the endemic equilibrium is considered in, e.g., refs. 11–13. Fig. 1 shows sample paths of the proportion of infectious individuals for the stochastic SIR model (together with the trajectory obtained from the ODE), illustrating fizzle, burnout, and persistence.

Fig. 1. Sample paths of the stochastic SIR model (Fig. 2) and the ODE (Eq. 5) showing fizzle, burnout, and persistence. (A) The frequencies of susceptible, infectious, and removed individuals in the ODE (symbols indicate the critical points of the curve of the corresponding colour). Dashed lines indicate the endemic equilibrium (Eq. 8) of the deterministic model (Eq. 5). (B) The proportion of infectious individuals as a function of time. The boundary layer—inside which we approximate the stochastic dynamics with a birth-death process—is shaded in yellow, and the first point at which the deterministic trajectory enters the boundary layer is indicated with a heavy yellow dot. (C) Probability density of the time to extinction, estimated from 106 realizations of the stochastic process. The vertical lines show the time τδ (Eq. 76) for which the probability of fizzle after τδ is less than δ (the lines correspond to δ=10−4 and 10−6). (D and E) Trajectories in the susceptible-infectious phase plane with the nullclines; the vertical scale is linear in (D) and logarithmic in (E). The thin black curve is the boundary of the deterministically accessible region, defined by S+I=1. Small yellow dots along trajectories are spaced by one time unit (the mean infectious period).

The problem of epidemic burnout has been of ongoing interest (4, 14, 15), e.g.,

“The question ‘will the agent go extinct after the first outbreak?’ cannot be answered within the context of a deterministic description. So we would like to be able to switch back to a stochastic description at the end of the epidemic outbreak. While it is well known how to calculate the probability of extinction from a branching process in a constant environment..., it seems difficult to do so when environmental quality (from the point of view of the agent, i.e., the presence of susceptibles!) is improving linearly at a certain rate.” (14, p. 42)

and has been previously approached via perturbation methods (16, 17) and by hybrid analytical-numerical approaches (18):

van Herwaarden (16) (henceforth vanH) starts from a large population diffusion approximation to the Markov chain formulation of the SIR model (Model). Under the assumption that the individual mortality rate is low, a highly accurate approximation to the solution of the infinite-population limit SIR ODEs is obtained, which is in turn used to estimate the point of entry to a boundary layer where the number of infectious individuals is very small. In the boundary layer, the backward equation for the diffusion approximation† is tractable and is used to obtain an analytical approximation to the burnout probability [vanH, Eq. (5.13)], which requires the numerical evaluation of an integral. It is, to quote Diekmann and Heesterbeek (14, p. 42), “an ingenious piece of work,” although it is challenging to interpret for non-experts.

By contrast, Meerson and Sasorov (17) (henceforth MS) retain the discrete population model. They estimate the probability of extinction as the probability of reaching the state with only one infective individual (weighted by the expected number of returns to this state‡) times the probability that a single infective recovers before transmitting to any other individuals. They approximate this probability by the product of the expected total time (summed over multiple returns) in the state with a single infectious individual and the disease recovery rate (which is the rate of going extinct given that there is only one infectious individual). The time in the single-infective state is characterized by linear equations obtained by integrating the forward equations (see, e.g., ref. 19, §14.2) for all transient states over all time, for which an approximate solution is found via a WKB ansatz (see, e.g., ref. 22, Chapter 10) in the large population limit (i.e., a diffusion approximation is introduced implicitly). Under these assumptions, the burnout probability is shown to decay exponentially in the population size, with a constant of proportionality that is approximated analytically in the parameter regime where the initial exponential growth rate of infectious individuals greatly exceeds the per capita turnover rate (equivalent to β−γ≫μ in our formulation below). While providing coarser estimates than vanH, this approach yields a deterministic approximation to the most probable trajectory to pathogen extinction via a Hamiltonian formalism (see, e.g., refs. 23 and 24, Exercise 5.7.36). Like vanH, the approximation of MS involves an integral that cannot be evaluated analytically and presents a non-trivial numerical problem due to singularities in the integrand.

More recently, after identifying discrepancies between the analytical results of vanH and MS and the results of simulations, especially at smaller values of the expected population size n, (18) introduced a computational approach that scales as O(n2). As in the previous approaches, Ballard et al. (18) use the solution of the SIR ODEs—now evaluated numerically and summed with a higher-order Gaussian correction (8, Theorem 11.2.3)—to identify the point of entry into a boundary layer, where a simplified form of the Markov chain is then simulated to estimate the probability of burnout.

The approximations of vanH and MS are summarized in §2.3 of ref. 18. We compare the performance of these approximations with that of an analytical approximation that we have derived in the spirit of the quote from ref. 14 above. Like vanH and ref. 18, we use the SIR ODEs to approximate the stochastic SIR trajectories outside a boundary layer. Then, inside the boundary layer, we use a time-inhomogeneous birth-and-death process that approximates the true stochastic dynamics more accurately than the diffusion approximation of vanH (in Boundary Layer Independent Estimates, we obtain the expression from vanH as an approximation to ours). Our approach is simpler and more intuitive than the diffusion approximation, and—in contrast to all previous work—we obtain fully analytical expressions that are numerically stable and can be computed without recourse to numerical evaluation of integrals. Our approach yields expressions for the probability of persistence after any number of epidemic waves and is also more amenable to generalizations than singular perturbation analysis of diffusion approximations; indeed, while we do not discuss the matter in detail here, the boundary-layer diffusions of vanH correspond to large population approximations for the branching processes we consider here (similar to limits in refs. 25 and 26).

Approach and Analysis

Model.

We consider the spread of an infectious disease in a discrete population in which births balance deaths on average, so there is a well-defined expected population size n. We consider a sequence of models indexed by n, and for the nth model denote by Sn(t), In(t) and Rn(t) the numbers of individuals at time t who are susceptible, infectious, and removed, respectively. The total population size is[1] Nn(t)=Sn(t)+In(t)+Rn(t).

Births and immigration of new susceptible individuals occur at constant rate μn, while deaths occur at per capita rate μ, independent of disease status. Thus, at every time t, we have [2] E[Nn(t)]=n,

where the expectation is taken over realizations of the stochastic process. Infectious individuals recover at rate γ, and new infections occur according to the law of mass action in a well-mixed population, i.e., at rate[3] βSn(t)In(t)n.

Since the demographic and epidemiological rates depend only on the state of the system at the current time, our sequence is an ensemble of Markov chain models (indexed by the expected total population size n).

Following a common convention in probability theory, we use upper case for functions and lower case for indices and the values of functions at a given time. We index the functions by expected population size because we need to consider the limit of the sequence of functions as n→∞, whereas we use subscripts on function values to specify time, e.g., s0=Sn(0).

The model structure is indicated in a compartmental transfer diagram in Fig. 2, and the nature and rates of each type of event are summarized in Table 1.

Fig. 2. Compartmental model for an SIR epidemic with vital dynamics. Labels on the arrows correspond to individual jump rates between states. For simplicity, the model is defined so that births/immigrations on average balance deaths, so that the expected total population size (E[Nn(t)]=n) is fixed.

Table 1. Event types in the stochastic SIR model

Event type	Rate	Transitions	
Birth/immigration	μNn	Sn→Sn+1	
Transmission	βSnIn/Nn	Sn→Sn−1, In→In+1	
Recovery	γIn	In→In−1, Rn→Rn+1	
Susceptible death	μSn	Sn→Sn−1	
Infectious death	μIn	In→In−1	
Removed death	μRn	Rn→Rn−1	

Deterministic Approximation.

In the limit of large population size, the stochastic SIR model (Fig. 2 and Table 1) is well-approximated by deterministic ODEs. More precisely, writing[4] Xn=Snn,Yn=Inn,Zn=Rnn,

in the limit n→∞, the frequencies (Xn(t),Yn(t),Zn(t)) converge [almost surely on finite time intervals (8)] to the solution (X(t),Y(t),Z(t)) of the ODEs, [5a] dXdt=μ(1−X)−βXY,

[5b] dYdt=(βX−γ−μ)Y,

[5c] dZdt=γY−μZ.

Formally, to make this connection, one must be careful to have a sensible relationship between the initial conditions for the stochastic processes and the initial conditions for the ODEs. For example, given an initial state (X(0),Y(0),Z(0)) for the ODEs, if one takes[6] (Xn(0),Yn(0),Zn(0))=1n(nX(0),nY(0),nZ(0)),

then the theorem applies. More generally, one must choose initial conditions (Xn(0),Yn(0),Zn(0)) for the stochastic processes such that the limits limn→∞Xn(0), etc. exist, and one must take these limits as initial conditions for the ODEs (see Theorem 11.2.1 in ref. 8, p. 456); Example B on p. 453 of Ethier & Kurtz (8) illustrates how the SIR model without demography relates to the hypotheses of the theorem, and Chapter 5 in ref. 27 provides a pedagogical introduction to Kurtz’s results in the context of epidemic models.

The trajectories of the deterministic SIR model (Eq. 5) always converge to a globally asymptotically stable (GAS) equilibrium point, which can be shown via a combination of the Poincaré Bendixson Theorem and Dulac’s criterion (28) or via a Lyapunov function (29). The nature of the asymptotic state is determined by the basic reproduction number (the expected total number of new infections caused by a single infective individual introduced into a naive population),[7] R0=βγ+μ.

If R0≤1, then the GAS fixed point is the disease free equilibrium, (x,y)=(1,0), whereas if R0>1, then all solutions converge—either via damped oscillations or monotonically—to an endemic equilibrium,[8] (x⋆,y⋆)=1R0,ε1−1R0,

where[9] ε=μγ+μ

gives the mean infectious period as a fraction of the mean host lifetime. Our analysis requires that ε is small but not too small (1n≪ε≪1), which is true for a wide variety of common acute immunizing infections (Table 2). The upper bound (ε≪1) is essential so we can justify perturbation expansions in ε. The lower bound (1n≪ε) is equivalent to nε≫1, which ensures that the number of infectives at equilibrium (ny⋆) is substantially greater than 1 (from Eq. 8, ny⋆∼nε). boundarylayerThe ODEs continue to provide a good approximation to the epidemic dynamics until the prevalence y (the proportion of hosts that are infectious) becomes small; we take “small” to mean that y is less than the equilibrium prevalence y⋆ (Eq. 8). Thus, we take the boundary layer—within which the dynamics must be treated stochastically—to be the region of the phase plane where y<y⋆ (in Boundary Layer Independent Estimates, we also give approximations independent of the specific choice of boundary layer).

Table 2. Representative parameters for acute immunizing infections (and HIV for comparison)

Disease	R0	Tlat
[days]	Tinf
[days]	ε×103	Source	
Measles	17	8	5	0.71	(4)	
Pertussis	17	8	14	1.2	(4)	
Mumps	12	15	6	1.1	(4)	
Chickenpox	11	10	5	0.82	(4)	
COVID-19 (Delta)	6.8	5.8	14	1.1	(39)	
Rubella	6.5	10	7	0.93	(4)	
Scarlet fever	5.5	1.5	18	1	(4)	
Smallpox	4.5	15	7	1.2	(40)	
COVID-19 (ancestral)	3	3.7	14	0.97	(39)	
HIV	2.2	87	270	19	(41)	
Influenza (1918)	1.8	2	2.5	0.25	(4, 42)	
Ebola	1.6	9.3	7	0.89	(43)	
Pneumonic plague	1.3	4.3	2.5	0.37	(44)	
The basic reproduction number (R0), mean latent period (Tlat), and mean infectious period (Tinf) are taken from the cited sources. The dimensionless parameter ε is defined in Eq. 9 in terms of the recovery rate (γ) and birth-death rate (μ) in the SIR model. We associate 1/γ with the mean generation interval of the SEIR model, i.e., 1/γ=Tlat+Tinf (45, 46), set μ=0.02/year to mimic human birth and death rates, and compute ε=μ/(γ+μ). Where original sources present a range, we have listed the midpoint. Many of the estimates come from Anderson and May (4) [R0 is taken from Table 4.1 (4, p. 70); the mean latent and infectious periods come from Table 3.1 (4, p. 31)]. All the diseases listed in this table are shown in Fig. 5.

The need to analyze the dynamics differently within the boundary layer is especially clear if we consider the introduction of a single infectious individual into a fully susceptible population. If R0>1, then in the ODE system (Eq. 5) Y(t) will deterministically increase, whereas in the stochastic model, Yn(t) will fizzle with probability 1/R0 (30); i.e., the ODE (Eq. 5) fails to capture the dynamics of the stochastic model (Fig. 2) when there are few infectives. We therefore use a birth-and-death process to approximate the dynamics of the number of infectious hosts when that number is small (in contrast, susceptibles can be assumed to remain sufficiently abundant that we can always use the deterministic approximation X(t)).

Birth-and-Death Process Heuristic.

New infections occur at rate[10] βSn(t)nIn(t)=βXn(t)In(t)≈βX(t)In(t),

while the number of infectious hosts decreases by one due to recovery or death at rate[11] (γ+μ)In(t).

When there are few infectious hosts (In(t)<ny⋆), we approximate In(t) by a linear birth and death process with time-inhomogeneous per capita rates b(t) and d(t), where[12] b(t)=βX(t),

[13] d(t)=γ+μ.

Note that when X(t) equals x⋆ (the classical herd immunity threshold), b(t) = d(t), and the birth and death process transitions from subcritical to supercritical. Unlike in models without demography, the birth of new susceptible individuals ensures that a population will eventually cross the herd immunity threshold. Therefore, even if the number of infectious hosts initially declines it can eventually grow exponentially, if the infection survives until X(t)>x⋆.

We can estimate the survival probability for this branching process, and thus, the persistence probability, using the following result.

Theorem 1 [Kendall (31)]. Let K(t) be a birth and death process with time-inhomogeneous per-capita birth rate b(t) and death rate d(t). The probability of eventual extinction starting from one individual at time 0 is

[14] q=1+1∫0∞e−∫0t[b(s)−d(s)]dsd(t)dt−1.

The extinction probability starting from k individuals is qk.

Consequently, the probability of indefinite persistence (a branching process will either go extinct or grow indefinitely), starting from k individuals at time 0, is [15] P{K(∞)>0}=1−qk.

To complete our persistence probability estimate, we need an expression for the proportion susceptible at time t (X(t) in Eq. 12). As suggested visually by the example shown in Fig. 1, inside the boundary layer (y<y⋆), both the deterministic and the stochastic trajectories spend most of their time at prevalences much lower than y⋆ [note the log scale in the subfigures (B) and (E)]. Consequently, we can approximate X(t) by solving (Eq. 5) with Y(0)=0. Thus, we set[16] dXdt≈μ(1−X),

and solve this approximate equation as if it were exact to obtain[17] X(x0,t)≈1−(1−x0)e−μt.

Here, x0 is the fraction susceptible at the initial time t=0, and we write X(x0,t) to emphasize the dependence on the initial state. We also write q(x0) for the value of q in Eq. 14 obtained by taking b(t)=βX(x0,t).

We first apply this branching process approximation to a population at the disease-free equilibrium (DFE). Thus, we set x0=1 in (Eq. 17), which yields X(1,t)≡1; hence, we have a time-homogeneous branching process in this case, and the integral in (Eq. 14) is easily evaluated and yields q(1)=1R0=x⋆. Considering a small number of initially infective individuals, In(0)=k, we recover the classical expression for the establishment probability (30), that is, the probability that the pathogen does not fizzle:[18] pk=1−x⋆k.

We now use Kendall’s q (Eq. 14) to compute the burnout probability. Assuming that the pathogen does not fizzle, the number of infectious hosts will rapidly exceed ny⋆ individuals,§ at which point the densities of both susceptible and infective hosts are well approximated by the ODEs (Eq. 5). To proceed, we need a formula for the fraction of hosts that are susceptible when the trajectory enters the boundary layer at the end of an epidemic; we denote this fraction xin to emphasize that it refers to the susceptible proportion upon entry into the boundary layer (the point (xin,y⋆) is indicated by a heavy yellow dot in Fig. 1). In ref. 34, assuming ε is small,¶ we derive an approximate expression for the fraction susceptible, X(y,xi), as a function of the fraction infectious (y) and the initial fraction susceptible (xi). Using that approximation, we have[19] xin=X(y⋆,xi)≈−x⋆W0(−R0xie−R0(xi−y⋆))    +εeR0y⋆(E1(R0y⋆)−E1(R0y¯0)).

Here, W0 denotes the principal branch of the Lambert W-function# (35), E1(x)=∫x∞e−ttdt is the exponential integral function (36, §8.2.1) and y¯0 is the peak prevalence in the limit ε→0 (i.e., μ→0), i.e., the maximum fraction infectious in the SIR model without vital dynamics,[21] y¯0=xi−x⋆(1+ln(xi/x⋆)).

(See e.g., ref. 37 for a derivation of y¯0.) Taking xi=1 corresponds to the invasion of a novel pathogen into an epidemiologically naive population (i.e., at the DFE). Later (Subsequent Epidemic Waves), we give an iterative scheme for xi,j, an “effective initial fraction susceptible” that—substituted for xi in Eqs. 19 and 21—gives the fraction susceptible at the end of the jth epidemic wave after invasion at the DFE. We compare our approximation of xin for xi=1 to the value obtained by numerically integrating the SIR ODEs (Eq. 5) in Fig. 3 and discuss its domain of applicability below (The Domain of Applicability of the Approximation (Eq. 19) to xin).

Fig. 3. Susceptible proportion (xin) upon entry into the boundary layer (y<y⋆). (A) xin as a function of R0 (Eq. 7). (B) xin as a function of ε (Eq. 9). The exact value of xin (obtained by numerically solving the SIR ODEs (Eq. 5)) is shown with solid curves, our approximation (Eq. 19) is shown with dashed curves, and the approximation of vanH is shown with dotted curves. Based on Eq. 42, the minimum R0 for which our approximation of xin (Eq. 19) is valid is ≈e2ε (i.e., 1.02027 for ε=0.01 and 1.0020027 for ε=0.001).

If we now take t=0 to be the end of a major epidemic, i.e., the time when the infectious host density falls below y⋆ and x0=xin, then the density of infectious hosts is small, and the density of susceptible hosts is well approximated by X(xin,t) (we are preparing a rigorous treatment of these results; here, we will content ourselves with showing that our analytical results closely match the results of individual-based simulations). We can thus estimate the conditional burnout probability—i.e., the probability of burnout conditional on not fizzling—by[22] q(xin)ny⋆.

and the conditional persistence probability by

[23] 1−q(xin)ny⋆.

Below (Computing the Epidemic Burnout Probability), we compute an exact expression for q(xin), [24a] q(xin)=1+εz−aezga,z−1

[24b] where  z=R0ε(1−xin),

[24c] and  a=R0ε1−x⋆.

Here, g denotes the lower incomplete gamma function‖ ( 36, §8.2.1); we use the nonstandard notation g to avoid confusion with our recovery rate parameter γ. Below (Asymptotics for Small ε), we derive an approximation for q(xin) that is extremely accurate for small values of ε:[25] q(xin)≈1+12πε(R0−1)(az)aez−a−1.

We emphasize that this expression is elementary and numerically stable.

Thus, the burnout probability—i.e., the probability of not fizzling (Eq. 18) but disappearing after an epidemic—is[26] pkq(xin)ny⋆,

where n is the expected total population size, y⋆ is the equilibrium prevalence (Eq. 8), q is the probability of eventual extinction (under post-epidemic conditions) starting from one infectious individual (Eq. 14), and k is the initial number of infectious individuals. Our exact expression for q(xin) is given in (Eq. 24). Similarly, the persistence probability—i.e., the probability of not fizzling (pk) and then not burning out after a first epidemic (Eq. 23)—is[27] P1(R0,ε,n,k)=pk1−q(xin)ny⋆.

More generally, the probability of persisting beyond the mth epidemic wave is[28] Pm(R0,ε,n,k)=pk∏j=1m1−q(xin,j)ny⋆,

where[29] xin,j=X(y⋆,xi,j)

(see Eq. 19 and Subsequent Epidemic Waves). For biologically reasonable values of ε, R0, and n, we find that the difference between P1(R0,ε,n,k) and Pm(R0,ε,n,k) is negligible, because q(xin,j)≪1 for j≥2. Intuitively, because the troughs between epidemics get shallower and shallower, an invading disease that survives burnout is almost certain to persist through many more cycles.

Thus, in Results, we focus on burnout after the initial epidemic when a novel disease invades a fully susceptible population. There, we use our accurate, numerically stable, and computationally efficient approximation for q(xin) (Eq. 25), obtained via Eqs. 19 and 37a, to compute the probability of burnout.

Results

Fig. 4 shows that our analytical approximation for the persistence probability (Eq. 27) agrees very well with the same probability estimated from large numbers of simulations. The probability is shown as a function of the basic reproduction number (R0) with fixed mean infectious period (ε=0.01). The panels differ only in the underlying expected population size (ranging from n=104 to 107). For each value of R0, the simulation-based persistence probability was estimated from 107 individual-based stochastic realizations of the model (Fig. 2 and Table 1). Note that ε=0.01 corresponds to an infectious period that is 1% of the average lifetime, far longer than is realistic for most acute immunizing infections; however, our approximation only improves for smaller ε. We use ε=0.01 in Fig. 4 so that discrepancies between the simulations and analytical results are visible.

Fig. 4. Persistence probability as a function of the basic reproduction number R0, for population sizes ranging from n=104 to 107. The vertical scale is linear in the left column and logarithmic in the right column; the horizontal scale is logarithmic (in R0−1) in all panels. (The horizontal axis range is from R0−1=164=0.015625 to 64, but our approximation is valid only for R0−1≳0.02027; see Eq. 42.) The initial state is (Sn(0),In(0),Rn(0))=(n−1,1,0). The mean infectious period as a fraction of mean lifetime is ε=0.01, which is unrealistically long for most infections (Table 2), but the agreement between the analytical approximation (Eq. 27) and numerical simulations (Stochastic simulation algorithm) is better for smaller ε (we use a large value of ε so that discrepancies are visible). In addition to our analytical approximation (Eq. 27), we show the semi-analytical approximations of Meerson and Sasorov [MS (17)] and van Herwaarden [vanH (16)]. The thin red curve shows the probability of not fizzling, 1−1R0.

Our simple approximation for Kendall’s q (Eq. 25) allows us to easily and quickly explore the conditional and unconditional probability of pathogen extinction across the entire range of biologically plausible values of R0 and ε. Fig. 5 shows a contour plot of the persistence probability (this graph would have required years of computer time to produce from simulations). As was observed previously (18, 38), Fig. 5 indicates that the burnout probability is non-monotone in R0 for ε≲0.016. In this range of ε, the probability of persistence is lowest for basic reproduction numbers in the range 2≲R0≲2.57, and increases rapidly with increasing R0. Below (Maximizing the Probability of Burnout), we compute a linear approximation to the value of R0 at which the persistence probability is minimized. The upper limit of 2.57 for the persistence-minimizing range of R0 is the limit as ε→0 in Eq. 63; Fig. 5 shows that this linear approximation performs very well over the range where the persistence probability is non-monotonic. Less intuitively, the persistence probability increases for small R0 (below the red curve in Fig. 5) as R0 decreases to one. We note, however, that except for very large expected population size n, the secondary peak in the persistence probability—which occurs for 1<R0≲2—remains small (cf.Fig. 4), except for pathogens with extremely long infectious periods. Fig. 5 also suggests that for fixed R0, the probability of persistence always increases with increasing ε, which we confirm analytically below (The Burnout Probability is a Decreasing Function of ε). Note that β=R0(γ+μ)=R0γ1−ε, so varying ε while holding R0 constant simultaneously varies the infectious period and contact rate by a factor of O(ε).

Fig. 5. Probability of persistence after a large epidemic (P1, Eq. 27) as a function of basic reproduction number (R0) and mean infectious period as a proportion of mean lifetime (ε), for population size n=106. The initial state is assumed to be a single individual introduced into a fully susceptible population (In(0)=k=1, Sn(0)=n−k). Positions for the red dots for infectious diseases of humans are from Table 2 [to avoid text overlap, measles is shifted up by 1 to 18, pertussis down by 1 to 16, and COVID-19 (Delta) up by 0.6 to 7.4]. The solid red curve shows the local minimum of persistence probability, and the dotted red line shows the analytical approximation (Eq. 63) to the local minimum.

Discussion

The problem of infectious disease persistence following a major epidemic (4, p. 20); (14, p. 42); (47, p. 451); 9, 15, 38) is important for identifying characteristics of pathogens that can successfully invade, and is related to the notion of a “critical community size” required for a disease to persist in the long term (30).

Given sufficient computing resources, it is possible to estimate the persistence probability for a given model from large numbers of stochastic, individual-based simulations. The gray curves in Fig. 4 show this probability estimated from simulations of the SIR model. Fig. 4 also shows the probability estimated using previous analytical methods (16, 17) (blue and orange curves) and our approximation (black curves). All three analytical approaches yield similar results,** and differences in the estimated probabilities can be seen only on a logarithmic scale in the limit as R0→1+ (e.g., for R0≲1.05 in Fig. 4), where all of these approximations†† are technically invalid: in a stochastic, finite population model, as R0→1+ there is no phase during which the deterministic model is a good approximation, and the distinction between fizzle, burnout, and fadeout breaks down (48). Analysis of the limit R0→1+ could improve understanding of the process of eradication as the magnitude of control measures is increased, but for the burnout problem on which we focus here, the limit R0→1+ is of limited interest.

While our approximation agrees closely with previous work (16, 17) for ranges of R0 that are biologically relevant, there are several important theoretical and practical advantages of our approach; our analysis

is simpler and easier to understand, since it is based directly on the underlying stochastic process rather than on a diffusion approximation (and is consequently easier to apply to models that are more complex than the SIR model considered here);

yields fully analytical approximations that are numerically stable, unlike the previous analytical approaches (16, 17), which depend on non-trivial numerical integrations with singular integrands;

predicts the persistence probability after an arbitrary number of epidemic waves.

We expand on these points below.

We have obtained useful analytical estimates (Eqs. 24, 25, 27, and 28) of the SIR epidemic burnout and persistence probabilities in a well-mixed population, via a hybrid use of ODEs when prevalence is high and time-dependent branching processes when prevalence is low. As noted after Eq. 28, the probability of burning out in each subsequent epidemic trough after persisting through the first is negligibly small for the SIR model.

Our time-dependent branching process approach (Birth-and-Death Process Heuristic) also yields analytical results that are more amenable to computation than previous approximations (16, 17). Our application of Laplace’s method to approximate the integral in Kendall’s q (Eq. 14) is particularly useful. Eq. 25 for the conditional burnout probability provides a fully analytical formula—not requiring the numerical evaluation of integrals as in previous approaches (16, 17)—that can be evaluated without numerical instabilities and agrees very well with numerical simulations across a wide range of biologically plausible values of R0 and ε. The convenience and speed of our simple analytical expression for the persistence probability (Eq. 27) also allows us to obtain results for larger population sizes than are tractable via hybrid numerical methods (18) and facilitates efficient exploration of more of the parameter space (though with less accuracy at smaller population sizes).

As is suggested visually by Fig. 5, and proved below (The Burnout Probability is a Decreasing Function of ε), the persistence probability increases with infectious period (ε) across all values of R0. For any given infectious period, one viable life history strategy for persistence is a high R0 (dark gray shading in Fig. 5). In addition to this high R0 strategy, for a limited range of longer infectious periods (ε≲0.016), there is a second life-history strategy that promotes persistence: R0 close to but greater than one. We use our analytical results to compute a linear approximation to the value of R0>1 at which the burnout probability is maximized (see Eq. 63 in Maximizing the Probability of Burnout). This approximation shows excellent agreement with the numerical results over the range of ε for which the secondary peak exists and the burnout probability is numerically distinguishable from 1 (in Fig. 5, the dotted red curve is the approximation and the solid red curve is the numerically computed exact value). Intriguingly, with the exception of the ancestral strain of SARS-CoV-2—which has been replaced by variants with much higher R0—the endemic infectious diseases of humans listed in Table 2 roughly divide into high and low R0 strategies.

These life history strategies can be interpreted in terms of the herd immunity threshold, x=x⋆=1R0, i.e., the minimum proportion susceptible at which the epidemic can grow from a small number of infections. When R0 is large, the herd immunity threshold x⋆ is low, allowing the fraction susceptible to rapidly reach the threshold. When R0 is low, there is a larger reservoir of susceptible hosts at the end of the first major epidemic, which reduces the wait until the herd immunity threshold is crossed. In either case, a longer infectious period (larger ε) allows the pathogen to “wait out” the period of herd immunity. This non-monotonicity of the burnout probability as a function of R0 was previously observed (18, 38), and the maximum burnout probability was conjectured to occur for R0≈3 (38) or R0=2 (18). We have shown that, in fact, the value of R0 at which the probability is maximized is a decreasing function of ε (solid red curve in Fig. 5). The probability-maximizing R0 varies from R0≃2.57 for ε→0 (Eq. 63) to R0≃2 for ε≃0.016; for larger ε, the persistence probability increases monotonically with R0.

These results also have evolutionary implications: reduced virulence may be associated with longer infectious periods (e.g., if fewer hosts die while infectious), thereby reducing the probability of burnout. This suggests a mechanism—distinct from the population genetics/weak selection arguments presented in ref. 32—that could explain how in finite populations, natural selection may favour strains with reduced virulence while maximizing R0: strains that achieve higher R0 by increasing their infectious period are more likely to persist than those that achieve higher R0 by increasing transmissibility. A model of multi-strain competition is necessary to test this hypothesis (we consider the special case of selection for vaccine escape in ref. 49).

While the qualitative inferences we have made from analysis of the stochastic SIR model are suggestive of general processes, and—as we have observed above—could have interesting implications, further research is needed to determine whether they really do generalize broadly. Most acute immunizing infections afflicting human populations have short infectious periods and moderate R0 values, and with these constraints, our analysis of the stochastic SIR model indicates that extinction of the pathogen at the end of the first major epidemic is almost certain in a well-mixed population.

Fig. 5 makes clear that the stochastic SIR model is insufficient on its own to explain pathogen persistence; it is essential to consider additional mechanisms, e.g., waning immunity or antigenic evolution resulting in effective loss of host immunity (50), rescue effects in a meta-population (51, 52), long-lived carrier infections (see ref. 53 for a recent survey), or zoonotic reservoirs (50, 54).

Multi-type or non-Markovian birth-and-death processes (55, 56), combined with more complicated compartmental models or renewal equation models with more general generation intervals (46) may allow our approach to be extended to models incorporating, e.g., latent periods and asymptomatic and carrier infections, or greater or lesser variability in infectious periods. A more difficult problem is to consider pathogen persistence in a meta-community of linked sites (52, 57), or other structured populations, rather than a well-mixed population. Smaller local community sizes tend to make local extinction more likely, whereas asynchrony in epidemic dynamics could allow pathogens to reinvade following a local extinction (58). Are these processes adequate to plausibly explain the persistence of pathogens? Is the existence of low/high R0 strategies generic, or an artifact of the SIR compartmental model? Are longer infectious periods always favourable for pathogen persistence? These questions suggest avenues for future work.

Materials and Methods

Computing the Epidemic Burnout Probability.

To apply Kendall’s q (Eq. 14) to the problem of epidemic burnout, we need to compute the integral [30a] I(xin)=∫0∞e−∫0t[βX(xin,s)−(γ+μ)]ds(γ+μ)dt

[30b] =∫0∞e−∫0τ[R0X(xin,σγ+μ)−1]dσdτ,

where, in the second line, we use the mean duration of infection (1/(γ+μ)) as the time unit and write σ=(γ+μ)s, τ=(γ+μ)t. Recalling Eq. 17, we can write[31] X(σ)≡Xxin,σγ+μ=1−(1−xin)e−εσ,

and hence[32] X′(σ)=ε(1−xin)e−εσ=ε(1−X(σ)).

Now, to evaluate the inner integral in Eq. 30a, we make a change of variables, using x=X(σ) as the variable of integration:[33] ∫0τR0Xxin,σγ+μ−1dσ=∫X(0)X(τ)[R0x−1]1dxdσdx=∫xinX(τ)R0x−1ε(1−x)dx=−R0ε(X(τ)−xin)−R0ε1−x⋆ln1−X(τ)1−xin.

Changing variables in a similar way, we have [34a] ∫0Te−∫0τ[R0X(xin,σγ+μ)−1]dσdτ

[34b] =∫0TeR0ε(X(τ)−xin)+R0ε1−x⋆ln1−X(τ)1−xindτ

[34c] =∫xinX(T)eR0ε(x−xin)1−x1−xinR0ε1−x⋆dxε(1−x).

We are interested in the probability of ultimate extinction, which corresponds to taking the limit as T→∞, or, equivalently, X(T)→1, giving us [35a] I(xin)=∫xin1eR0ε(x−xin)1−x1−xinR0ε1−x⋆dxε(1−x)

[35b] =1εeR0ε(1−xin)R0ε(1−xin)−R0ε1−x⋆×∫0R0ε(1−xin)e−xxR0ε1−x⋆−1dx

[35c] =1εeR0ε(1−xin)R0ε(1−xin)−R0ε1−x⋆×gR0ε(1−x⋆),R0ε(1−xin),

where we recall g denotes the lower incomplete gamma function. Eq. 24 follows immediately.

Asymptotics for Small ε.

We may also write I(xin) (Eq. 35) as [36a] I(xin)=1ε∫xin111−xeR0εx−xin+(1−x⋆)ln1−x1−xindx

[36b] =1ε∫xin1h(x)eϕ(x)εdx,

for h(x)=11−x and ϕ(x)=R0x−xin+(1−x⋆)ln(1−x1−xin). Assuming ε is small, we can apply Laplace’s method (22, §6.4): provided xin≤x⋆, ϕ(x) has its maximum at x=x⋆, so [37a] I(xin)∼1ε2πε|ϕ″(x⋆)|h(x⋆)eϕ(x⋆)ε

[37b] =2πε(R0−1)eR0εx⋆−xin1−x⋆1−xinR0ε1−x⋆,

yielding Eq. 25.

Remark 1: Note that, since xin<x⋆,[38] 0 <−R0∫xfx⋆ln(1−t)dt =R0((x⋆−xin)+(1−x⋆)ln(1−x⋆)−(1−xin)ln(1−xin)) <R0((x⋆−xin)+(1−x⋆)ln(1−x⋆)−(1−x⋆)ln(1−xin)) =ϕ(x⋆),

so the Laplace approximation and thus the original integral (Eq. 35a) are both exponentially large in ε−1.

Subsequent Epidemic Waves.

In ref. 34, we derive an iterative scheme to compute “effective initial conditions” for every epidemic wave following initial disease invasion. Writing xi,j for the fraction susceptible at the start of the jth epidemic wave, we find our trajectory approximations agree very closely with the “exact” value obtained by solving the SIR ODEs (Eq. 5) numerically, starting from the DFE.

Setting xi,1=1, we iteratively obtain y¯0,j and xi,j+1 (Eq. 19) from xi,j by computing [39a] xf,j=−x⋆W0(E(−xi,j/x⋆)),

[39b] y¯0,j=xi,j−x⋆1+ln(xi,j/x⋆).

[39c] xi,j+1=1+(1−x⋆)W0(E(−1−xf,j1−x⋆)).

Note that xf,j and y¯0,j, are the final fraction susceptible (i.e., when the pathogen has gone extinct) and maximal fraction infectious, respectively, for the SIR model without vital dynamics (ε=0) with initial condition (xi,j,0+).

The Domain of Applicability of the Approximation (Eq. 19) to xin.

The refined trajectory approximation that yields Eq. 19 is derived in ref. 34 under the assumption that R0 is large. Despite this, we find that the approximation to xin obtained from it (Eq. 19) performs very well for all but values of R0 very close to 1 or very large values of ε>0 (Fig. 3). In particular, W0(x) is undefined for x<−e−1, so we must have[40] −R0e−R01−y⋆>−e−1,

or, expanding and rearranging using Eq. 8,[41] ε<1−lnR0R0−1=1+x⋆lnx⋆1−x⋆.

Alternately, we can find an approximate lower bound for R0,[42] R0>e2ε,

by observing that for any x≥1, 1−lnxx−1≤12lnx. To derive this latter inequality, note that both sides approach a limit of 0 as x→1, whereas[43] ddx1−lnxx−1−12lnx=1(x−1)2lnx−x2−12x.

Again, lnx−x2−12x vanishes at x=1, whereas[44] ddxlnx−x2−12x=−(x−1)22x≤0,

so lnx−x2−12x≤0 for x≥1, and thus, ddx1−lnxx−1−12lnx≤0 also, proving the desired inequality.

Boundary Layer Independent Estimates.

Thus far, we have computed the burnout probability via a specific, but arbitrary choice of boundary layer y⋆, and explicit solutions for xin, the fraction susceptible when first entering the boundary layer under the ODE approximation (Eq. 5). Here, we consider an alternative approach, using results from refs. 34 and 59 to implicitly characterize xin. In conjunction with Eq. 35, this allows us—at the cost of a small loss of precision—to give expressions for the extinction and persistence probabilities that are independent of the precise choice of threshold, provided the threshold is O(ε). In addition to being of interest in and of themselves, we use them to compute the value of R0 maximizing the burnout probability (The R0 Maximizing the Probability of Burnout) and also to show how one derives the result of ref. 16 as an approximation to Eq. 24 (Boundary Layer Independent Estimates).

In refs. 34 and 59, we use the method of matched asymptotic expansions (60, 61) to derive analytical approximations to the phase-plane trajectories of the SIR model with vital dynamics, i.e., expressions Y(x) and X(y) expressing the density of infectious hosts as a function of the density of susceptible hosts and vice versa. In the boundary layer, we obtain lowest- and first-order approximations to Y(x): the lowest-order approximation (34) is[45] Y(x)≈y¯01−xf1−xR0ε(1−x⋆)eR0ε(xf−x),

whereas the refined estimate (59) is[46] Y(x)≈1xf−1x⋆−xf1−xf1−xR0ε(1−x⋆)×e−R0ε(x−xf)+1xf−1−1Yxf1(1),

where[47] xf=−x⋆W0(−R0e−R0)

is the final size of the SIR epidemic without vital dynamics (62) and [48a] Yxf1(1)=∫xf1(x⋆t−1)1u−111−u+x⋆lnu−1xf−11t−xfdt

[48b] ≈1xf−x⋆x⋆−xflnxf−x⋆x⋆−xf1xf−1.

A very closely related expression (using μ rather than ε as the small parameter) is derived in ref. 16.

Recalling that, Y(xin)=y⋆, evaluating either of Eq. 45 or 46 at x=xin gives us a relation between xin, xf, and y⋆. From the former (Eq. 45), we have[49] 11−xinR0ε(1−x⋆)e−R0εxin≈y⋆y¯011−xfR0ε(1−x⋆)e−R0εxf,

whereas the latter (Eq. 46) gives us[50] 11−xinR0ε(1−x⋆)e−R0εxin≈y⋆(1−xf)x⋆xf−111−xfR0ε(1−x⋆)×e−R0εxf+1xf−1−1Yxf1(1).

Substituting Eq. 49 into the integrand in Eq. 35a and proceeding as above gives [51a] I(xin)≈y⋆y¯0∫xin11−x1−xfR0ε(1−x⋆)eR0ε(x−xf)dxε(1−x)

[51b] =1εy⋆y¯0R0ε(1−xf)−R0ε(1−x⋆)eR0ε(1−xf)×gR0ε(1−x⋆),R0ε(1−xin).

Set z′=R0ε(1−xf). Then, az (Eqs. 24b and 24c) and az′ are fixed, while as ε→0, a→∞ and g(a,z)∼Γ(a)−zae−z (see ref. 36, §8.11.6) and similarly for g(a,z′). Thus [52] (z′)−aez′(g(a,z′)−g(a,z))∼zz′aez′−z−1.

and the error in replacing xin by xf in the incomplete gamma function in Eq. 51b is equal to [53a] 1εy⋆y¯01−xin1−xfR0ε1−x⋆e−R0ε(xf−xin)

[53b]  =1εy⋆y¯0eR0ε1−x⋆ln1−xin1−xf−(xf−xin)

[53c]  =1εy⋆y¯0eR0ε1−x⋆ln1+xf−xin1−xf−(xf−xin)

[53d]  =1εy⋆y¯0eR0ε(xf−xin)xf−x⋆1−xf+O(ε2).

Both xf−xin and y⋆ are O(ε), whereas xf−x⋆1−xf is O(1), so this error is O(1). Thus, in absolute terms, the error is not small. However, as we observed above, I(xin) is exponentially large in ε−1, so the error is negligible relative to this leading term [indeed, replacing the incomplete gamma function by ΓR0ε(1−x⋆) produces a similarly negligible error]. We can also replace xin by xf in the Laplace approximation with negligible error:[54] I(xin)≈y⋆y¯02πε(R0−1)az′aez′−a.

Similarly, repeating the same argument using the higher-order expression, Eq. 50, gives [55a] I(xin)≈1εy⋆(1−xf)1R0xf−1(z′)−a×ez′+1xf−1−1Yxf1(1)g(a,z′)

[55b]  ≈y⋆(1−xf)1R0xf−12πε(R0−1)az′a×ez′−a+1xf−1−1Yxf1(1).

Now, we recall from Eqs. 14 and 22 that the burnout probability is [56a] q(xin)ny⋆=1+1I(xin)−ny⋆

[56b]  =e−ny⋆ln1+1I(xin)

[56c]  ≈e−ny⋆I(xin).

Above, we showed that I(xin) is exponentially large in ε−1 as ε→0 and, thus, that the error in making the last approximation (Eq. 56c) is exponentially small.

Substituting any of the expressions Eq. 51b, 54, 55a, or 55b for I(xin) in Eq. 56c, we see that the factors y⋆ cancel, giving us an approximate expression for the burnout probability that does not depend on the specific choice of threshold, only upon its order of magnitude, ε: [57a] q(xin)ny⋆≈e−nεy¯0(z′)−aez′g(a,z′)

[57b]  ≈e−ny¯0ε(R0−1)2πz′aaez′−a,

or [58a] q(xin)ny⋆≈ exp−nε(1−xf)1R0xf−1z′−aez′+1xf−1−1Yxf1(1)ga,z′

[58b]  ≈e−n(1−xf)1R0xf−1ε(R0−1)2πz′aaez′−a−1xf−1−1Yxf1(1),

respectively.

Remark 2: If in Eq. 58a we approximate ga,z′ by Γ(a) (i.e., if we approximate the integral up to z′≫1 by the integral over the whole real line, introducing an error of O(ε)), we obtain an expression for the burnout probability equivalent to that from ref. 16 (up to minor differences resulting from using different small parameters, μ and ε).

The R0 Maximizing the Probability of Burnout.

Using the simplified expression for the burnout probability (Eq. 57b), we can obtain an approximation to the value of R0 that maximizes the probability of burnout linear in ε which is highly accurate across the range of values of ε for which the burnout probability is non-monotone. Eq. 57b is minimized when[59] y¯0ε(R0−1)2π1−xf1−x⋆R0ε1−x⋆eR0εxf−x⋆

is maximized, or equivalently, when its partial derivative with respect to R0 is equal to zero. Computing the partial derivative and collecting terms of like order in ε, we seek R0 such that[60] 1εlnR0+W0−R0e−R0R0−1−1R0+lnR0R0(R0−1−lnR0)+R0−12=0.

An analytical closed-form solution does not appear to exist, but one can use a formal asymptotic series expansion R0=∑j=0∞rjεj to obtain a polynomial approximation in ε to arbitrarily large degree (here, we content ourselves with a linear approximation). Substituting this series into Eq. 60 and collecting terms of order ε−1 and order one, we obtain[61] lnr0+W0−r0e−r0r0−1−1r0=0,

[62] −1+(r02−r0+1)W0−r0e−r0r02(r0−1)(1+W0−r0e−r0)r1+r0−12+lnr0r0(r0−1−lnr0)=0.

We may solve Eq. 61 by Newton iteration to find the unique root r0=2.572629848, which we use to solve Eq. 62 to find r1=−27.71866282, giving us the linear approximation [63] arg maxR0>1  q(xin)ny⋆≈2.572629848−27.71866282ε.

We compare this linear approximation to the numerically determined minimum in Fig. 5.

The Burnout Probability is a Decreasing Function of ε.

In what follows, we show that ∂q(xin)∂ε≤0, from which we conclude that q(xin) is decreasing as ε increases, for all values of R0. Using Eq. 24, we have [64] ∂q(xin)∂ε=−q(xin)21ezz−ag(a,z)−ε∂∂ε[ezz−ag(a,z)]e2zz−2ag(a,z)2.

The first term in the large brackets on the right-hand side is always positive, so the result follows if one can show that ε∂∂ε[ezz−ag(a,z)]≤0. Applying the chain rule gives [65a] ε∂∂ε[ezz−ag(a,z)]

[65b]  =ε∂z∂ε∂∂z[ezz−ag(a,z)]+ε∂a∂ε∂∂a[ezz−ag(a,z)]

[65c]  =−(z+R0∂xin∂ε)∂∂z[ezz−ag(a,z)]−a∂∂a[ezz−ag(a,z)].

Recalling that g(a,z)=∫0zta−1e−tdt, the latter is equal to [66] −(z+R0∂xin∂ε)ezz−ag(a,z)1−az+1z−aezz−a∫0zta−1e−tlnt dt−g(a,z)lnz.

Integrating by parts in the rightmost term, this becomes [67] −(z+R0∂xin∂ε)ezz−ag(a,z)1−az+1z−aezz−a∫0zg(a,t)tdt=−R0∂xin∂εezz−ag(a,z)1−az+1z−ezz−a+1g(a,z)−1−ezz−a∫0z(atg(a,t)−azg(a,z))dt.

Now, g(a,z)≥0, whereas az=1−x⋆1−xin≤1, since xin<x⋆, so 1−az≥0 and, since g(a,z) is an increasing function of z, we have [68] ∫0z(atg(a,t)−azg(a,z))dt≥∫0zat−azg(a,t)dt≥0.

Thus, provided ∂xin∂ε≥0, ε∂∂ε[ezz−ag(a,z)]≤0, as required.

Finally, from Eq. 19, we see that[69] ∂xin∂ε=limε→0eR0y⋆(E1(R0y⋆)−E1(R0y¯0))≥0,

since y⋆≤y¯0 and E1(x) is a decreasing function of x.

Simulations

Simulations.

Stochastic simulation algorithm.

Exact realizations of the stochastic SIR model (Fig. 2 and Table 1) can be obtained using the standard Gillespie algorithm (63, 64). If we denote the various event rates ai (e.g., a1=μn, etc.), then the total event rate is a=∑iai. The time to the next event is drawn from an exponential distribution with mean 1/a, and the event is taken to be of type i with probability ai/a. This algorithm scales with expected population size n and is prohibitively slow when running large numbers of simulations with n≳105. We therefore used the adaptive τ-leaping approximation (65), as implemented in the adaptivetau R package (66). The key idea in this approach is to identify, at any point of the simulation, a time τ over which the various event rates can be considered approximately constant, and then determine the number of events of each type that can be expected over this time interval. We then “leap forward” by time τ rather than treating events individually.

Estimating the Required Number of Simulations.

To determine the number of simulations required to estimate the epidemic burnout probability to a given accuracy, we use the central limit theorem. Suppose we run m independent simulations. Let [70] 1i=1if theithsimulation ends in burnout,and0otherwise.

Then, the law of large numbers (67, §6) tells us that [71] limm→∞1m∑i=1m1i=E[11]=q,

where q is Kendall’s q (Eq. 14). Consequently, qm=1m∑i=1m1i is an unbiased estimator (67, p. 483) of q. Let Δ=q−qm be the error in our estimates. Then, the central limit theorem (67, §27) tells us that mΔ=1m∑i=1m(1i−q) converges to a normal distribution with the same variance as 11−q, i.e., mΔ converges in distribution to a normal random variable with variance [72a] σ2=E[(11−q)2]=E[112−2q11+q2]

[72b]  =(12·q+02·(1−q))−2q(1·q+0·(1−q))+q2

[72c]  =q(1−q)≤14,

where the inequality in (Eq. 72c) follows because 0≤q≤1. In particular, for large m, the expected squared error is E[Δ2]≲14m, and thus, to have E[Δ2]≤δ, we perform at least m=⌈14δ⌉ runs.

Fizzle vs. Epidemic Burnout.

To efficiently distinguish fizzles from epidemic burnout, we use Eq. 15 to estimate a time τδ (measured in units of the mean infectious period 1/(γ+μ)) such that the probability is less than δ that, starting from k infectious individuals and xi=1−kn, a sample path in which infective individuals are still present at time τδ eventually fizzles. Let Tk be the (random) time of fizzle starting from k individuals. Then, using our birth-and-death process approximation, [73a] P{Tk>t}     =P{In(t)>0∣In(0)=k}

[73b] ≈1−1+1∫0te−∫0s[β−(γ+μ)]du(γ+μ)ds−k

[73c] =1−1+11R0−11−e−(β−γ−μ)t−k.

Now, because fizzle is not a certainty, [74] limt→∞P{Tk>t}=1−x⋆k>0.

To determine τδ, we condition on eventual fizzle to estimate its time of occurrence:

[75a] P{Tk>t∣Tk<∞}=P{Tk>t}−P{Tk=∞}P{Tk<∞}

[75b] ≈1R0k−1+11R0−11−e−(β−γ−μ)t−k1R0k

[75c] =1−1+1R0R0−11−e−(β−γ−μ)t−k.

Solving for P{Tk>τδ∣Tk<∞}=δ yields[76] τδ=1R0−1ln(1−δ)−1k−1R0(1−δ)−1k−1.

Choosing a suitably small δ, we assume that any sample path in which infective individuals are still present at τδ will not fizzle.

This project was partially supported by the CNRS International Emerging Actions (IEA) grant “Structured Populations, Epidemics, and Control Strategies (SPECS).” D.J.D.E., J.D., and B.M.B. were supported by the Natural Sciences and Engineering Research Council of Canada (NSERC). We thank David Champredon, Michelle deJonge, Sarah Drohan, Karsten Hempel, Chai Molina, Irena Papst, and Dora Rosati for their contributions to the preliminary work that led to ref. 68.

Author contributions

All authors contributed to the conception of the research; T.L.P. and D.J.D.E. carried out the analysis, simulations, and figure creation, with input from B.M.B. and J.D.; T.L.P. and D.J.D.E. wrote the manuscript; and all authors revised the manuscript and approved the final version.

Competing interests

The authors declare no competing interest.

Data, Materials, and Software Availability

All study data are included in the main text. Our open-source R package, which we used to create our figures, is available at https://github.com/davidearn/burnout.

This article is a PNAS Direct Submission.

*While “fade-out” (or “fadeout”) is commonly used to describe this extinction, e.g., ref. 4, §2.3, we find it conceptually useful to follow ref. 9 in distinguishing between extinction after a first major epidemic versus that occurring after multiple epidemics, and reserve the term fadeout for the latter.

†See, e.g., refs. 19 or 20 for a discussion of the forward and backward diffusion equations; ref. 21 is an excellent introduction to boundary-layer methods for Markov chains.

‡In practice, there is negligible probability of returning to the state with one infective after an excursion to a state with many infectives.

§More precisely, for any y<y¯0 (Eq. 21), conditional on not fizzling, the probability that In(t) hits 0 before hitting yn is exponentially small in n with exponential rate depending on y [for a rigorous demonstration see ref. 32 (Supplementary Information §8.2); ref. 33 gives explicit higher-order terms for the SIS model].

¶In ref. 34, we use ϵ=ε/R0 rather than ε as the small parameter, because using ϵ leads to simpler expressions (see, e.g., (Eq. 24) below). Here, however, we analyze the dependence of our expressions on the epidemiologically relevant parameters R0 and ε and have re-written expressions from ref. 34 accordingly.

#If E(z)=zez, Lambert’s W-function W(z) (35) solves the “left-sided” inverse relation E(W(z))=z. This equation has countably many solutions, each corresponding to branches Wi of the W-function; we will need the two real branches, W0, which maps [−1e,∞) to [−1,∞), and W−1, which maps [−1e,0) to (−∞,−1]. For these two branches, Wi is a partial “right-sided” inverse function for E(z): [20] W−1(E(z))=zifz≤−1W0(E(z))=zifz≥−1.

‖ g(a,z)=∫0zta−1e−tdt is proportional to the cumulative distribution function for the gamma distribution. We use this fact to compute g(a,z) accurately in our burnout R package, mentioned in footnote **.

**We have implemented all three approximations in an open-source R package, which we used to create our figures. The package is available at https://github.com/davidearn/burnout.

††Differences between our approximation and those of refs. 16 and 17 as R0→1+ arise at least in part because they use μ rather than ε as the small parameter, and consequently predict persistence for β/γ>1 rather than β/(γ+μ)>1.
==== Refs
1 F. Brauer, C. Castillo-Chavez, Z. Feng, Mathematical Models in Epidemiology (Springer, 2019), vol. 32 .
2 M. S. Bartlett, Measles periodicity and community size. J. R. Stat. Soc., Ser. A 120 , 48–70 (1957).
3 W. London, J. A. Yorke, Recurrent outbreaks of measles, chickenpox and mumps. I. Seasonal variation in contact rates. Am. J. Epidemiol. 98 , 453–468 (1973).4767622
4 R. M. Anderson, R. M. May, Infectious Diseases of Humans: Dynamics and Control (Oxford University Press, Oxford, UK, 1991).
5 D. J. D. Earn, P. Rohani, B. M. Bolker, B. T. Grenfell, A simple model for complex dynamical transitions in epidemics. Science 287 , 667–670 (2000).10650003
6 K. Hempel, D. J. D. Earn, A century of transitions in New York City’s measles dynamics. J. R. Soc. Lond. Interface 12 , 20150024 (2015).
7 D. Mollison, Dependence of epidemic and population velocities on basic parameters. Math. Biosci. 107 , 255–287 (1991).1806118
8 S. N. Ethier, T. G. Kurtz, Markov Processes: Characterization and Convergence (John Wiley and Sons, New York, NY, 1986).
9 J. Dushoff, “Incorporating stochasticity in simple models of disease spread” in Modeling Paradigms and Analysis of Disease Transmission Models, A. B. Gumel, S. Lenhart, Eds. (American Mathematical Society, 2010), vol. 75.
10 P. Jagers, Stabilities and instabilities in population dynamics. J. Appl. Prob. 29 , 770–780 (1992).
11 O. A. Van Herwaarden, J. Grasman, Stochastic epidemics: Major outbreaks and the duration of the endemic period. J. Math. Biol. 33 , 581–601 (1995).7608639
12 I. Nåsell, On the time to extinction in recurrent epidemics. J. R. Stat. Soc., Ser. B 61 , 309–330 (1999).
13 A. Kamenev, B. Meerson, Extinction of an infectious disease: A large fluctuation in a nonequilibrium system. Phys. Rev. E 77 , 061107 (2008).
14 O. Diekmann, J. A. P. Heesterbeek, Mathematical Epidemiology of Infectious Diseases: Model Building, Analysis and Interpretation (John Wiley & Sons, New York, NY, 2000).
15 T. Britton , Five challenges for stochastic epidemic models involving global transmission. Epidemics 10 , 54–57 (2015).25843384
16 O. A. van Herwaarden, Stochastic epidemics: The probability of extinction of an infectious disease at the end of a major outbreak. J. Math. Biol. 35 , 793–813 (1997).9269737
17 B. Meerson, P. V. Sasorov, WKB theory of epidemic fade-out in stochastic populations. Phys. Rev. E 80 , 041130 (2009).
18 P. G. Ballard, N. G. Bean, J. V. Ross, The probability of epidemic fade-out is non-monotonic in transmission rate for the Markovian SIR model with demography. J. Theor. Biol. 393 , 170–178 (2016).26796227
19 S. Karlin, H. M. Taylor, A Second Course in Stochastic Processes (Academic Press, San Diego, CA, 1981).
20 C. W. Gardiner, Handbook of Stochastic Methods for Physics, Chemistry and the Natural Sciences (Springer, Berlin/Heidelberg, Germany, 2004).
21 J. Grasman, O. A. van Herwaarden, Asymptotic Methods for the Fokker–Planck Equation and the Exit Problem in Applications (Springer, Berlin/Heidelberg, Germany, 1999).
22 C. M. Bender, S. A. Orszag, Advanced Mathematical Methods for Scientists and Engineers (McGraw-Hill, New York, NY, 1978).
23 R. Graham, T. Tél, Existence of a potential for dissipative dynamical systems. Phys. Rev. Lett. 52 , 9–12 (1984).
24 A. Dembo, O. Zeitouni, Large Deviations Techniques and Applications, Applications of Mathematics (Springer, New York, NY, ed. 2, 1998), vol. 38.
25 W. Feller, “Diffusion processes in genetics” in Proceedings on the Second Berkeley Symposium on Mathematical Statistics and Probability, July 31–August 12, 1950, J. Neyman, Ed. (University of California Press, Berkeley, CA, 1951), pp. 227–246.
26 J. Lamperti, The limit of a sequence of branching processes. Z. Wahrsch. Verw. Gebiete 7 , 271–288 (1967).
27 H. Andersson, T. Britton, Stochastic Epidemic Models and their Statistical Analysis (Springer, New York, NY, 2000), vol. 151 .
28 H. W. Hethcote, H. W. Stech, P. van den Driessche, “Periodicity and stability in epidemic models: A survey” in Differential Equations and Applications in Ecology Epidemics, and Population Problems, S. N. Busenberg, K. L. Cooke Eds. (Elsevier, 1981), pp. 65–82.
29 A. Korobeinikov, G. C. Wake, Lyapunov functions and global stability for SIR, SIRS, and SIS epidemiological models. Appl. Math. Lett. 15 , 955–960 (2002).
30 M. S. Bartlett, “Deterministic and stochastic models for recurrent epidemics” in Proceedings of the Third Berkeley Symposium on Mathematical Statistics and Probability (1956), vol. 4, pp. 81–108.
31 D. G. Kendall, On the generalized “Birth-and-Death’’ process. Ann. Math. Stat. 19 , 1–15 (1948).
32 T. L. Parsons, A. Lambert, T. Day, S. Gandon, Pathogen evolution in finite populations: Slow and steady spreads the best. J. Royal Soc. Interface 15 , 20180135 (2018).
33 T. L. Parsons, Invasion probabilities, hitting times, and some fluctuation theory for the stochastic logistic process. J. Math. Biol. 77 , 1193–1231 (2018).29947947
34 T. L. Parsons, D. J. D. Earn, Analytical approximations for the phase plane trajectories of the SIR model with vital dynamics. [Preprint] (2023). https://cnrs.hal.science/hal-04178969 (Accessed 8 August 2023).
35 R. M. Corless, G. H. Gonnet, D. E. G. Hare, D. J. Jeffrey, D. E. Knuth, On the Lambert W function. Adv. Comput. Math. 5 , 329–359 (1996).
36 F. W. J. Olver, D. W. Lozier, R. F. Boisvert, C. W. Clark, NIST Handbook of Mathematical Functions (Cambridge University Press, 2010).
37 H. W. Hethcote, The mathematics of infectious diseases. SIAM Rev. 42 , 599–653 (2000).
38 D. J. D. Earn et al., “The puzzling persistence of invading pathogens” in Poster presentation at 2014 Ecology and Evolution of Infectious Diseases Conference “Multi-Scale Mechanisms of Disease Emergence and Control” (Colorado State University, 2014).
39 Z. Levine, D. J. D. Earn, Face masking and COVID-19: Potential effects of variolation on transmission dynamics. J. R. Soc. Interface 19 , 20210781 (2022).35506215
40 O. Krylova, “Predicting epidemiological transitions in infectious disease dynamics: Smallpox in historic London (1664–1930),” PhD, McMaster University, Canada (2011).
41 T. D. Hollingsworth, R. M. Anderson, C. Fraser, HIV-1 transmission, by stage of infection. J. Infect. Dis. 198 , 687–693 (2008).18662132
42 C. E. Mills, J. M. Robins, M. Lipsitch, Transmissibility of 1918 pandemic influenza. Nature 432 , 904–906 (2004).15602562
43 Z. Y. Wong, C. M. Bui, A. A. Chughtai, C. R. Macintyre, A systematic review of early modelling studies of Ebola virus disease in West Africa. Epidemiol. Infect. 145 , 1069–1094 (2017).28166851
44 R. Gani, S. Leach, Epidemiologic determinants for modeling pneumonic plague outbreaks. Emerg. Infect. Dis. 10 , 608–614 (2004).15200849
45 O. Krylova, D. J. D. Earn, Effects of the infectious period distribution on predicted transitions in childhood disease dynamics. J. R. Soc. Lond. Interface 10 , 20130098 (2013).
46 D. Champredon, J. Dushoff, D. J. D. Earn, Equivalence of the Erlang SEIR epidemic model and the renewal equation. SIAM J. Appl. Math. 78 , 3258–3278 (2018).
47 A. A. King, S. Shrestha, E. T. Harvill, O. N. Bjørnstad, Evolution of acute infections and the invasion-persistence trade-off. Am. Nat. 173 , 446–455 (2009).19231966
48 I. Nåsell, “The threshold concept in stochastic epidemic and endemic models” in Epidemic Models: Their Structure and Relation to Data, D. Mollison, Ed. (Cambridge University Press, 1995), pp. 71–83.
49 S. Gandon, A. Lambert, T. Day, T. L. Parsons, The speed of vaccination rollout and the risk of pathogen adaptation. medRxiv [Preprint] (2022). 10.1101/2022.08.01.22278283 (Accessed 13 August 2022).
50 D. J. D. Earn, J. Dushoff, S. A. Levin, Ecology and evolution of the flu. Trends Ecol. Evol. 17 , 334–340 (2002).
51 B. T. Grenfell, J. Harwood, (Meta)population dynamics of infectious diseases. Trends Ecol. Evol. 12 , 395–399 (1997).21238122
52 D. J. D. Earn, P. Rohani, B. T. Grenfell, Persistence, chaos and synchrony in ecology and epidemiology. Proc. R. Soc. Lond. B 265 , 7–10 (1998).
53 K. Hampson, D. Haydon, Persistent pathogens and wildlife reservoirs. Science 374 , 35–36 (2021).34591640
54 D. T. Haydon, S. Cleaveland, L. H. Taylor, M. K. Laurenson, Identifying reservoirs of infection: A conceptual and practical challenge. Emerg. Infect. Dis. 8 , 1468–1473 (2002).12498665
55 C. J. Mode, Multitype Branching Processes: Theory and Applications (Elsevier, New York, NY, 1971).
56 P. Jagers, Branching Processes with Biological Applications (Wiley, London, UK, 1975).
57 C. McCluskey, D. J. D. Earn, Attractivity of coherent manifolds in metapopulation models. J. Math. Biol. 62 , 509–541 (2011).20425115
58 D. J. D. Earn, S. A. Levin, P. Rohani, Coherence and conservation. Science 290 , 1360–1364 (2000).11082064
59 T. L. Parsons, D. J. D. Earn, Refined asymptotic approximations for the phase plane trajectories of the SIR model with vital dynamics. [Preprint] (2023). https://cnrs.hal.science/hal-04178983 (Accessed 8 August 2023).
60 R. E. O’Malley Jr., Singular Perturbation Methods for Ordinary Differential Equations (Springer, 1991), vol. 89 .
61 J. Kevorkian, J. D. Cole, Multiple Scale and Singular Perturbation Methods (Springer-Verlag New York Inc., New York, NY, 1996).
62 J. Ma, D. J. D. Earn, Generality of the final size formula for an epidemic of a newly invading infectious disease. Bull. Math. Biol. 68 , 679–702 (2006).16794950
63 D. T. Gillespie, A general method for numerically simulating the stochastic time evolution of coupled chemical reactions. J. Comput. Phys. 22 , 403–434 (1976).
64 D. T. Gillespie, Exact stochastic simulation of coupled chemical reactions. J. Phys. Chem. 81 , 2340–2361 (1977).
65 D. T. Gillespie, Stochastic simulation of chemical kinetics. Annu. Rev. Phys. Chem. 58 , 35–55 (2007).17037977
66 P. Johnson, adaptivetau: Tau-Leaping Stochastic Simulation. R package v 2.2-3 (2019) https://CRAN.R-project.org/package=adaptivetau.
67 P. Billingsley, Probability and Measure, Wiley Series in Probability and Statistics (John Wiley & Sons, Hoboken, NJ, Anniversary ed., 2012).
68 D. J. D. Earn, B. M. Bolker, J. Dushoff, T. L. Parsons, Burnout: a package for computing infectious disease burnout and persistence probabilities. R package version 0.0.1. Github. https://github.com/davidearn/burnout (Accessed 13 August 2023).
