Stochasticity and the limits to confidence when estimating R_0 of Ebola and other emerging infectious diseases

Dynamic models - often deterministic in nature - were used to estimate the basic reproductive number, R_0, of the 2014-5 Ebola virus disease (EVD) epidemic outbreak in West Africa. Estimates of R_0 were then used to project the likelihood for large outbreak sizes, e.g., exceeding hundreds of thousands of cases. Yet fitting deterministic models can lead to over-confidence in the confidence intervals of the fitted R_0, and, in turn, the type and scope of necessary interventions. In this manuscript we propose a hybrid stochastic-deterministic method to estimate R_0 and associated confidence intervals (CIs). The core idea is that stochastic realizations of an underlying deterministic model can be used to evaluate the compatibility of candidate values of R_0 with observed epidemic curves. The compatibility is based on comparing the distribution of expected epidemic growth rates with the observed epidemic growth rate given"process noise", i.e., arising due to stochastic transmission, recovery and death events. By applying our method to reported EVD case counts from Guinea, Liberia and Sierra Leone, we show that prior estimates of R_0 based on deterministic fits appear to be more confident than analysis of stochastic trajectories suggests should be possible. Moving forward, we recommend including a hybrid stochastic-deterministic fitting procedure when quantifying the full R_0 CI at the onset of an epidemic due to multiple sources of noise.


Introduction
The SEIRD model of Ebola virus disease (EVD) dynamics, introduced by Legrand and colleagues [15], considers the transitions among susceptible, exposed, infectious, recovered and dead (but unburied) individuals. Variants of this core model have been utilized to estimate the basic reproductive number, R 0 , of EVD from incidence and cumulative case data in the 2014-5 epidemic in West Africa (e.g., [1,10,7,27,21,22,16,18,23]). In some instances these estimates include 95% confidence intervals. For example, one of the first mathematical epidemiology papers published in response to the ongoing epidemic in W. Africa estimated the R 0 of EVD to be: 1.51 (95% CI 1.50-1.52) in Guinea; 1.59 (1.57-1.60) in Liberia; and 2.53 (2.41-2.67) in Sierra Leone [1]. The confidence limits appear over-confident, particularly for Guinea and Liberia.
Subsequently, the World Health Organization (WHO) Ebola Response Team analyzed a model with pre-and post-death transmission and estimated R 0 to be: 1.71 (1.44-2.01) for Guinea; 1.83 (1.72-1.94) for Liberia; and 2.02 (1.79-2.26) for Sierra Leone [27]. Similarly, a follow-up model for EVD in the Montserrado region of Liberia estimated R 0 to be 2.49 (2.38-2.60) [16]. Comparison of these case studies shows that models can provide incompatible inferences with non-overlapping CIs despite similar data. Obviously, differences will arise in model structure, data range and quality, and the treatment of noise. Yet, are even these more recent studies over-confident about the precision of R 0 estimates for EVD? Aaron King and colleagues posed a similar question and cautioned that reported confidence intervals (CIs) of R 0 for EVD are likely too narrow whenever they neglect stochasticity in the disease transmission process [14]. Such over-confidence is further heightened by fitting deterministic models to cumulative case curves (CCCs). A CCC is a monotonically increasing function of time, representing the total number of individuals reported to have become infected during an outbreak. Individual time points within CCCs are not independent, and so "error" in fitting deterministic models can appear artefactually low. Such low error in fits -if not properly accounted for -can lead to a misleading interpretation of overly narrow CIs for R 0 . Instead King et al. [14] recommend a stochastic-based fitting procedure to incidence case curves (ICCs), to account for observation noise and process noise into estimates of R 0 and its associated CIs.
Here, we propose a hybrid stochastic-deterministic approach to address similar issues. The key differentiating theoretical feature of our approach from that proposed by King et al. [14], is that we address uncertainty in generation-interval distributions and its effect on R 0 . The generation-interval distribution, g(a), is the normalized fraction of secondary cases caused by an infectious individual at "age" a since infection. In the case of EVD, there was -and is -uncertainty with respect to transmission event times, including that of post-death transmission [25,20].
We apply our hybrid approach retrospectively to estimate the CIs for R 0 for EVD in the summer of 2014, coincident with the release of the first projections for the potential size of the outbreak. We conclude that many early estimates of the CIs for R 0 for EVD were almost certainly over-confident. The over-confidence arose, at least in part, because estimates did not account for the effects of process noise and for uncertainty in disease generation times. More generally, we explain how models of disease transmission can be adapted to our approach -by using a principled filtering method to identify ensembles of simulated stochastic trajectories "compatible" with a single, measured epidemic outbreak. Overall, we propose a flexible framework that is easily adaptable to the specific dynamics of a particular emerging epidemic while remaining simple to computationally implement.

Measuring Compatible Epidemic Trajectories given Stochasticity
Dynamic models provide a means to estimate the overall intensity of an outbreak from measurements of infectious case data. The intensity is usually quantified in terms of the basic reproduction number, R 0 ,-the total number of secondary cases caused by a single infectious individual in an otherwise naive population. An alternative metric of intensity is the growth rate of disease incidence, r 0 , which is the reciprocal of the characteristic time of exponential growth, τ c . These quantities are closely related to the doubling time τ 2 = τ c ln 2 = ln 2 r 0 . These two metrics of disease outbreak intensity, R 0 and r 0 , are related but not equivalent; they are linked by the generation-interval distribution, g(a). For example, in a simple SIR model then g(a) = 1 T e −a/T where the average infectious period is T . As is well known, the theoretically expected epidemic growth rate is r 0 = (R 0 − 1)/T . If the generation-interval distribution is known, an estimate of r 0 implies a unique corresponding estimate of R 0 , e.g., it is R 0 = 1 + r 0 T in the case of a SIR model. Conversely, when the generation-interval distribution is uncertain, then a range of values of R 0 may be compatible with a given rate of increase in disease incidence -in particular, an increase in the estimated mean generation-interval T would lead to an increase in estimated R 0 . This identifiability problem linking the measured value of r 0 and the estimated value of R 0 hampers efforts to estimate the potential scope of the disease outbreak over the longterm [25,6].
Uncertainty in generation intervals is one of many challenges in estimation of R 0 . Other important sources of uncertainty include the intrinsic stochasticity of the epidemic, the case-reporting process and uncertainty in the structural mode of disease transmission. Ignoring or under-estimating uncertainty can lead to over-confident estimates of R 0 . The premise of our approach is that a measured time series of infectious case data represents a single trajectory from an ensemble of stochastic trajectories given a set of underlying and unknown parameters. The observed trajectory of disease incidence need not grow exponentially at theoretically expected values. Rather there are many trajectories that could appear statistically indistinguishable from the observed epidemic outbreak. We term these trajectories: "Compatible Epidemic Trajectories given Stochasticity" or COMETS. The expected variation in the growth rates of COMETS is quantifiable and can be used to bound the confidence of an R 0 estimate from a single stochastic trajectory.
The hybrid stochastic-deterministic approach to estimate the CIs of R 0 from synthetic data relies on an inverse problem approach to regression-based fits. The approach involves the following steps: • A range of deterministic models is considered that vary in disease-associated parameters, including R 0 .
• For each model and fixed parameter set within the range, we simulate an ensemble of stochastic trajectories and utilize a metric to compare the simulated trajectories to the case data. The metric is the regression-based estimate of the characteristic time, τ c , given either a CCC or an Incidence Case Curve (ICC).
• A value of R 0 is defined to be compatible with a data set if the value of τ c inferred from the data,τ c , lies in the middle portion of the distribution of values of τ c inferred from the ensemble of simulations associated with that value,τ c . Throughout this paper, we focus on the middle 95% of the distribution, and use the symbolsˆandt o denote estimates from empirical data and simulations, respectively.
• We estimate CIs associated with R 0 by identifying the range of deterministic models that can yield dynamics compatible with the case data, i.e., by identifying COMETS.
This series of steps can be adopted to any epidemic context. Moreover, the well-known problems with inferring CIs for R 0 based on the quality of model fits to a single CCC do not arise via the COMETS approach. In the next section we illustrate why the method is relevant to EVD by highlighting the extent to which process noise drives uncertainty in the distribution ofτ c of stochastic epidemic trajectories in a SEIRD framework.

Stochastic trajectories of epidemic outbreaks have substantial variation in epidemic growth rates
Stochastic variation in the timing of discrete transmission events can generate variation in estimates of the realized growth rate of an epidemic, r 0 . The amount of variation depends on the particular disease model, disease parameters, and the time-scale over which r 0 is estimated. As an example inspired by EVD, we consider disease dynamics based on a SEIRD model: We assume that the average latent period is T E = 11 days, the average infectious period is T I = 6 days, there is a f = .7 chance of an infected individual ultimately dying and there is an average time of T D = 4 days before burial. The total population is set at N = 10 6 . We set the disease transmission rates to be β D = 0.20 and β I = 0.25, meaning that a fraction ρ D = R 0 (dead) R 0 = 0.25 of transmission is attributable to post-death transmission (see A.1). In a deterministic framework, epidemics that obey such a model given those parameters should increase with a characteristic time of τ c = 21 days [24,25]. Note, that our choice of models here and throughout the paper inherently alter the estimates of the confidence intervals. Thus we caution that identifying an appropriate underlying model structure is a key but separate issue to the analysis process.
We simulate an ensemble of 10 4 stochastic realizations of this SEIRD model beginning with a single infectious individual (see A.2). We consider only realizations that produce at least 50 total cases (approximately 58% of all realizations). This threshold acts as a trigger, from which point we track the time series of incidence, accumulating cases at exact daily intervals. We do not account for reporting error in this example. We estimate the growth rate by fitting an incidence curve (CCC or ICC) from the trigger time t 0 until t f = t 0 + 2τ c . For example, in the case of τ c = 21 then the fitting period is 42 days in duration. The measured epidemic growth rate,r 0 , for a given trajectory is the slope of the best-fit line to log-transformed censused time series assuming errors are Poisson distributed. Phenomenologically these errors are attributed to observation noise. Thus our presented model only implicitly includes observation noise. That is we do not introduce external noise to our data representing errors in reporting. Note that the results depend on the choices for t 0 and τ c . Hence, the modeler should choose t 0 such that the error in the fitting line to the case data is satisfactoraly small. Meanwhile for an emerging epidemic, τ c should be dictated by the total time elapsed from the onset. The demonstration of our methodology does not rely on our specific choices in this paper.
The left panel of Figure 1 shows five examples of stochastic epidemic trajectories. The time from the index case until the point at which the epidemic "takes off" (i.e., reaches a threshold number of cases), is highly variable. The upper right panel of Figure 1 shows that, even after take-off, there is substantial variation in growth rate between different epidemic realizations. The trajectories have been time shifted so that t = 0 is now defined as the day where the trigger population is reached or surpassed. There is over 40% variation above and below the number of infectious cases at t = 42 days. The bottom right panel of Figure 1 shows the variation in characteristic times estimated from realized data, within an ensemble of epidemic trajectories. The median characteristic time is 20.73 days with a CI of [17.0, 30.1] days. In the appendix, we perform a sensitivity analysis of our method by considering different trigger thresholds. There, the distribution of the characteristic time widens as the cumulative trigger population decreases (see Figure 4). This exploratory analysis reveals substantial variability in the growth rate of epidemic trajectories. Quantifying the uncertainty in the characteristic time of epidemic outbreaks in an ensemble is the basis for estimating the CI for R 0 .

The expected uncertainty in R 0 for an EVD-like outbreak as estimated from a single stochastic trajectory
Here, we estimate the uncertainty in R 0 for an EVD-like outbreak with SEIRD dynamics. As in the previous section, the fraction of post-death transmission is assumed to be ρ D = 0.25. Then, β I and β D are varied to yield different expected characteristic times for the case counts, τ c (see analytic relationships in A.1). For each parameter set, we simulate 10 4 stochastic trajectories with the same censusing conditions as in the previous section, conditioned on the fact that the epidemic continues throughout the sampling period. Variation in the measuredτ c given a range of theoretically expected values of τ c is shown in the left panel of Figure 2. This figure is the basis for identifying COMETS and, in turn, for estimating the CIs for R 0 at the onset of the epidemic. The conventional way to interpret Figure 2 is as a "forward problem" -as introduced in the previous section. In the forward problem approach, Figure 2 depicts variation in the measured characteristic time of stochastic epidemic trajectories,τ c , given variation in the theoretically expected τ c . The variation is summarized as a probability distribution p(τ c |τ c ). The central 95% of the distribution of p(τ c |τ c ) covers a range whose relative magnitude increases with τ c .
The alternative way to interpret Figure 2 is as an inverse problem. In the inverse problem approach, a measurement of any givenτ c from a single trajectory is compatible with a range of theoretically expected τ c values. For example, given measurement of τ c = 20 days, then the associated CI of compatible τ c is obtained by scanning horizontally (the red line) for intersections with the forward problem distributions (the black lines). This intersection is estimated by linear interpolation of the 2.5% and 97.5% CI for p(τ c |τ c ) across different values of τ c . In this example, we estimate the CI of τ c to be 16.6-27.9 days given a single measurement ofτ c = 20 days.
The CIs for τ c from the case data is directly translated to CIs for R 0 by utilizing a generating function approach [24]. This method involves determining the moment generating function associated with the distribution of the age of secondary infections resulting from a single infectious individual in an otherwise susceptible population. This distribution varies with τ c and specifies the initial deterministic dynamics. Converting the τ c CI to R 0 CI is shown on the right panel of Figure 2. For the mock case data withτ c = 20 we obtain R 0 = 2.10 with a CI of 1.80 − 2.45. In summary, these uncertainty ranges arise solely due to process noise and were identified by quantifying COMETS in the SEIRD model framework.

Confidence intervals for R 0 for the EVD outbreak in W. Africa
We now apply the hybrid stochastic-deterministic approach to estimate the confidence limits in R 0 for early outbreak dynamics of EVD in Guinea, Liberia and Sierra Leone. The present method differs from that of the previous section in one key way. Previously, we varied rates of transmission to yield models with different theoretical τ c . Here, we vary rates of transmission as well as the proportion of post-death transmission, ρ D . Due to identifiability issues, measured values ofτ c can correspond to different R 0 depending on the underlying parameters of the model [25,6]. Hence, we systematically vary R 0 between ensembles rather than first varying τ c and then transforming these into ranges of R 0 . We construct the ensembles in a pseudo-Bayesian framework. Specifically, we account for uncertainty in the fraction of infections transmitted by the deceased by considering a uniform distribution of values of ρ D between 0.1 − 0.4 (see A.3). In this way, we account for uncertainty in the distribution of times to secondary transmissions in addition to uncertainty from process noise. Additionally, the total population for each country is set to match census data (N = 4.3 * 10 6 for Liberia, N = 6.1 * 10 6 for Sierra Leone, and N = 11.8 * 10 6 for Guinea).
Case data was obtained from the WHO [26]. We use cumulative counts and consider dynamics in each country given a start date at which at least 50 cumulative infections have occurred and a final date of September 7, 2014. We choose this time period to reflect the onset of the epidemic across all three countries. The measured characteristic time,τ c is obtained by fitting an exponential to the CCC for each country after the trigger population is reached. The lowest value of cumulative case counts above 50 for each country is also considered the trigger population for the stochastic simulations. We consider a SEIRD model with a gamma distributed exposed period with n = 2 classes in accordance with previous analysis [27]. The ensemble of simulated stochastic trajectories are conditioned on the epidemic remaining throughout the census period, as was the case with the EVD case data. The fits to the stochastic trajectories are subject to the same conditions as the case data.
Resulting CIs of the characteristic time for Liberia, Sierra Leone, and Guinea are shown in the left column of Figure 3. The characteristic time measured from the case data specifies the location of the red line along the y-axis. For each value of theoretical R 0 we obtain a distribution of measured τ c from 10 4 simulated stochastic trajectories. The confidence intervals for the reproductive number are determined by where the red line intersects the lower and upper limits of central 95% characteristic time distribution for each theoretical R 0 . Country-specific estimates of R 0 and associated CIs are shown in Table 2. Again, these CIs represent the combined uncertainty of estimating R 0 from a single stochastic trajectory given uncertainty in the distribution of times to secondary infection. Hence, they represent a lower bound of uncertainty in R 0 given additional uncertainty arising from observation noise. These CIs are larger than many early estimates, e.g., [1].
In the middle (right) column of Figure 3 we project forward the uncertainty in the cumulative (incident) case counts based on our CIs. We project forward by 8 weeks to show the range of expected case counts over time assuming no further control measures are implemented. The case data lies within the projection for Guinea and Sierra Leone, but are substantially lower for Liberia, likely because transmission parameters had already changed substantially, due to behavior change, control measures, or both, or due to structural differences in the epidemic process [4]. Overall, we expect large uncertainty in projected case counts during the exponential phase of the epidemic, due to process noise and the generation-interval uncertainty alone.

Discussion
A stochastic implementation of an epidemic leads to significant variation in the realized time series of infectious individuals, due to "process noise" [13] -the stochastic sequence of discrete events in which individuals become infected, infect others, and eventually recover or die. The difference between trajectories predicted in deterministic vs. stochastic models has been observed and studied for decades in other disease contexts [8,12,2,17]. Here, we explored the practical consequences of such differences when trying to infer epidemiological parameters, including R 0 , at the early stages of an epidemic, and showed that stochastic variation constrains the extent to which confidence limits in R 0 can be narrowed. In practice, time series that are used for fitting dynamical epidemiological models are drawn from the early stages of an epidemic. At such early stages, a single trajectory may be well fit by a single exponential rate from which R 0 can be estimated. Yet, the trajectories within an ensemble generated given the same underlying parameters will be fit by a distribution of exponential rates. Hence, the 95% CIs for an estimate of R 0 when using a deterministic model are significantly broadened due to process noise.
This general limitation of fitting deterministic models applies to the study of EVD. Multiple groups have proposed model variants of EVD dynamics, fit deterministic models to case data, and then used such fits to estimate R 0 , including associated CIs (e.g., [1,16,22,27]). The time range we used for fitting includes the period from onset as defined by the time with greater than 50 cumulative infections to early September 2014. This range spans from the onset of the epidemic to reported cumulative case counts of nearly 700 in Guinea and over 1500 in both Sierra Leone and Liberia. Using a stochastic implementation of a SEIRD model, we inferred R 0 to be 1.24 for Guinea (1.04-1.42 95% CI), 2.06 for Liberia (1.93-2.27 95% CI), and 1.71 for Sierra Leone (1.40-1.82 95% CI). The estimates and ranges reflect three assumptions: (i) the model structure; (ii) the prior assumptions about the exposed and infectious periods, time to bury and probability of death; (iii) the availability and quality of epidemic case data. In particular, the second assumption involves varying the distribution of time to secondary infection. These CIs denote limits to inference in this class of model imposed by the nature of the data available: a single stochastic epidemic outbreak. The original estimates of CIs from [1] are evidently too narrow, yet even later estimates should be revisited using the hybrid approach proposed here.
The present method is complementary to alternative, profile-likelihood based approaches [14]. We obtain distributions for R 0 by fitting to ensembles obtained by stochastically simulating a range of deterministic models. We marginalized our ensembles so that the trajectories in the ensembles were obtained from underlying models with the same theoretical τ c . Previous work used particle filtering techniques to obtain a distribution based on fitting to the data [14]. This approach estimates uncertainty due to fitting deterministic models to a single, stochastic trajectory. In our case, we find that differences between fitting to cumulative or incident data are negligible, so long as the quality of individual fits within an ensemble is not used as the primary source of information on the CIs (see A.4). We note that our methodology of estimating COMETS is similar to the method of "plausible parameter sets" advocated for use in estimating disease-associated parameters during outbreaks [5]. Our goal is not to replace any previous methods but instead to provide an alternative flexible and computationally simple framework to estimate confidence intervals for an emerging epidemic.
The hybrid approach comes with certain precautions. First off, this method is does not dictate the underlying model structure. Thus the relevance of the results relies on how accurately the actual dynamics of the epidemic are captured by the chosen underlying model. Of course, this is an issue whenever modeling actual epidemics and should be addressed prior to implementing our method. Regardless of specific modeling choices, nonlinear models often generally display sloppiness such that many combinations of parameters have little effect on certain system dynamics [11]. The identifiability problem [25,6] also applies here (see A.5) , so that the strength of fit does not necessarily exclude a range of compatible mechanisms. Implementing the current, hybrid approach should be straightforward to adapt to both well-mixed and spatially-explicit models of disease transmission. For example, a recent analysis of spatial data suggested that apparent exponential dynamics is a result of aggregation of local epidemics [3]. Yet, we expect that high-performance implementations of stochastic models will be required for spatially explicit projections of outbreak sizes and associated CIs for R 0 at the early stages of an epidemic [19].
In summary, we should have expected to have less confidence in estimates of CIs of R 0 at the outset of a EVD epidemic given process noise. We hope that the current method, similar in intent to that of King et al. [14], provides an accessible route for estimating realistic CIs for R 0 in epidemics. In practice, the 95% confidence intervals in R 0 estimated from stochastic model fits will be broader than that estimated from deterministic model fits to cumulative case data. Estimates of CIs using the current hybrid approach represent a lower bound of uncertainty due to stochastic sources of noise. As an epidemic continues and the number of infected individuals increases, observation noise contributes a relatively larger proportion of uncertainty as compared to process noise [17]. Remaining realistic about the limits to confidence in model fits should also be incorporated into public health practice, e.g., when projecting the necessary scope of intervention based on "optimal" fits [18]. We encourage the academic, governmental, and non-governmental public health communities to consider incorporating unavoidable uncertainty into their decision making pipelines when responding to emergent disease outbreaks.

Process
Reaction rate/probability Pre-death infection

A.2 Stochastic simulations of epidemic outbreaks
Stochastic realizations of the SEIRD model are simulated using the Gillespie framework [9], given the "reaction" events in the following table: Processes transition individuals who are susceptible (S), exposed (E), infected (I), deceased (D), recovered (R), and buried (B). The total population is fixed at N = S + E + I + R + B. Epidemics are initiated with one infectious individual in an otherwise susceptible population. Mathematically, the initial state is, in the respective ordering, y = (N 0 − 1, 0, 1, 0, 0) at t = 0. The total rate of outbreak-associated events is r tot = 6 i=1 r i . The time until next event is determined randomly such that δt ∼ − log χ rtot where χ is a uniformly distributed number between 0 and 1. In this way, the time between events follows an exponential distribution with rate r tot . Then, the probability of each event is r i /r tot . After selecting an event and updating the discrete number of individuals, the reaction rates are recalculated and the process continues. The same framework can be extended to include multiple n subclasses within the exposed class, to capture the peaked nature of the exposed period (centered around 11 days for EVD). Trajectories are complete when the epidemic dies out because there are no more infectious individuals. In the present context, we are interested in those trajectories that do not die out before the end of the simulation time.

A.3 Pseudo-Bayesian approach for uncertainty in ρ D
The fraction of infections of EVD in west Africa due to transmissions from the deceased has been estimated to be between .01-0.3 [27]. However an estimate from a prior EVD outbreak in Congo is even higher [15]. To address this treat ρ D as uncertain with a uniform prior distribution between 0.1-0.4. We are interested with models that have a deterministic R 0 . For each R 0 , we randomly, uniformly sample ρ D from the prior distribution. Each sample of R 0 and ρ D determine the growth rate r 0 = τ −1 c which in turn determines the infection parameters, β I and β D .
An ensemble of 10 4 trajectories are obtained for each value R 0 with ρ D sampled uniformly from [. 1 .4] for each trajectory. The resulting distribution of measured characteristic timesτ c can be interpreted as a marginal distribution across ρ D . Even with a fixed R 0 , as ρ D varies the secondary infection distribution changes. Hence, the marginal distribution across ρ D can also be interpreted as a marginal distribution across the time to secondary infection distributions that all correspond to the same reproductive number of the disease in the population τ c .

A.4 Comparing results from CCC and ICC data
There are a number of potential pitfalls of using CCCs rather than ICCs. Yet, King and colleagues in their critique of CCCs noted that the summary statistics of the epidemic growth for SEIR models were functionally equivalent when fitted to an ensemble of CCCs and ICCs given the same underlying disease parameters [14]. The difference, as they pointed out, was how to interpret the quality of the fits in inferring CIs. Similarly, here we find that the resulting distributions, p(τ c |τ c ), measured using either the CCC or ICC are similar but not equivalent ( Figure 5). Differences in the resulting CIs are minor, but the median of the ICC distribution is skewed larger than the theoretical value. These differences scale up to the overall analysis, but with minor effect on the τ c and R 0 CIs. The resulting τ c CIs from our synthetic data are 17.0-30.1 given fits to cumulative data and 16.6-28.5 given fits to incident data. The resulting CIs are approximately the same size. This contrasts with prior results from profile likelihoods in which the CIs as inferred from CCCs are contained within the CIs as inferred from ICCs [14].
We emphasize that our method does not utilize the quality of any individual fit to generate the CIs for R 0 , for precisely the reasons cautioned by King and colleagues [14]. Estimating the growth rate, r 0 , from linear regression is uncertain due to the fitting procedure itself. Associated with each estimate of r 0 are confidence intervals. We compare the distributions of errors due to fitting associated with the CCCs and ICCs of our simulated data in Figure 6. The errors associated with ICCs are larger than those associated with CCCs. It remains an open question as to whether/when CCCs rather than ICCs are preferable when leveraging regression fitting to estimate r 0 given process noise alone, rather than observation noise.

A.5 Identifiability problem persists when inferring relative fraction of post-death transmission from stochastic trajectories
The measured growth rate of an epidemic,r 0 =τ −1 c can be used to infer the basic reproductive number, R 0 . Using a generating function approach, it can be shown that a . [24] where g(a) is the normalized fraction of all secondary cases caused by an infectious individual at "age" a since infection. A range of values of R 0 may be compatible with a single estimate ofr 0 [25]. This uncertainty is a consequence of an identifiability problem given uncertainty in the relative fraction of transmission events that could be attributed to post-death transmission. Here we ask: what is the variation in the growth rate, r 0 , and characteristic time, 1/r 0 , compatible with varying fraction of post-death transmission, 0 < ρ D < 1. Figure 7 shows the measured variation in the characteristic time, 1/r 0 for 5000 ensembles for three different expected growth rates r 0 = 1/14, 1/21 and 1/28 days −1 . In each case, we varied ρ D from 0 to 1 in increments of 0.2. As in the prior section, we find that the characteristic time of an epidemic can vary substantially for a fixed value of r 0 . Here, we also expect that a range of mechanistic models can all yield the same expected characteristic time. As is evident, the identifiability problem highlighted in prior analyses of deterministic models [25,6] also applies in the case of stochastic models. For SEIRD models, the expected value of the basic reproductive number increases with ρ D . Therefore, efforts to constrain estimates of R 0 from EVD case data will be subject to inherent variability due, in part, to uncertainty in mechanism, e.g., the relative fraction of post-death transmission, and process, i.e., stochastic outbreak dynamics.

B Acknowledgments
The work was funded by a grant to JSW from the Burroughs Wellcome Fund and from the Army Research Office grant #W911NF-14-1-0402. We would like to acknowledge Luis F. Jover for reviewing the manuscript.      : Identifiability problem persists when fitting SEIRD models to stochastic data.The three scenarios correspond to cases where the characteristic time, 1/r 0 =14, 21 and 28 days. The realized epidemic growth rates of stochastic trajectories are measured given variation in ρ D from 0 to 1 in increments of 0.1. Circles denote the median characteristic time while triangles denote the 95% confidence intervals from an ensemble of 10 3 simulations per condition.