Duration of a minor epidemic.

Disease outbreaks in stochastic SIR epidemic models are characterized as either minor or major. When ℛ 0 < 1 , all epidemics are minor, whereas if ℛ 0 > 1 , they can be minor or major. In 1955, Whittle derived formulas for the probability of a minor or a major epidemic. A minor epidemic is distinguished from a major one in that a minor epidemic is generally of shorter duration and has substantially fewer cases than a major epidemic. In this investigation, analytical formulas are derived that approximate the probability density, the mean, and the higher-order moments for the duration of a minor epidemic. These analytical results are applicable to minor epidemics in stochastic SIR, SIS, and SIRS models with a single infected class. The probability density for minor epidemics in more complex epidemic models can be computed numerically applying multitype branching processes and the backward Kolmogorov differential equations. When ℛ 0 is close to one, minor epidemics are more common than major epidemics and their duration is significantly longer than when ℛ 0 ≪ 1 or ℛ 0 ≫ 1 .

Introduction

Public health intervention and control strategies are designed to shorten the course of an epidemic. This is often achieved through reduction of the basic reproduction number to a value below the critical threshold of one. However, for close to one, the duration of an epidemic is increased. This property is demonstrated in the following study of the duration of minor epidemics in stochastic epidemic models.

It is well-known for deterministic SIR epidemic models that if , the number of infected cases decreases over time, whereas if , the number of cases increases. However, the outcome differs in stochastic SIR epidemic models. In 1955 Whittle characterized two different types of epidemics in stochastic SIR models as minor or major (Whittle, 1955). Only minor epidemics occur when and either minor or major occur when . Minor epidemics are distinguished from major ones as generally being of shorter duration and with substantially fewer cases. Whittle derived the well-known formulas for the probability of a minor epidemic and for the probability of a major epidemic, given there are initially i infected individuals and . Analytical expressions for the duration of these minor outbreaks are not well known. In this investigation, analytical formulas are derived that approximate the probability density, the mean, and the higher-order moments for a minor epidemic of continuous-time Markov chain (CTMC) SIR, SIS and SIRS epidemic models. In particular, it is shown that the mean duration T of a minor epidemic given iswhere β is the transmission rate, γ is the recovery rate, and .

The probability of a minor or a major epidemic has been studied in more complex disease settings than the CTMC SIR epidemic model by approximating the dynamics with a simple birth-death process, the backward Kolmogorov differential equations, and a multitype branching process approximation (e.g. Allen, 2015; Allen, 2017; Allen & Lahodny, 2012; Allen & van den Driessche, 2013; Athreya & Ney, 1972; Griffiths, 1972; Griffiths, 1973). The probability estimates are generally better for large population sizes and large values of (e.g. Allen, 2015; Allen, 2017; Allen & Lahodny, 2012; Allen & van den Driessche, 2013). Most theoretical research in stochastic models has focused on the duration of major SIR epidemics and quasistationarity in SIS models when conditioned on nonextinction (e.g. Artalejo, 2012; Daley & Gani, 1999; Hernandez-Ceron, Chavez-Casillas, & Feng, 2015; Hernández-Suárez & Castillo-Chavez, 1999; Kryscio & Lefévre, 1989; Norden, 1982; Nåsell, 1996; Nåsell, 1999; Nåsell, 2001; van Doorn, 1991). In this investigation, the emphasis is on the duration of minor epidemics, conditioned on extinction.

In the following section, the simple birth-death process and branching process theory are used to derive the probability density and analytical formulas for the moments for time to extinction. In Section 3, the results from the simple birth-death process are applied to minor epidemics in CTMC SIR and SIS models. Numerical results show good agreement between the simulations and analytical estimates, given the population size is large. The increase in duration when is close to one is demonstrated in these models and in Section 4 as well for the SEIR model with a latent stage.

Simple birth-death process

The simple birth-death process is a well-known approximation for many population processes (see e.g. Allen, 2010; Athreya & Ney, 1972; Bailey, 1975; Daley & Gani, 1999; Novozhilov, Karev, & Koonin, 2006; Sehl, Zhou, Sinsheimer, & Lange, 2011; Whittle, 1955). The mean of the simple birth-death process is an exponential growth model,where B and D are the per capita birth and death rates, respectively. We summarize briefly the simple birth-death process and the expressions that lead to the probability of extinction and the time to extinction.

Let denote the discrete random variable for the population size in a time-homogeneous process with transition probabilities defined asThe infinitesimal transition probabilities for a small period of time can be expressed asThe expressions in (2) lead to the backward Kolmogorov differential equations

An important assumption in the simple birth-death process is that births and deaths occur independently of each other. Therefore, the probability of extinction at time t beginning from i individuals, , can be written as . Derivation of the expression for follows from generating functions and branching process theory (Athreya & Ney, 1972; Bailey, 1975; Daley & Gani, 1999; Harris, 1963):

An outline of this derivation is given in Appendix A. The probability of ultimate extinction follows directly from (4):

If , extinction occurs with probability one, but if , extinction occurs with probability and the population grows without bound with probability . The time to extinction is finite only in the case . If the mean time to extinction is infinite while the probability of extinction is one. If , the mean time to extinction is also infinite. We are interested in the probability of hitting zero and the time it takes to hit zero, that is, population extinction when . Therefore, we condition on the time to extinction being finite and only consider those sample paths that hit zero (see Fig. 1).

Fig. 1

Five sample paths for the simple birth-death process are plotted along with the exponential growth model . Three sample paths have finite time to extinction, whereas two sample paths approach infinity. Close-ups of the sample paths are graphed to time 20 on the right. Parameter values are and , with initial condition . For the simple birth-death process, probability of extinction is and the mean and the standard deviation for time to extinction, conditioned on extinction, are and .

Probability density function for time to extinction

Let denote the time to extinction given . In the case , extinction occurs with probability one, so that the cumulative distribution function (cdf) of is given by in (4). In the case , extinction occurs with probability . Therefore, conditioned on extinction the cdf of is . The cdf with is equivalent to the cdf with by interchanging the roles of B and D, i.e., . Differentiating with respect to t, leads to the pdf for time to extinction:where and . In the case and , has a Pareto type II distribution (also known as a Lomax distribution) with infinite mean and higher-order moments (Lomax, 1954).

Graphs of the pdf for parameter values , and are plotted for initial conditions in Fig. 2. The tails of the distribution become “fat” as B approaches one (see Table 1). Some new analytical expressions for the mean and higher-order moments for the time to extinction are computed in the next section.

Fig. 2

Plots of the pdf for , , and show that the tails of the distribution become “fat” as B approaches one.

Table 1

The probability that time to extinction is greater than five or ten time units given X(0) = 3 and D = 1 (See Fig. 2.).

B	ℙ(T_3>5)	ℙ(T_3>10)
2	0.011	6.81×10^−5
1.5	0.0843	6.75×10^−3
1.25	0.207	0.0518
1.1	0.325	0.143
0.9	0.347	0.156
0.75	0.249	0.0642
0.5	0.123	0.0101
0.1	0.0297	3.33×10^−4

Mean and higher-order moments for time to extinction

The mth moment of can be computed directly from the integralwhere is one of the two forms in (5). A closed form analytical expression for the mean when the initial size is is given byKarlin and Tavaré (Karlin & Tavaré, 1982) derived the preceding expression in the case when applied to a population genetics model with mutation.

In the following cases, we assume that . This assumption is not restrictive since scaling time in (3), , leads to the backward differential equation for expressed in terms of a single parameter . In the following, we keep the notation B instead of recalling that .

Explicit formulas for the moments are expressed in terms of special functions, polylogarithm functions of order n, , . The polylogarithm function of order n can be expressed as a power series in x (Apostol et al., 2010):or defined recursively asFor order one, the polylogarithm has a closed form expression, The first and higher-order moments for the duration until extinction are defined in terms of these special functions.

The mean time to extinction for and iswhere are the harmonic numbers, , etc. If the summation is undefined, then it is set to zero. Extensions to higher-order moments when and arefor . The preceding expressions for the moments agree with (6), when , and .

The higher-order moments , , involve several polylogarithm functions. (The explicit expressions for the case and are given in Appendix B.) Recently, Singh et al. (Singh, Schneider, & Myers, 2014) derived formulas for the moments of but for a birth-death process with killing (Karlin & Tavaré, 1982).

The preceding formulas are used to compute some values for the probability of extinction, , mean , and standard deviation in Table 2 for .

Table 2

The probability of extinction, the mean, and the standard deviation for the time until extinction, conditioned on extinction, for the simple birth-death process with X(0) = 3 and D = 1.

B	ℙ_ext	μ	σ
0.7	1.00	3.50	3.12
0.9	1.00	5.66	6.92
0.95	1.00	7.28	10.93
1.05	0.864	7.05	10.73
1.1	0.751	5.34	6.71
1.5	0.296	2.22	1.92

Application to epidemic models

Near the disease-free equilibrium (DFE), the dynamics of the classic differential equations for SIR, SIS, and SIRS epidemics are approximated by the differential equation for infected individuals (Appendix C):Parameter β is the transmission rate, γ is the recovery rate, and . Setting and , application of the theory from the preceding simple birth-death approximation leads to formulas for the duration of a minor epidemic. It is well-known from Whittle's work in 1955 (Whittle, 1955) that the probability of a minor epidemic, , or a major epidemic, , given , areSee Fig. 3. The formulas for the mean and the standard deviation, equations (6), (8), for the duration of a minor epidemic are not well-known:and for ,

Fig. 3

Plot of the probability of a minor outbreak and a major outbreak as a function of for , .

In the following sections, for large population sizes, numerical computations for the CTMC SIR () and CTMC SIS ) models show good agreement with the preceding formulas.

Stochastic SIR epidemic model

To define the stochastic SIR epidemic model, let , , and be discrete random variables for the number of susceptible, infected (and infectious), and immune individuals in a homogeneously mixed population with constant total population size N. The CTMC SIR model is a bivariate process that depends on two discrete random variables with transition probabilities,In the following numerical examples, .

The probability density function for the duration of a minor epidemic in the CTMC SIR epidemic model (12) is plotted as a probability histogram in Fig. 4. All simulations are run until and the time recorded, . The probability histogram shows a mixture of two densities when : the duration of a minor epidemic with probability and the duration of a major outbreak with probability . The density for is the curve overlayed in red on the probability histograms in Fig. 4. The figures illustrate the long duration until disease eradication for as opposed to .

Fig. 4

Probability histograms, computed from sample paths in the CTMC SIR model (12) show a mixture of two densities when , the probability of a minor or major epidemic. The curve overlayed in red is a graph of the minor epidemic , where is given in equation (5). Parameter values are , , and with initial conditions , , and .

The mean and standard deviation for a minor epidemic are computed numerically in the CTMC SIR model (12) and they are compared to the analytical formulas derived from the simple birth-death process, (10), (11), (20) (Fig. 5). In the numerical simulations of the CTMC SIR epidemic model for , the sample paths for infected individuals are divided into two sets, those that hit zero before reaching a predefined outbreak level (minor epidemic) and those that reach the value (major epidemic). The outbreak level is set at . For , the duration is computed as but for , the duration is . If hits zero before , then it is assumed to be a minor epidemic but if hits before zero, then it is assumed to result in a major epidemic. The values for , mean, and standard deviation computed from the CTMC SIR simulations and the corresponding approximations from the simple birth-death process are given in Table 3 for the case in Appendix D.

Fig. 5

The analytical formulas for the mean and standard deviation, from equations (10), (11), (20) for a minor epidemic are graphed as a function of for three initial conditions . The corresponding mean and standard deviation for duration of a minor epidemic are computed numerically from sample paths of the CTMC SIR model (12) (+for ; x for ; for ). Parameter values are , , and .

Stochastic SIS model

The dynamics of the CTMC SIS model can be expressed in terms of one discrete random variable , where .

Numerical computations are performed for the CTMC SIS model (13) similar to those for the CTMC SIR model. The difference is that in the CTMC SIS model when , all simulations are run until either hits zero or reaches a predefined outbreak level . The level is either the endemic equilibrium value from the deterministic model: or . In the following numerical examples, .

Probability histograms for the duration of minor or major epidemics in the CTMC SIS model (13) are graphed in Fig. 6. For , the duration of a minor epidemic is , whereas for , the duration is . The curve overlayed in red is the density , where is defined in (5).

Fig. 6

Probability histograms computed from sample paths of the CTMC SIS model show a mixture of two densities when , the probability of a minor epidemic and the probability of reaching the endemic equilibrium. The overlay curve in red is a graph of . Parameter values are , and with initial conditions and

Graphed in Fig. 7 are the mean and standard deviation for duration of a minor epidemic. The curves are from the formulas (10), (11), and (20) in the birth-death process. The calculated values for the CTMC SIS model (12) (similar to Fig. 5) assume a predefined outbreak level of . For , the duration of a minor epidemic is computed as but for , the duration is computed as . If hits 30 before zero, then is assumed to reach an endemic level, but if hits zero before 30, it is assumed to be a minor epidemic and part of the density . Similar to Fig. 5, the numerical values computed for the CTMC SIS model show good agreement with the estimated values from the simple birth-death approximation.

Fig. 7

The mean and standard deviation, from equations (10), (11), (20), for the time to extinction in the simple birth-death process are graphed as a function of for three initial conditions . Calculated values for the CTMC SIS model (12) are computed from sample paths, where the threshold for reaching an endemic equilibrium is set at (+ for , x for and for ). Parameter values are , , and .

SEIR epidemic model

An SEIR epidemic model includes a latent or an exposed stage, E, where individuals are infected but not yet infectious. The E and I stages in the deterministic model with δ the transition rate from E to I has the form:Linearization of this system near the DFE, , yields , whereand (van den Driessche & Watmough, 2002).

In the corresponding CTMC SEIR model, the transition probabilities include three random variables, , since :The multitype branching process approximation is restricted to E and I. The transition probabilities for the branching process are given as (14) but with and

With no births and deaths, the probability of a minor epidemic is the same as for a CTMC SIR or SIS model, equation (9) and Fig. 3. The parameter δ for transition from exposed to the infectious stage does not affect this probability but does affect the duration of a minor epidemic. Although an analytical formula is not available, the pdf for a minor epidemic can be computed numerically from the solution of the branching process approximation and the backward Kolmogorov differential equations (Appendix E). In Fig. 8, the pdf for the duration a minor or major epidemic are graphed for several values of . For , the duration is much greater as compared to and also, as compared to the CTMC SIR model in Fig. 4. The approximate mean and standard deviation for a minor epidemic are calculated numerically and summarized in Table 5 in Appendix E.

Fig. 8

Probability histograms, computed from sample paths in the CTMC SEIR model show a mixture of two densities when , the probability of a minor or major epidemic. The curve overlayed in red is a graph of the minor epidemic, computed numerically from the branching process approximation and backward Kolmogorov differential equations in Appendix E. Parameter values are , , and with initial conditions , , and .

Conclusion

The birth-death process is well-known in the literature and has been applied to many different biological applications, e.g., (Allen, 2010; Bailey, 1975; Daley & Gani, 1999; Novozhilov et al., 2006; Whittle, 1955). In the birth-death process, either the population size hits zero or approaches infinity. However, conditioned on a finite time to extinction, new analytical formulas are obtained for the mean and higher-order moments for time to extinction when . These analytical formulas can be applied to stochastic SIS, SIR, and SIRS epidemic models, specifically, formulas (10), (11), (20). For large population sizes, the analytical estimates agree with the simulations from the SIR and SIS epidemic models (Fig. 5, Fig. 7). In addition, minor epidemics play a major role in the dynamics when is close to one (Fig. 4, Fig. 6). The mean duration of a minor epidemic increases when is close to one, as predicted by the analytical formula (10). Additional stages prior to the infected stage, such as a latent stage in the SEIR epidemic model prolong the duration of a minor epidemic (Fig. 8). These results have applications to other infectious disease models and to in-host models (Ciupe & Heffernan, 2017).

Emerging diseases, where the value of is close to one, are likely to have minor epidemics with undetected cases of prolonged duration, allowing time for the pathogen to adapt to the host. This can be of special concern in the case of zoonotic infectious diseases, where the infection may be continually re-introduced into human hosts from animal reservoirs or from an intermediate host (Blumberg & Lloyd-Smith, 2013).

Table 3

Estimates for probability of a minor or major epidemic and the mean and standard deviation for duration of a minor epidemic in the CTMC SIR model (hit zero prior to reaching I(t) = 30), , Mean, and STD, are computed from 106 sample paths. These values are compared to the analytical values , μ and σ in the birth-death process. Parameters are N = 5000, γ = 1, and with initial condition I(0) = 3.

ℛ_0	ℙ_0	ℙ_ext	Mean	μ	STD	σ
0.7	1	1	3.49	3.50	3.09	3.12
0.75	1	1	3.80	3.82	3.54	3.58
0.8	1	1	4.19	4.23	4.12	4.23
0.85	1	1	4.73	4.80	4.97	5.21
0.9	1	1	5.45	5.66	6.19	6.92
0.95	1	1	6.55	7.28	8.02	10.93
1.05	0.852	0.864	5.96	7.05	7.29	10.73
1.1	0.764	0.751	5.46	5.34	6.70	6.71
1.15	0.674	0.658	4.75	4.42	5.67	5.00
1.2	0.590	0.579	4.08	3.82	4.58	4.03
1.25	0.521	0.512	3.55	3.39	3.75	3.38
1.3	0.462	0.455	3.16	3.05	3.15	2.93

Table 4

Estimates for probability of a minor epidemic or reaching an endemic level and the mean and standard deviation for duration of a minor epidemic in the CTMC SIS model (hit zero prior to reaching 30), , Mean and STD, are computed from sample paths. These values are compared to the analytical values , μ. and σ in the birth-death process. Parameters are and with initial condition .

N	I¯	ℛ_0	ℙ_0	ℙ_ext	Mean	μ	STD	σ
2000	98	1.05	0.833	0.864	5.50	7.05	6.48	10.73
1000	48	1.05	0.840	0.864	5.60	7.05	6.58	10.73
500	24	1.05	0 0.859	0.864	5.88	7.05	7.10	10.73
2000	182	1.1	0.745	0.751	4.94	5.34	5.63	6.71
1000	91	1.1	0.753	0.751	5.08	5.34	5.87	6.71
500	45	1.1	0.772	0.751	5.36	5.34	6.30	6.71
2000	667	1.5	0.298	0.296	2.23	2.22	1.93	1.92
1000	333	1.5	0.300	0.296	2.24	2.22	1.95	1.92
500	167	1.5	0.305	0.296	2.28	2.22	2.00	1.92

Table 5

The probability of a minor epidemic and the mean μ and standard deviation σ for the duration of a minor epidemic in the SEIR epidemic given initial conditions or for parameter values and .

ℛ_0	ℙ_ext	(E(0),I(0))=(0,3)		(E(0),I(0))=(3,0)
		μ	σ	μ	σ
0.9	1.00	10.48	13.53	11.63	13.47
0.95	1.00	13.82	21.44	14.96	21.38
1.05	0.864	11.82	19.84	12.81	19.80
1.1	0.751	7.84	12.20	8.71	12.20
1.5	0.296	1.325	2.96	1.69	3.25
2	0.125	0.353	1.307	0.517	1.614