Early estimates of epidemic final sizes

ABSTRACT

Early in a disease outbreak, it is important to be able to estimate the final size of the epidemic in order to assess needs for treatment and to be able to compare the effects of different treatment approaches. However, it is common for epidemics, especially of diseases considered dangerous, to grow much more slowly than expected. We suggest that by assuming behavioural changes in the face of an epidemic and heterogeneity of mixing in the population it is possible to obtain reasonable early estimates.

KEYWORDS:

• Epidemics
• final size estimates

AMS 2010 SUBJECT CLASSIFICATION:

• 92D30

1. Introduction

It has been standard practice in analysing disease outbreaks to formulate a dynamical system as a deterministic compartmental model, then to use observed early outbreak data to fit parameters to the model, and finally to analyse the dynamical system to predict the course of the disease outbreak and to compare the effects of different management strategies. In general, such models predict an initial stochastic stage (while the number of infectious individuals is small), followed by a period of exponential growth. Measurement of this early exponential growth rate is an essential step in estimating contact rate parameters for the model.

However, instances have been noted where the growth rate of an epidemic is clearly slower than what this approach predicts. One of the earliest examples [4] concerns the growth of HIV/AIDS in the United States, and a possible explanation might be the mixture of short-term and long-term contacts. This could be a factor in other diseases where there are repeated contacts in family groups and less frequent contacts outside the home. Some suggested explanations for slower than expected growth include metapopulation models with spatial structure including cross-coupling and mobility, clustering in spatial structure, dynamic contacts, agent-based models with differences in infectivity and susceptibility of individuals, and reactive behavioural changes early in a disease outbreak.

In Liberia, a country with a population of close to 4,000,000, the 2014 epidemic of the Ebola virus produced fewer than 12,000 cases although estimates of the reproduction number from early data would have predicted more than 1,5000,000 cases. Such large discrepancies led to growth models that are not mechanistic but are very useful for fitting data and making short-term predictions. One example is a discrete model [5] that has been applied to Ebola data in [6]. Another approach has been to formulate a so-called general growth model of the form $C^{\prime }\left(t\right)=rC(t{)}^p,$ where $C\left(t\right)$ is the number of disease cases occurring up to time t, and p, $0\le p\le 1$, is a ‘deceleration of growth’ parameter [3,12]. These models have been very successful for fitting epidemic growth, but they are not useful for predicting the final size of an epidemic. The parameter p represents the effect of behavioural responses during the initial growth phase of the epidemic but since without some additional term to model the decreasing phase of the epidemic it cannot be incorporated into a final size relation.

An approach that may be more helpful in estimating the final size of an epidemic from early data is the suggestion of [2] that behavioural change may decrease the basic reproduction number in a severe epidemic.

It has been observed that in many disease outbreaks, the phenomenon of superspreading events, situations in which a small fraction of the population causes more than its share of disease cases, has been significant [9]. For example, superspreading events were important in the SARS outbreak of 2002–2003 [10] in which there was a cluster of at least 125 cases apparently infected by a single index patient as well as another cluster of perhaps 300 cases. For a given basic reproduction number, the number of disease cases is generally fewer if there is superspreading than if the population mixing is homogeneous.

In the 2014–2015 West Africa Ebola epidemic, there is also evidence of superspreading [8]. It is suggested in [7] that commonly roughly $20\mathrm{\%}$ of the individuals in a population are responsible for $80\mathrm{\%}$ of the disease cases.

In this note, we suggest that these two effects may give useful early estimates of the final size of an epidemic.

2. Behavioural changes during an epidemic

It is natural to expect that when a disease outbreak begins there may be behavioural response to attempt to protect against infection. This is particularly evident if the disease is considered to be very dangerous, as has been the case in Ebola virus outbreaks, and much less evident with diseases not considered to be dangerous, such as seasonal influenza outbreaks. A model described in [2] for the Ebola virus outbreak of 2014 in West Africa used data from the first 40 weeks of the outbreak to estimate disease model parameters and then assumed a response leading to a decrease by a factor of 0.6 in the reproduction number. This assumption led to a very good fit of the model to continuing disease data and the ultimate final size for the epidemic.

We suggest that assuming a simple model of SIR or SEIR type for an emerging disease outbreak, even though it may fail to capture many of the properties of the disease, together with an assumption of a decrease in the reproduction number after an initial stage, may give a good start on obtaining a reasonable estimate of the epidemic final size. For simplicity, we will assume that the effect of the measures causing this reduction is immediate, even though it would be more realistic to assume a reduction carried out over a given time interval. Clearly, we would expect the reduction factor in the reproduction number and the time of the reduction to be sensitive parameters.

Example 2.1

The 2014 Ebola virus in Liberia has been modelled by a system that includes some specific aspects of the Ebola virus [1,2]. At the beginning of an outbreak of an emerging or re-emerging disease, the specific aspects of the disease may not yet be well understood. In order to make an early estimate of the final size of the outbreak of a disease serious enough that there may be substantial behavioural change in response, we suggest using a generic SEIR model incorporating such response to obtain an initial estimate. We use an SEIR model for the Ebola outbreak in Liberia (1) $\begin{array}{rl}{S}^{\prime }& =-\frac{\beta }{N}SI,\\ E^{\prime }& =\frac{\beta }{N}SI-\eta E,\\ I^{\prime }& =\eta E-\alpha I,\end{array}$(1) with parameter values [2] $\beta =0.139,\alpha =1/10.4,\sigma =0.1,{\mathcal{R}}_0=1.4456,$ and [1] $N=3,787,000,I\left(0\right)=35,$ giving estimates consistent with those of [2].

Numerical simulations give $S\left(280\right)=3,751,655$, indicating 35,345 cases in the first 40 weeks of the outbreak. It is estimated in [2] that after 40 weeks the overall contact rate had dropped to 60 % of its original value. Then, multiplying ${\mathcal{R}}_0$ by a factor 0.6 and using $S\left(280\right)$ as an initial value, we use the final size relation for an SEIR model to estimate $S_{\mathrm{\infty }}$. This final size relation now takes the form $\log \frac{S\left(280\right)}{S_{\mathrm{\infty }}}={\mathcal{R}}_0\left[1-\frac{S_{\mathrm{\infty }}+R\left(280\right)}{N}\right],$ because there are removed members of the population at time t=280. Here, $R\left(280\right)=24,540$. We obtain a final size of $S_{\mathrm{\infty }}=3,694,772$, giving a total number of disease cases over the course of the epidemic as 92,228. This is comparable to but considerably larger than the reported number of 10,678 disease cases in Liberia.

3. Heterogeneity in epidemic models

In order to describe a simple superspreading epidemic model, we consider an SEIR epidemic model in a population divided into two groups. Group 1 members (superspreaders), forming a fraction ρ of the total population, have a contact rate $\sigma a$ and Group 2 members have a contact rate a, with $\sigma >1$. Mixing between groups is proportionate. We take $\sigma =16,\rho =0.2$ in order to conform to the suggestion in [7] that commonly roughly $20\mathrm{\%}$ of the individuals in a population are responsible for $80\mathrm{\%}$ of the disease cases.

We assume

• Total population size is N.

  Group 1 population size is $N_1=\rho N$.

  Group 2 population size is $N_2=(1-\rho )N$.

• Group i is divided into susceptibles ($S_i$), exposed members ($E_i$), infectives ($I_i$) and recovered members ($R_i$).

• Mixing between groups is proportionate. Fraction of contacts of susceptibles with groups 1 and 2 respectively are $p_1=\frac{\rho \sigma }{\rho \sigma +(1-\rho )},p_2=\frac{1-\rho }{\rho \sigma +(1-\rho )}.$

• Infections in each group recover at rate α.

• There are no demographic effects (births, deaths, migration) on the population.

Then $\begin{array}{rl}{S}_1^{\prime }& =-\sigma a{S}_1\left[\frac{\rho \sigma }{\rho \sigma +(1-\rho )}\frac{I_1}{N_1}+\frac{1-\rho }{\rho \sigma +(1-\rho )}\frac{I_2}{N_2}\right]\\ & =-\sigma a{S}_1\left[\frac{\sigma }{\rho \sigma +(1-\rho )}\frac{I_1}{N}+\frac{1}{\rho \sigma +(1-\rho )}\frac{I_2}{N}\right]\\ & =-\sigma a{S}_1\frac{Q}{N}\end{array}$ with $Q=\frac{\sigma I_1+I_2}{\theta },\theta =\rho \sigma +(1-\rho ).$ Similarly, $S_2^{\prime }=-a{S}_2\frac{Q}{N}.$ Thus the full model is (2) $\begin{array}{rl}{S}_1^{\prime }& =-\sigma a{S}_1\frac{Q}{N},\\ S_2^{\prime }& =-a{S}_2\frac{Q}{N},\\ E_1^{\prime }& =\sigma a{S}_1\frac{Q}{N}-\eta E_1,\\ E_2^{\prime }& =a{S}_2\frac{Q}{N}-\eta E_2,\\ I_1^{\prime }& =\eta E_1-\alpha I_1,\\ I_2^{\prime }& =\eta E_2-\alpha I_2.\end{array}$(2) The parameter values corresponding to the SEIR model (1(1) $\begin{array}{rl}{S}^{\prime }& =-\frac{\beta }{N}SI,\\ E^{\prime }& =\frac{\beta }{N}SI-\eta E,\\ I^{\prime }& =\eta E-\alpha I,\end{array}$(1) ) are (3) $\begin{array}{rl}\sigma & =16,\rho =0.2,a=1/10.4,\alpha =0.0107,\eta =0.1,N=3,787,000\\ S_1\left(0\right)& =757,390,E_1\left(0\right)=2,I_1\left(0\right)=8\\ S_2\left(0\right)& =3,029,570,E_2\left(0\right)=4,I_2\left(0\right)=32.\end{array}$(3) To calculate the reproduction number using the next generation method [11], we write $F=\left[\begin{array}{cccc}0& \frac{\rho {\sigma }^{2}a}{\theta }& 0& \frac{\sigma \rho }{\theta }\\ 0& 0& 0& 0\\ 0& \frac{(1-\rho )\sigma a}{\theta }& 0& \frac{(1-\rho )a}{\theta }\\ 0& 0& 0& 0\end{array}\right],V=\left[\begin{array}{cccc}\sigma & 0& 0& 0\\ -\sigma & \alpha & 0& 0\\ 0& 0& \sigma & 0\\ 0& 0& -\sigma & \alpha \end{array}\right].$ Thus the next generation matrix is $F{V}^{-1}=\left[\begin{array}{cccc}\frac{\rho {\sigma }^{2}a}{\alpha \theta }& \frac{\rho {\sigma }^{2}a}{\alpha \theta }& \frac{\rho \sigma a}{\alpha \theta }& \frac{\rho \sigma a}{\alpha \theta }\\ 0& 0& 0& 0\\ \frac{(1-\rho )\sigma a}{\alpha \theta }& \frac{(1-\rho )\sigma a}{\alpha \theta }& \frac{a(1-\rho )}{\alpha \theta }& \frac{a(1-\rho )}{\alpha \theta }\\ 0& 0& 0& 0\end{array}\right].$ This matrix has the same non-zero eigenvalues as the $2\times 2$ matrix $K_L=\left[\begin{array}{cc}\frac{\rho {\sigma }^{2}a}{\alpha \theta }& \frac{\rho \sigma a}{\alpha \theta }\\ \frac{a(1-\rho )}{\alpha \theta }& \frac{a(1-\rho )}{\alpha \theta }\end{array}\right]$ (the next generation matrix with least domain). The matrix $K_L$ has determinant zero. Thus the non-zero eigenvalue is the trace of this matrix, and ${\mathcal{R}}_0=\frac{a}{\alpha \theta }[{\sigma }^2\rho +(1-\rho \left)\right]=\frac{a}{\alpha }\frac{{\sigma }^2\rho +(1-\rho )}{\sigma \rho +(1-\rho )}.$ From the model (2(2) $\begin{array}{rl}{S}_1^{\prime }& =-\sigma a{S}_1\frac{Q}{N},\\ S_2^{\prime }& =-a{S}_2\frac{Q}{N},\\ E_1^{\prime }& =\sigma a{S}_1\frac{Q}{N}-\eta E_1,\\ E_2^{\prime }& =a{S}_2\frac{Q}{N}-\eta E_2,\\ I_1^{\prime }& =\eta E_1-\alpha I_1,\\ I_2^{\prime }& =\eta E_2-\alpha I_2.\end{array}$(2) ), we see that $(S_1+E_1{)}^{\prime }=-\eta E_1,(S_1+E_1+I_1{)}^{\prime }=-\alpha I_1,$ from which we deduce that $E_1\to 0,I_1\to 0$, and $S_1\left(t\right)\to S_1\left(\mathrm{\infty }\right)\ge 0$. Integration gives $\alpha {\int }_0^{\mathrm{\infty }}{I}_1\left(t\right)\mathrm{d}t=N_1-S_1\left(\mathrm{\infty }\right).$ Similarly, we see that $(S_2+E_2{)}^{\prime }=-\eta E_2,(S_2+E_2+I_2{)}^{\prime }=-\alpha I_2,$ from which we deduce that $E_2\to 0,I_2\to 0$, and $S_2\left(t\right)\to S_2\left(\mathrm{\infty }\right)\ge 0$.

Integration gives $\alpha {\int }_0^{\mathrm{\infty }}{I}_1\left(t\right)\mathrm{d}t=N_1-S_1\left(\mathrm{\infty }\right)$ and $\alpha {\int }_0^{\mathrm{\infty }}{I}_2\left(t\right)\mathrm{d}t=N_2-S_2\left(\mathrm{\infty }\right).$ Integration of the equation for $S_1$ in the model gives $\log \frac{S_1\left(0\right)}{S_1(\mathrm{\infty }}=\frac{\sigma a}{N}{\int }_0^{\mathrm{\infty }}Q\left(t\right)\mathrm{d}t.$ We combine the expressions $Q=\frac{\sigma I_1+I_2}{\theta }$ and $\alpha {\int }_0^{\mathrm{\infty }}\sigma I_1\left(t\right)+I_2\left(t\right)\mathrm{d}t=\sigma [N_1-S_1(\mathrm{\infty }\left)\right]+[N_2-S_2(\mathrm{\infty }\left)\right]$ to give $\log \frac{S_1\left(0\right)}{S_1\left(\mathrm{\infty }\right)}=K_{11}\left[1-\frac{S_1\left(\mathrm{\infty }\right)}{N_1}\right]+K_{12}\left[1-\frac{S_2\left(\mathrm{\infty }\right)}{N_2}\right].$ Since $S_1^{\prime }/S_1=\sigma S_2^{\prime }/S_2$, we may carry out analogous calculations to give $\log \frac{S_2\left(0\right)}{S_2\left(\mathrm{\infty }\right)}=K_{21}\left[1-\frac{S_1\left(\mathrm{\infty }\right)}{N_1}\right]+K_{22}\left[1-\frac{S_2\left(\mathrm{\infty }\right)}{N_2}\right],$ with the matrix K given by $K=\left[\begin{array}{cc}{K}_{11}& K_{12}\\ K_{21}& K_{22}\end{array}\right]=\left[\begin{array}{cc}\frac{\rho {\sigma }^{2}a}{\alpha \theta }& \frac{(1-\rho )\sigma a}{\alpha \theta }\\ \frac{\rho \sigma a}{\theta \alpha }& \frac{(1-\rho )a}{\alpha \theta }\end{array}\right].$ We note that K is the transpose of the next generation matrix with least domain $K_L$. We now have the final size relation (4) $\begin{array}{rl}\log \frac{S_1\left(0\right)}{S_1\left(\mathrm{\infty }\right)}& =K_{11}\left[1-\frac{S_1\left(\mathrm{\infty }\right)}{N_1}\right]+K_{12}\left[1-\frac{S_2\left(\mathrm{\infty }\right)}{N_2}\right]\\ \log \frac{S_2\left(0\right)}{S_2\left(\mathrm{\infty }\right)}& =K_{21}\left[1-\frac{S_1\left(\mathrm{\infty }\right)}{N_1}\right]+K_{22}\left[1-\frac{S_2\left(\mathrm{\infty }\right)}{N_2}\right].\end{array}$(4) We note that ${\mathcal{R}}_0=K_{11}+K_{22}.$ Because $S_1^{\prime }/S_1=\sigma S_2^{\prime }/S_2$, we have $\log \frac{S_1\left(0\right)}{S_1\left(\mathrm{\infty }\right)}=\sigma \log \frac{S_2\left(0\right)}{S_2\left(\mathrm{\infty }\right)},$ which implies (5) $S_1\left(\mathrm{\infty }\right)=\frac{S_1\left(0\right)}{S_2(0{)}^{\sigma }}{S}_2(\mathrm{\infty }{)}^{\sigma }.$(5) We could substitute (5(5) $S_1\left(\mathrm{\infty }\right)=\frac{S_1\left(0\right)}{S_2(0{)}^{\sigma }}{S}_2(\mathrm{\infty }{)}^{\sigma }.$(5) ) into the final size relation (4(4) $\begin{array}{rl}\log \frac{S_1\left(0\right)}{S_1\left(\mathrm{\infty }\right)}& =K_{11}\left[1-\frac{S_1\left(\mathrm{\infty }\right)}{N_1}\right]+K_{12}\left[1-\frac{S_2\left(\mathrm{\infty }\right)}{N_2}\right]\\ \log \frac{S_2\left(0\right)}{S_2\left(\mathrm{\infty }\right)}& =K_{21}\left[1-\frac{S_1\left(\mathrm{\infty }\right)}{N_1}\right]+K_{22}\left[1-\frac{S_2\left(\mathrm{\infty }\right)}{N_2}\right].\end{array}$(4) ) to obtain a single final size equation.

We simulate the superspreader model (2(2) $\begin{array}{rl}{S}_1^{\prime }& =-\sigma a{S}_1\frac{Q}{N},\\ S_2^{\prime }& =-a{S}_2\frac{Q}{N},\\ E_1^{\prime }& =\sigma a{S}_1\frac{Q}{N}-\eta E_1,\\ E_2^{\prime }& =a{S}_2\frac{Q}{N}-\eta E_2,\\ I_1^{\prime }& =\eta E_1-\alpha I_1,\\ I_2^{\prime }& =\eta E_2-\alpha I_2.\end{array}$(2) ) for $0\le t\le 280$, using the parameter values and initial conditions (3(3) $\begin{array}{rl}\sigma & =16,\rho =0.2,a=1/10.4,\alpha =0.0107,\eta =0.1,N=3,787,000\\ S_1\left(0\right)& =757,390,E_1\left(0\right)=2,I_1\left(0\right)=8\\ S_2\left(0\right)& =3,029,570,E_2\left(0\right)=4,I_2\left(0\right)=32.\end{array}$(3) ), obtaining $S_1\left(280\right)=746,830,R_1\left(280\right)=7350,S_2\left(280\right)=3,026,915,R_2\left(280\right)=1870.$ The number of disease cases for $0\le t\le 280$ is $\left[S_1\right(0)-S_1(280\left)\right]+\left[S_2\right(0)-S_2(280\left)\right]=10,560+2655=13,215.$ Next, we use the final size relations (4(4) $\begin{array}{rl}\log \frac{S_1\left(0\right)}{S_1\left(\mathrm{\infty }\right)}& =K_{11}\left[1-\frac{S_1\left(\mathrm{\infty }\right)}{N_1}\right]+K_{12}\left[1-\frac{S_2\left(\mathrm{\infty }\right)}{N_2}\right]\\ \log \frac{S_2\left(0\right)}{S_2\left(\mathrm{\infty }\right)}& =K_{21}\left[1-\frac{S_1\left(\mathrm{\infty }\right)}{N_1}\right]+K_{22}\left[1-\frac{S_2\left(\mathrm{\infty }\right)}{N_2}\right].\end{array}$(4) ) using the compartment sizes for t=280 as initial values. This yields $S_1\left(\mathrm{\infty }\right)=729,142,S_2\left(\mathrm{\infty }\right)=3,022,384.$ The number of new disease cases after t=280 is 17,718+4531=22,249 , and the total number of cases is 13,215+22,249=35,464. This is still larger than the actual number of observed cases, but considerably closer to the actual number of observed cases than the estimate provided by (1(1) $\begin{array}{rl}{S}^{\prime }& =-\frac{\beta }{N}SI,\\ E^{\prime }& =\frac{\beta }{N}SI-\eta E,\\ I^{\prime }& =\eta E-\alpha I,\end{array}$(1) ).

The factor σ in the model (2(2) $\begin{array}{rl}{S}_1^{\prime }& =-\sigma a{S}_1\frac{Q}{N},\\ S_2^{\prime }& =-a{S}_2\frac{Q}{N},\\ E_1^{\prime }& =\sigma a{S}_1\frac{Q}{N}-\eta E_1,\\ E_2^{\prime }& =a{S}_2\frac{Q}{N}-\eta E_2,\\ I_1^{\prime }& =\eta E_1-\alpha I_1,\\ I_2^{\prime }& =\eta E_2-\alpha I_2.\end{array}$(2) ) represents the degree of superspreading. To indicate the extent to which the epidemic final size would depend on the value of the parameter σ, we used the final size relation (4(4) $\begin{array}{rl}\log \frac{S_1\left(0\right)}{S_1\left(\mathrm{\infty }\right)}& =K_{11}\left[1-\frac{S_1\left(\mathrm{\infty }\right)}{N_1}\right]+K_{12}\left[1-\frac{S_2\left(\mathrm{\infty }\right)}{N_2}\right]\\ \log \frac{S_2\left(0\right)}{S_2\left(\mathrm{\infty }\right)}& =K_{21}\left[1-\frac{S_1\left(\mathrm{\infty }\right)}{N_1}\right]+K_{22}\left[1-\frac{S_2\left(\mathrm{\infty }\right)}{N_2}\right].\end{array}$(4) ) for the model (without incorporating behavioural response).

With $\sigma =16$, so that 20% of the population is responsible for 80% of the disease cases, we obtain $S_1\left(\mathrm{\infty }\right)=338,577,S_2\left(\mathrm{\infty }\right)=2,880,893,$ so that there are 567,490 disease cases overall.

With $\sigma =8$, so that 20% of the population is responsible for two-thirds of the disease cases, we obtain $S_1\left(\mathrm{\infty }\right)=323,658,S_2\left(\mathrm{\infty }\right)=2,724,124,$ so that there are 739,118 disease cases overall.

With $\sigma =4$, so that 20% of the population is responsible for half of the disease cases, we obtain $S_1\left(\mathrm{\infty }\right)=286,593,S_2\left(\mathrm{\infty }\right)=2,376,116,$ so that there are 1,124,251 disease cases overall.

This suggests that superspreading is significant in decreasing the number of disease cases corresponding to a given value of the basic reproduction number.

4. Conclusion

It is useful to have some idea early in a disease outbreak how serious the outbreak will turn out to be. If the disease is viewed as a serious threat, there will surely be substantial behavioural changes as the disease develops. The standard approaches to estimating a disease outbreak using initial exponential growth rates may give predictions which are so large as to be meaningless. We suggest that it may be possible to use relatively simple SIR or SEIR models without detailed assumptions about the particular disease being studied but assuming a decrease in the contact rate and heterogeneity of mixing to obtain more plausible estimates. While a decrease in the contact rate would probably be continuous in reality, we approximate it by a discontinuous decrease at a fixed time. Of course, the results will depend strongly on the size of the assumed decrease and the time at which it is assumed to take place.

Disclosure statement

No potential conflict of interest was reported by the author.

Additional information

Funding

This research was supported by the Natural Sciences and Engineering Council of Canada [grant no. RGPIN-2016-03706].