Estimating the transmission potential of supercritical processes based on the final size distribution of minor outbreaks.

Use of the final size distribution of minor outbreaks for the estimation of the reproduction numbers of supercritical epidemic processes has yet to be considered. We used a branching process model to derive the final size distribution of minor outbreaks, assuming a reproduction number above unity, and applying the method to final size data for pneumonic plague. Pneumonic plague is a rare disease with only one documented major epidemic in a spatially limited setting. Because the final size distribution of a minor outbreak needs to be normalized by the probability of extinction, we assume that the dispersion parameter (k) of the negative-binomial offspring distribution is known, and examine the sensitivity of the reproduction number to variation in dispersion. Assuming a geometric offspring distribution with k=1, the reproduction number was estimated at 1.16 (95% confidence interval: 0.97-1.38). When less dispersed with k=2, the maximum likelihood estimate of the reproduction number was 1.14. These estimates agreed with those published from transmission network analysis, indicating that the human-to-human transmission potential of the pneumonic plague is not very high. Given only minor outbreaks, transmission potential is not sufficiently assessed by directly counting the number of offspring. Since the absence of a major epidemic does not guarantee a subcritical process, the proposed method allows us to conservatively regard epidemic data from minor outbreaks as supercritical, and yield estimates of threshold values above unity. Crown

Introduction

The basic reproduction number, R 0 is a summary measure of the transmission potential of infectious diseases, interpreted as the average number of secondary cases generated by a typical primary case throughout its entire course of infection in a fully susceptible population (Diekmann and Heesterbeek, 2000). Because of the well-known threshold property (Anderson, 1991), given R 0<1 we observe only minor outbreaks in which the infectious agent goes extinct before affecting a substantial fraction of the population (Diekmann and Heesterbeek, 2000). Given R 0>1, an introduction of a single infected individual into a population could potentially lead to a major epidemic, which does not naturally decline until the pool of susceptible individuals falls below a critical proportion, and in which substantial fraction of the population thus experiences infection (though, as we will discuss in the present study, it must be remembered that R 0>1 can also result in minor outbreaks). Accordingly, R 0 has played a key role in identifying optimal strategies for infectious disease control (Keeling and Rohani, 2007), including the determination of vaccination coverage required to eradicate an infectious disease (Anderson and May, 1991).

Various statistical methods have been employed to estimate R 0 based on empirical data from epidemic or endemic diseases. These include maximum likelihood method and the use of estimators obtained from biomathematical and theoretical studies (Dietz, 1993, Farrington et al., 2001, Li et al., 2011). Recently, the estimation of R 0 has mainly focused on the analysis of early exponential (linear) growth rates (Chowell et al., 2007, Vynnycky et al., 2007). This type of analysis requires the prior observation of an actual major epidemic. However, it is frequently the case for rare and exotic diseases that we rarely (or never) observe a major epidemic prior to an actual large-scale epidemic event (e.g. severe acute respiratory syndrome (SARS) from 2002 to 2003).

Also, an observed single epidemic curve represents only a single sample path profile (i.e. a single stochastic realization) of all possible sample paths (Lindsey, 2001, van den Broek and Nishiura, 2009), and thus, a model and estimate for R 0 explicitly accounting for demographic stochasticity are needed (Becker, 1989, Andersson and Britton, 2000). Branching process approximation is one stochastic modeling approach. The discrete-type branching process is referred to as Galton–Watson process (Harris, 1963), and the process with infection–age structure is referred to as Crump–Mode–Jagers (CMJ) process (Jagers, 1975, Mode and Sleeman, 2000). Other than the branching processes, a martingale model has been employed to estimate R 0 and its uncertainty based on a single sample of the final size of a major epidemic (i.e. based on a single observation of cumulative incidence) (Becker, 1989). For theoretical purposes, not only the analysis of single samples of the final size, but also final size distributions of major epidemics have been well analyzed through the so-called Sellke-construction (Andersson and Britton, 2000, Ball and O'Neill, 1999, Ball and Clancy, 1993, Nishiura et al., 2011a). While major epidemics have attracted modelers’ interest, little attention has been paid to the final size distribution of minor outbreaks. Only a small number of studies have analyzed the final size distribution of minor outbreaks adopting the key assumption that a subcritical process (i.e. R 0<1) resulted in the observation (Becker, 1974). Classical mortal branching process models have been applied to surveillance data for diseases such as measles and avian influenza (Farrington et al., 2003, Ferguson et al., 2004). However, the final size distribution of minor outbreak data given a supercritical process (i.e. R 0>1) have not been explicitly examined to date.

The purpose of the present study is to estimate R 0 for a supercritical process based on the final size distribution of minor outbreaks. We illustrate the estimation method with application to primary pneumonic plague. In the next section, we begin with a description of our practical motivations.

Materials and methods

Pneumonic plague

This study considers the human-to-human transmission of primary pneumonic plague, a fatal clinical presentation of bacterial infection with Yersinia pestis. Presently, primary pneumonic plague is rarely observed; though the form of bubonic plague transmitted by flea bites (in humans and rodents) is sometimes reported (Pollitzer, 1954). When infected fleas bite humans, transmission of the bubonic plague with fever and lymph node swelling is possible. If Y. pestis gets into the blood stream, septicemic plague occurs. Systemic circulation of Y. pestis can also progress to the “secondary pneumonic plague”, a lung infection spreadable through respiratory droplet and from which an epidemic of primary pneumonic plague can occur. However, pneumonic plague is not transmitted by aerosols (i.e. not by droplet nuclei) but by droplet infection (Kool, 2005). In the absence of early treatment, the case fatality ratio of pneumonic plague can reach 100% (Pollitzer, 1954). Due to the high virulence of Y. pestis, industrialized countries across the world have identified its potential for use in biological warfare (Inglesby et al., 2000). For these reasons, knowledge of the human-to-human transmission potential of pneumonic plague is a critical question among experts.

Several studies have analyzed historical outbreak data of primary pneumonic plague, estimating the reproduction number R of human-to-human transmission (i.e. the average number of secondary human cases directly caused by a single human primary case of pneumonic plague). Gani and Leach (2004) fitted a geometric distribution directly to observed distributions of secondary cases produced by a single primary case based on eight minor outbreak data sets showing transmission networks (i.e. a network showing who acquired infection from whom based on contact tracing). They found a mean and variance of 1.3 and 3.1, respectively. However, minor outbreak data are limited to outbreaks that decline to extinction and include substantial numbers of primary cases that did not generate any secondary transmissions. To avoid potential underestimation of R , Nishiura and Kakehashi (2005) focused on time-dependency in the transmissibility and Nishiura et al. (2006) estimated the initial value of the reproduction number (i.e. theoretical R at time 0 of an epidemic) as ranging from 2.8 to 3.5. Lloyd-Smith et al. (2005) took into account the effectiveness of public health interventions, fitting a branching process model with a negative binomially distributed offspring distribution to a transmission network, estimating R as 1.32 (90% confidence interval (CI): 1.01–1.61) and dispersion parameter as k=1.37 (90% CI: 0.88–3.53). They identified that published outbreak data tended to start with a superspreading event (implying that an R of 2.8–3.5 in Nishiura et al. (2006) may be overestimated), and that the offspring distribution of pneumonic plague is not overdispersed. Other published studies on this subject include an analysis of minor outbreak data by directly counting the number of offspring (Egan, 2010) and a reanalysis of major epidemic data from Manchuria, 1911 in which the absence of bubo among pneumonic plague cases was partly noted (Nishiura, 2006). Overall, the mean number of human secondary cases per single human primary case (i.e. R ) is likely to be greater than unity (because we observed an actual major epidemic in Manchuria) with a relatively small variance in the offspring distribution. Since historical data for pneumonic plague are scarce, the finding of R being slightly greater than unity should be validated using several other methods.

The present study examined historical minor outbreak data (see Table 1 for a summary). We include outbreaks known to be the pneumonic plague (through bacteriological confirmation) and involving two or more cases (i.e. outbreaks with at least one secondary transmission event). These limits are necessary as the frequency of bubonic plague is much higher than that of pneumonic plague in literature, and it was not considered feasible to collect reports with only one pneumonic plague case. Indeed, the reporting of no offspring is usually underestimated (Lloyd-Smith, 2007). We also ignored transmission events from animals to humans (e.g. dog-to-human transmission). Since an explicit estimation of the effectiveness of interventions is too challenging from only small minor outbreak data, and because previous studies documented extinction events or a decline in the number of cases before implementing control measures (Egan, 2010, Nishiura et al., 2006), we ignore the influence of public health interventions on the transmission dynamics. It should be noted that a few published studies explicitly took into account the timing and effectiveness of interventions (Gani and Leach, 2004, Lloyd-Smith et al., 2005, Egan, 2010). As bacteriological confirmation was a prerequisite for inclusion, most of the outbreaks studied occurred in the 20th century. Included are 18 minor outbreaks and 1 major epidemic. Among the 18 minor outbreaks, there was mean of 14.7 cases (standard deviation: 13.3), and a median of 11 cases (quartiles: 3–21). All of these outbreaks started with only one index case.

Table 1

Reported outbreaks of primary pneumonic plague since the 20th century.

Year	Location	Total number of casesa	Original references	Secondary use of the data in modeling analyses
1906/07	Shotley, UK	8	Bulstrode (1927), Van Zwanenberg (1970), Black and Black (2000)	Egan (2010)
1907	Seattle, USA	5	Anonymous (1955)	Gani and Leach (2004)
1910	Freston, UK	4	Bulstrode (1927)	Egan (2010)
1910/11	Manchuria, China	5009	Temporary Quarantine Section, Kanto Totokufu (1912)	Nishiura (2006)
1918	Erwarton, UK	2	Van Zwanenberg (1970), Black and Black (2000)	Egan (2010)
1919	Oakland, USA	13	Kellogg (1920)	Gani and Leach (2004), Nishiura and Kakehashi (2005)
1924	Los Angeles, USA	35	Dickie (1926), Viseltear (1974)
1939	Ecuador	18	Murdock (1940)	Gani and Leach (2004), Nishiura and Kakehashi (2005)
1946	Mukden, China	39	Tieh et al. (1948)	Gani and Leach (2004), Nishiura et al. (2006)
1946	Rangoon	16	Wynne-Griffith (1948)	Gani and Leach (2004), Nishiura and Kakehashi (2005)
1957	Madagascar	42	Brygoo and Gonon (1958)	Gani and Leach (2004), Nishiura et al. (2006)
1993	Zambia	3	McLean (1995)	Gani and Leach (2004)
1997	Madagascar	18	Ratsitorahina et al. (2000)	Gani and Leach (2004)
1998	Ecuador	12	Gabastou et al. (2000)
2002	India	30	Joshi et al. (2009)
2002	India	10	Gupta and Sharma (2007)
2004	Uganda (#1)b	2	Begier et al. (2006)
2004	Uganda (#2)b	2	Begier et al. (2006)
2006	Uganda	6	Anonymous (2009)

The total number of cases includes index case.

There were two independent clusters of pneumonic plague cases.

Not only to validate published findings of transmission potential, but we used the minor outbreak data to overcome the scarcity of data on pneumonic plague. The occurrence of a major epidemic in Manchuria ( Fig. 1) possibly suggests that R >1. Using the pneumonic plague as an example, we aim to demonstrate a method that uses minor outbreak data for an assessment of the reproduction number in a supercritical process. Where possible, reproduction numbers should be estimated using major epidemic data. However, bites from infected fleas also cause human plague. Therefore, major epidemic data that involved many cases without confirmation of the absence of bubo and are not independent of flea biting may not be appropriate for the estimation of R , because the time evolution of cases may be influenced by other modes of transmission (e.g. those arising directly or indirectly from bubonic plague cases). Also, dead rodents were observed during the 1911 epidemic in Manchuria, indicating an involvement of bubonic plague (Temporary Quarantine Section, Kanto Totokufu, 1912). To limit our interest to human-to-human transmission, we do not use R 0 to refer to the transmission potential of pneumonic plague and hereafter consistently use R .

Fig. 1

The distribution of final epidemic size. The final size distribution of a homogeneously mixing population usually exhibits a bimodal shape where the first peak represents minor outbreaks and the second peak corresponds to the mode of major epidemics. Arrows indicate the available evidence for the final size of pneumonic plague (18 minor outbreaks and 1 major epidemic).

Probability generating functions of supercritical branching process

Assuming R >1, we considered a statistical model for analyzing the final size distribution of minor outbreaks. We used a discrete-time Galton–Watson (GW) process. This has been used to estimate a subcritical threshold parameter based on a mortal branching process model (Becker, 1974, Farrington et al., 2003). The model can also be regarded as an embedded age-dependent CMJ process (Jagers, 1975, Mode and Sleeman, 2000). Let Y be the integer-valued random variable representing the total size of an epidemic (including the index case), described by the embedded GW process with X cases in generation n. We have

Let p represent the probability that the total number of cases is y, i.e., Pr(Y=y). Also, let π be the probability of extinction. Since R is assumed as greater than 1, the final size Y is infinite with probability (1−π). In other words, we have

Let the probability mass function of the offspring distribution ξ be f . As a well-known property of the GW process (Harris, 1963), the probability generating function (pgf) h(s) is written asfor s∈[0,1]. With an assumption R >1, we assume that R =h′(1)>1. If the offspring distribution follows a negative binomial distribution with dispersion parameter k, the probability of extinction is calculated by

With a gamma-distributed infectious period and constant infectiousness, k in Eq. (4) could be interpreted as being equivalent to an inverse of the coefficient of variation of the infectious period (Anderson and Watson, 1980, Britton and Lindenstrand, 2009). Poisson and geometric distributions are special cases of the negative binomial distribution with k→∞ and k=1, respectively. The derivatives of h(s), which are required later, areandSimilarly, the total outbreak size is also considered, i.e.Here we consider g(s) as if it were the probability generating function (though it is not, because g(1)=π≤1 (Yan, 2008)). Eq. (7) is useful to determine the probability p . For n=1, 2, 3, …, the sequence of partial sums is considered:

For s∈[0.1], letbe the pgf of the random variable Y . Since Y 1=1+X 1, we have (Mode and Sleeman, 2000):Moreover, let {Y ; k=1, 2, …, X 1} be a collection of independently and identically distributed copies of Y . For every n≥1,leading us to find thatBy applying the dominated convergence theorem, the pgf of the random variable Y satisfiesfor s∈[0,1].

The final size distribution of minor outbreaks

Since R >1, we have Eq. (2), and thus, the conditional distribution of the final size which starts with one primary case (conditioned on minor outbreak) is

Closed-form solution does not always exist for g(s). For the negative binomial distribution, Eq. (13) is rewritten as

For y=1,

For y≥2, the distribution is recursively calculated as (Yan, 2008)

These can be numerically checked using (5), (6) for y=2 and 3. When k→∞, Eq. (17) agrees with the Borel–Tanner distribution (Farrington et al., 2003). When k=1, the distribution is geometric:as derived by Mode and Sleeman (2000) and confirmed by Yan (2008).

Given a total of m complete observations, each with information of the final size n , including an index case, and assuming that k is known, the likelihood function to estimate R isHowever, we excluded all observations of y=1 (i.e. index cases with no secondary case were excluded) and the observed distribution is thus truncated. The likelihood should now read

The probability of extinction π depends on both R and k, but assuming that k is known for now, π is written as a function of R . We consider five different values of k being 1/3, 1/2, 1, 2 and 3, because a previous study analyzing transmission networks of pneumonic plague estimated k at 1.37 (90% CI: 0.88–3.53) (Lloyd-Smith et al., 2005). Five possible k correspond to a plausible range of the dispersion parameter, and also account for the fact that overdispersed k<1 are more likely to be missed in empirical observations. Provided that k=1/2, 1 or 2, we havefor k=1/2,as is well known for k=1 (i.e. the case of geometric distribution), andfor k=2. To compare likelihood values to assess goodness of fit, we used Akaike Information Criterion (AIC). In the case of k=1, it is a conventional geometric offspring distribution with only one unknown parameter (R ). In the other four scenarios (k≠1), we need the parameters R and k, and both were regarded as unknown.

Estimation of the reproduction number using two other methods

Focusing on the data in Table 1 we used two other methods to estimate R . The first was to estimate R based on the observed probability of extinction, assuming that the dispersion parameter of the offspring distribution is known. Among 19 samples in Table 1, only one major epidemic was reported, and thus, the probability of extinction was 18/19. Using this, we estimated R from the relationship π=18/19. It should be emphasized that this method required us to adopt an unrealistic assumption that the data in Table 1 fully represented the total occurrences of pneumonic plague outbreaks (including immediate extinctions and thus the final size distribution is not truncated).

Second, we estimated R from the major epidemic data. Fig. 2 shows an epidemic curve in Manchuria, 1911. We estimated the intrinsic growth rate r, which is also referred to as the Malthusian growth rate. Assuming that generation time follows a gamma distribution with mean μ=5.66 days and coefficient of variation v=0.644 (where the standard deviation σ satisfies σ=μv=3.65 days), the reproduction number was calculated using the following estimator derived from the Euler–Lotka equation (Wallinga and Lipsitch, 2007):

Fig. 2

Epidemiological dynamics of primary pneumonic plague in Manchuria, China from 1910 to 1911: (A) daily number of plague notifications in and around the area under direct political control of the Japanese empire (n=5009). (B) Probability density function of the serial interval for primary pneumonic plague. The mean (and standard deviation) serial interval has been estimated at 5.66 (3.65) days (Nishiura, 2006).

This argument permits us to replace r by the right-hand side of

We estimated R (or r) based on a pure birth process. Given observations of the cumulative number of cases on day 0, 1, …, t, C(0), C(1), …C(t), we have

which can be used to construct a likelihood function for r:

Detailed derivation of the likelihood is given elsewhere (Bailey, 1964, Nishiura et al., 2009). The likelihood function (27) was used for estimating R . Profile likelihood was used to derive 95% confidence intervals. Since the generation time of pneumonic plague has yet to be fully quantified (Nishiura, 2010, Klinkenberg and Nishiura, 2011), we investigated the sensitivity of R to different μ ranging from 3.66 to 7.66 days, corresponding to empirically observed quartiles of 3–7 days (Nishiura, 2006).

Results

Estimate of R based on the sample probability of extinction

Assuming that the samples in Table 1 (n=19) represented the total presentations of a pneumonic plague case into a fully susceptible population, the probability of epidemic would be calculated as 100(1−π)=5.26%. If the number of secondary cases per single primary case is geometrically distributed, we get

For other dispersion parameters of the negative-binomially distributed offspring distribution (k=1/3, 1/2, 2 and 3), R is calculated using a solution of (4) (e.g. see (21), (22), (23)) and estimated as 1.12, 1.08, 1.04 and 1.04, respectively. As the offspring distribution becomes overdispersed, the estimate of R increases due to the fixed nature of the data used.

Estimate of R based on the final size distribution of minor outbreaks

Table 2 shows the estimates of R based on the final size of minor outbreaks. As mentioned above, we assumed k as a known quantity and considered possible integer values of k (or an inverse of integer values) in realistic ranges. When the offspring distribution is geometric, R is 1.16 (95% CI: 0.97, 1.38). When overdispersed with k=1/3 and 1/2, the maximum likelihood estimates of R were 1.22 and 1.19, respectively. When less dispersed with k=2 and 3, R was estimated at 1.14 and 1.13, respectively. This qualitative pattern occurs because empirical data are fixed and only the dispersion parameter is varied. In terms of AIC, the geometric distribution appeared to yield the best fit, consistent with Lloyd-Smith et al. (2005). The estimate based on the sample probability of extinction was not significantly deviated from the above-mentioned crude estimate.

Table 2

Estimates of the threshold parameter of primary pneumonic plague based on minor outbreak data (n=18).

Dispersion parameter (k)		Reproduction number (95% CI)a	AICb
1/3	Overdispersed	1.22 (0.95, 1.60)	139.0
1/2	↑	1.19 (0.96, 1.49)	138.8
1	Geometric distribution	1.16 (0.97, 1.38)	136.3
2	↓	1.14 (0.98, 1.32)	137.7
3	Less dispersed	1.13 (0.98, 1.30)	137.5

CI, Confidence intervals, derived from profile likelihood.

AIC, Akaike Information Criterion.

Estimate of R based on growth data from a major epidemic

Fig. 3A shows the results from conditionally fitting exponential growth to the major epidemic data. Assuming that the mean generation time is 5.66 days, the reproduction number is estimated as 2.16 (95% CI: 2.06–2.28) which is significantly greater than the estimates of R based on minor outbreak data. As the mean generation time varied from 3.66 to 7.66 days, the maximum likelihood estimate of R also varied from 1.59 to 2.96 (Fig. 3B). The lower 95% CI of R with the shortest assumed mean generation time (3.66 days) was 1.54 and was significantly greater than estimates based on minor outbreak data, except for the unrealistically overdispersed offspring distribution (k≤1/3).

Fig. 3

Transmission dynamics of primary pneumonic plague in Manchuria, China from 1910 to 1911: (A) observed (white) and predicted (gray) daily numbers of primary pneumonic plague cases. Assuming exponential growth, the daily growth rate (and the 95% confidence intervals) was estimated at 0.16 (0.15, 0.17) per day. Day 1 corresponds to 1 January 1911. (B) The sensitivity of the reproduction number to different mean lengths of the generation time ranging from 3.65 to 7.65 day. Circles represent maximum likelihood estimates, while the whiskers extend to the upper and lower 95% confidence intervals based on the profile likelihood.

Discussion

The present study proposed a method for the estimation of the threshold parameter of a supercritical process, based on the final size distribution of minor outbreaks. The method was applied to pneumonic plague, which rarely occurs and has limited availability of major epidemic data. To estimate R using the proposed method, the likelihood function must be conditioned on minor outbreak and include the probability of extinction in the denominator. Therefore, the dispersion parameter was assumed as known to rewrite the probability of extinction as a function of R (and eliminate the probability of extinction from the likelihood). For this reason, we examined the sensitivity of R to various dispersion parameters k of the negative binomial distribution. As a result, the geometric offspring distribution yielded R =1.16 and a less dispersed offspring distribution with k=2 suggested an R of 1.14. These estimates are in concordance with those reported from an analysis of transmission networks (Lloyd-Smith et al., 2005) and did not differ substantially from estimates calculated from directly counting the number of secondary cases per primary case from a transmission tree (Gani and Leach, 2004). Overall, the human-to-human transmission potential of plague is not very high, and is close to a critical level, in good agreement with field observations (Begier et al., 2006, Kool, 2005).

Given that an observed outbreak represents only a single sample path profile of all possible sample paths, there are certain diseases for which the occurrence of a major epidemic is rare, with most of the available data being from minor outbreaks. The present study shows that this can happen even when the underlying process is supercritical (R >1), but the probability of major epidemic is small. The pneumonic plague most likely illustrates this case. Even though most historical data arise from observed small outbreaks, what is more important in public health is the assessment of the transmission potential of a major epidemic, such as that occurring in Manchuria. However, large-scale outbreaks of plague cannot avoid contamination by the flea–rodent cycle, and thus, if our interest is limited to precise estimation of human-to-human transmission potential, large epidemics may not be useful. Clearly, pneumonic plague is a disease that required an estimation of the reproduction number using the proposed method. Our study indicates that, given only minor outbreaks, epidemiological assessment of the transmission potential is not sufficiently achieved by directly counting the number of offspring in the transmission tree (because the reproduction number to be estimated is likely to be below 1 given minor outbreaks, even in the case of supercritical process). As long as our objective is to assess the transmission potential, minor outbreak data need to be interpreted with caution adopting assumptions of both subcritical and supercritical processes.

Assuming the dispersion parameter as a known quantity is not regarded as a significant weakness of the proposed approach for a number of reasons. First, data from of transmission network can be used to derive a prior estimate and plausible range for the dispersion parameter k. Second, given a plausible range of k, we can assess the sensitivity of the reproduction number to different k, and check the robustness of our findings. Third, if assuming k as a known quantity is regarded as implausible (or impractical), one should be able to employ a special case, i.e. geometric, to the offspring distribution and (perhaps) conservatively estimate the reproduction number.

The reproduction number was estimated to be as high as 1.59–2.96 based on early growth data from the major epidemic in Manchuria, 1911. There are several explanations for this potential overestimation. First, although the majority of cases were believed as primary pneumonic plague in Manchuria (due to absence of bubo), Y. pestis was isolated from dead rodents at the time (Temporary Quarantine Section, Kanto Totokufu, 1912). The outbreak undoubtedly involved bubonic plague and thus transmission routes other than human-to-human transmission during the course of the epidemic. Consequently, the growth rate is greater than that based only on human-to-human transmission, illustrating the difficulty in estimating R from major epidemic data. Second, validity of R based on epidemic growth rates is vulnerable to various epidemiological characteristics including heterogeneity and mobility of host population (Mercer et al., 2011, Nishiura et al., 2010, Nishiura and Kamiya, 2011). Indeed, as described by Nishiura (2006), the data in Fig. 2A was extracted from the Japanese area only, probably ignoring dynamics in surrounding geographical areas in Manchuria. It remains that the exponential growth rate is useful epidemiological information, offering important insights into estimating R . However, we must be mindful of the relevant pitfalls in the interpretation of the growth rate, and perhaps validate the estimate using other methods, as outlined in the present study.

As the probability of extinction could be higher in heterogeneously mixing populations (Nishiura et al., 2011b), the proposed method could have underestimated R . Our model accounts for dispersion (Garske and Rhodes, 2008), but the negative binomially distributed offspring distribution represents only contact heterogeneity at individual levels. In reality, contact heterogeneity depends on the contact network structure in a population (Trapman, 2007), often elevating the probability of extinction as compared with models with approximate network structures. Moreover, we assumed that the index case was in the population from the beginning of its infectiousness, although the involvement of imported cases in secondary transmission tends to be delayed (Nishiura and Inaba, 2011, Nishiura and Roberts, 2010, Roberts and Nishiura, 2011). Thus, the potential underestimation of R should be considered a limitation of the proposed method.

Because of the need to interpret minor outbreak data carefully, we believe that the proposed method should always be chosen as one of estimation methods given minor outbreak data and would be very useful for estimating the transmission potential, especially for rare diseases. When we observe only minor outbreaks, it must be remembered that the absence of a major epidemic does not always guarantee that the pathogen's transmission potential is subcritical (Begier et al., 2006), and in such an instance, assessment of the transmission potential cannot rely only on directly counting the number of offspring. Thus, the proposed method allows us to conservatively regard the epidemic process leading to minor outbreaks as supercritical and produce the supercritical threshold parameter estimate. As long as we are not sure if the reproduction number is below unity, the investigation of minor outbreak data should always be accompanied by the proposed analysis. We believe that our estimation method has high potential for application to other data sets and we plan to extend the method to the analysis of minor outbreak data in heterogeneously mixing populations.