# Stochasticity and the limits to confidence when estimating ${\cal{R}}_{0}$
of Ebola and other emerging infectious diseases

Bradford P. Taylor  School of Physics, Georgia Institute of Technology,
Atlanta, GA, USA  Corresponding author<EMAIL_ADDRESS>Jonathan
Dushoff  Department of Biology and Institute for Infectious Disease Research,
McMaster University, Hamilton, Canada  Joshua S. Weitz  School of Physics,
Georgia Institute of Technology, Atlanta, GA, USA  School of Biology, Georgia
Institute of Technology, Atlanta, GA, USA

###### Abstract

Dynamic models - often deterministic in nature - were used to estimate the
basic reproductive number, ${\cal{R}}_{0}$, of the 2014-5 Ebola virus disease
(EVD) epidemic outbreak in West Africa. Estimates of ${\cal{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 ${\cal{R}}_{0}$,
and, in turn, the type and scope of necessary interventions. In this
manuscript we propose a hybrid stochastic-deterministic method to estimate
${\cal{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 ${\cal{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 ${\cal{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 ${\cal{R}}_{0}$ CI
at the onset of an epidemic due to multiple sources of noise.

##  1 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,
${\cal{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 ${\cal{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 ${\cal{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 ${\cal{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 ${\cal{R}}_{0}$ estimates for EVD?

Aaron King and colleagues posed a similar question and cautioned that reported
confidence intervals (CIs) of ${\cal{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
${\cal{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 ${\cal{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 ${\cal{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
${\cal{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 ${\cal{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.

##  2 Methods

###  2.1 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, ${\cal 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, $\tau_{c}$. These quantities are
closely related to the doubling time
$\tau_{2}=\tau_{c}\ln{2}=\frac{\ln{2}}{r_{0}}$. These two metrics of disease
outbreak intensity, ${\cal{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)=\frac{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}=({\cal{R}}_{0}-1)/T$. If the generation-
interval distribution is known, an estimate of ${r}_{0}$ implies a _unique_
corresponding estimate of ${\cal R}_{0}$, e.g., it is ${\cal{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 ${\cal{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 ${\cal R}_{0}$. This identifiability problem linking the
measured value of ${r}_{0}$ and the estimated value of ${\cal{R}}_{0}$ hampers
efforts to estimate the potential scope of the disease outbreak over the long-
term [25, 6].

Uncertainty in generation intervals is one of many challenges in estimation of
${\cal{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 ${\cal{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 ${\cal{R}}_{0}$ estimate from a
single stochastic trajectory.

The hybrid stochastic-deterministic approach to estimate the CIs of
${\cal{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 ${\cal{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, $\hat{\tau}_{c}$, given either a CCC or
an Incidence Case Curve (ICC).

  * •

A value of ${\cal{R}}_{0}$ is defined to be compatible with a data set if the
value of ${\tau}_{c}$ inferred from the data, $\hat{\tau}_{c}$, lies in the
middle portion of the distribution of values of $\tau_{c}$ inferred from the
ensemble of simulations associated with that value, $\tilde{\tau}_{c}$.
Throughout this paper, we focus on the middle $95\%$ of the distribution, and
use the symbols $\mbox{}~{}\hat{\mbox{}}$ and $\mbox{}~{}\tilde{\mbox{}}$ to
denote estimates from empirical data and simulations, respectively.

  * •

We estimate CIs associated with ${\cal{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 ${\cal{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 $\tilde{\tau}_{c}$ of stochastic epidemic trajectories in a
SEIRD framework.

###  2.2 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:

$\displaystyle\frac{dS}{dt}$ $\displaystyle=$
$\displaystyle-\beta_{I}SI/N-\beta_{D}SD/N,$ (1) $\displaystyle\frac{dE}{dt}$
$\displaystyle=$ $\displaystyle\beta_{I}SI/N+\beta_{D}SD/N-E/T_{E},$ (2)
$\displaystyle\frac{dI}{dt}$ $\displaystyle=$ $\displaystyle E/T_{E}-I/T_{I},$
(3) $\displaystyle\frac{dR}{dt}$ $\displaystyle=$ $\displaystyle(1-f)I/T_{I},$
(4) $\displaystyle\frac{dD}{dt}$ $\displaystyle=$ $\displaystyle
fI/T_{I}-D/T_{D}.$ (5)

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 $\beta_{D}=0.20$ and $\beta_{I}=0.25$, meaning that a
fraction $\rho_{D}=\frac{{\cal{R}}_{0}(\mathrm{dead})}{{\cal{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 $\tau_{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\tau_{c}$. For example, in the case of $\tau_{c}=21$ then
the fitting period is 42 days in duration. The measured epidemic growth rate,
$\hat{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 $\tau_{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,
$\tau_{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.

Figure 1: Stochastic realizations of a SEIRD epidemic include substantial
variability in epidemic growth rate. (Left) Each panel denotes a randomly
chosen trajectory for which the cumulative case count exceeded 50 (gray
period) followed by a 42 day measurement period (black period). The estimated
$\tilde{\tau}_{c}$ in the measurement period is denoted in the upper left of
each panel. (Top Right) Variation in cumulative case counts once a threshold
is reached. The results are for $5817$ epidemics where the measurement time
denotes the period after the first day the trigger was crossed. The solid blue
line denotes the median number of cases and the shaded region the central 95%
of simulations. (Bottom Right) Variation in the characteristic time estimated
from simulation data, $\tilde{\tau}_{c}$. In all cases, simulations correspond
to stochastic simulations of a SEIRD model in which $N=10^{6}$,
$\beta_{I}=0.25$, $\beta_{D}=0.2$, $\sigma=1/11$, $\gamma=1/6$, $f=0.7$ and
$\rho_{D}=0.25$. The theoretically expected characteristic time is
$\tau_{c}=21$ days. The median characteristic time is $\tilde{\tau}_{c}=20.73$
days.

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 ${\cal{R}}_{0}$.

###  2.3 The expected uncertainty in ${\cal{R}}_{0}$ for an EVD-like outbreak
as estimated from a single stochastic trajectory

Here, we estimate the uncertainty in ${\cal{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 $\rho_{D}=0.25$. Then, $\beta_{I}$ and
$\beta_{D}$ are varied to yield different expected characteristic times for
the case counts, $\tau_{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 $\tilde{\tau}_{c}$ given a range of theoretically expected values of
$\tau_{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
${\cal{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, $\tilde{\tau}_{c}$, given variation in the theoretically
expected $\tau_{c}$. The variation is summarized as a probability distribution
$p(\tilde{\tau}_{c}|\tau_{c})$. The central 95% of the distribution of
$p(\tilde{\tau}_{c}|\tau_{c})$ covers a range whose relative magnitude
increases with $\tau_{c}$.

Figure 2: Fundamental limits on inferring ${\cal{R}}_{0}$ from a single
stochastic epidemic trajectory. Variation in the CIs of theoretical (left)
$\tau_{c}$ and corresponding (right) ${\cal{R}}_{0}$ when using CCCs. See text
for details on methods and interpretation.

The alternative way to interpret Figure 2 is as an inverse problem. In the
inverse problem approach, a measurement of any given $\tilde{\tau}_{c}$ from a
single trajectory is compatible with a range of theoretically expected
$\tau_{c}$ values. For example, given measurement of $\tilde{\tau}_{c}=20$
days, then the associated CI of compatible $\tau_{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(\tilde{\tau}_{c}|\tau_{c})$
across different values of $\tau_{c}$. In this example, we estimate the CI of
$\tau_{c}$ to be 16.6–27.9 days given a single measurement of
$\tilde{\tau}_{c}=20$ days.

The CIs for $\tau_{c}$ from the case data is directly translated to CIs for
${\cal{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 $\tau_{c}$ and specifies the initial deterministic
dynamics. Converting the $\tau_{c}$ CI to ${\cal{R}}_{0}$ CI is shown on the
right panel of Figure 2. For the mock case data with $\hat{\tau}_{c}=20$ we
obtain ${\cal{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.

###  2.4 Confidence intervals for ${\cal{R}}_{0}$ for the EVD outbreak in W.
Africa

We now apply the hybrid stochastic-deterministic approach to estimate the
confidence limits in ${\cal{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 $\tau_{c}$. Here, we vary rates of
transmission as well as the proportion of post-death transmission, $\rho_{D}$.
Due to identifiability issues, measured values of $\hat{\tau}_{c}$ can
correspond to different ${\cal R}_{0}$ depending on the underlying parameters
of the model [25, 6]. Hence, we systematically vary ${\cal R}_{0}$ between
ensembles rather than first varying $\tau_{c}$ and then transforming these
into ranges of ${\cal 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 $\rho_{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, $\hat{\tau}_{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.

Figure 3: Estimated ${\cal{R}}_{0}$ CIs determined by the intersection of the
measured case data $\hat{\tau}_{c}$ with the measured $\hat{\tau}_{c}$ CIs
across a range of theoretical ${\cal{R}}_{0}$ for each country. (Left column)
Characteristic times estimated from case data are given by the height of the
red line on the y-axis and the intervals are determined by the intersection of
the red line with boundary of the CIs projected onto the x-axis. The blue line
refers to the median of the distributions. (Top) Guinea, (Middle) Liberia,
(Bottom) Sierra Leone. Table 2 shows the quantitative results for the CIs.
(Middle Column) Cumulative case data with projections based on ${\cal{R}}_{0}$
CIs. The blue line corresponds to the estimate based fitting the case data.
The gray area delineates the middle 95% of the projection. The green diamonds
are the reported cumulative case data. (Top) Guinea (Middle) Liberia and
(Bottom) Sierra Leone. (Right Column) Incident case data with projections
based on ${\cal{R}}_{0}$ CIs. The blue line corresponds to the estimate based
fitting the case data. The gray area delineates the middle 95% of the
projection. The green diamonds are the reported incident case data. (Top)
Guinea (Middle) Liberia and (Bottom) Sierra Leone.

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 ${\cal{R}}_{0}$ we obtain a distribution of
measured $\tau_{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 ${\cal{R}}_{0}$. Country-specific
estimates of ${\cal{R}}_{0}$ and associated CIs are shown in Table 2. Again,
these CIs represent the combined uncertainty of estimating ${\cal{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 ${\cal{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.

Country | $\mathcal{R}_{0}$ |  $\mathcal{R}_{0}$ CI
---|---|---
Liberia | $2.06$ | $1.93-2.27$
Sierra Leone | $1.71$ | $1.40-1.82$
Guinea | $1.24$ | $1.04-1.42$
Table 1: Confidence intervals of the SEIRD model given EVD case data in
Liberia, Sierra Leone, and Guinea. The ${\cal{R}}_{0}$ values are obtained
from fitting a Poisson regression to the case data. The CIs are calculated
based on Poisson fits to an ensembles of $10^{4}$ stochastic trajectories. The
models used to obtain the trajectories varied in uniformly in $\rho_{D}$
between [.1, .4].

##  3 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 ${\cal{R}}_{0}$, at the early
stages of an epidemic, and showed that stochastic variation constrains the
extent to which confidence limits in ${\cal{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
${\cal{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 ${\cal{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
${\cal{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 ${\cal{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 ${\cal{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
$\tau_{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
${\cal{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 ${\cal{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 ${\cal{R}}_{0}$
in epidemics. In practice, the 95% confidence intervals in ${\cal{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.

## References

  * [1] C. L. Althaus, Estimating the reproduction number of ebola virus (ebov) during the 2014 outbreak in west africa, PLoS Currents Outbreaks, 6 (2014).
  * [2] S. Cauchemez and N. M. Ferguson, Likelihood-based estimation of continuous-time epidemic models from time-series data: application to measles transmission in london, Journal of The Royal Society Interface, 5 (2008), pp. 885–897.
  * [3] G. Chowell and H. Nishiura, Transmission dynamics and control of Ebola virus disease (evd): a review, BMC Medicine, 12 (2014), p. 196.
  * [4] G. Chowell, L. Simonsen, C. Viboud, and Y. Kuang, Is west africa approaching a catastrophic phase or is the 2014 ebola epidemic slowing down? different models yield different answers for liberia, PLoS Currents Outbreaks, 6 (2014).
  * [5] J. M. Drake, R. Kaul, L. W. Alexander, S. M. O’Regan, A. M. Kramer, J. T. Pulliam, M. J. Ferrari, and A. W. Park, Ebola cases and health system demand in liberia, PLoS Biol, 13 (2015), p. e1002056.
  * [6] M. C. Eisenberg, J. N. Eisenberg, J. P. D’Silve, E. V. Wells, S. Cherng, Y.-H. Kao, and R. Meza, Forecasting and uncertainty in modeling the 2014-2015 Ebola epidemic in West Africa, http://arxiv.org/pdf/1501.05555v3.pdf, (2015).
  * [7] D. Fisman, E. Khoo, and A. Tuite, Early epidemic dynamics of the west african 2014 ebola outbreak: estimates derived with a simple two-parameter model, PLoS Currents Outbreaks, 6 (2014).
  * [8] G. J. Gibson and E. Renshaw, Estimating parameters in stochastic compartmental models using markov chain methods, Mathematical Medicine and Biology, 15 (1998), pp. 19–40.
  * [9] D. T. Gillespie, Exact stochastic simulation of coupled chemical reactions, J. Phys. Chem., 82 (1977), pp. 2340–2361.
  * [10] M. F. Gomes, A. Pastore y Piontti, L. Rossi, D. Chao, I. Longini, M. E. Halloran, and A. Vespignani, Assessing the international spreading risk associated with the 2014 west african ebola outbreak, PLoS Currents Outbreaks, 6 (2014).
  * [11] R. N. Gutenkunst, J. J. Waterfall, F. P. Casey, K. S. Brown, C. R. Myers, and J. P. Sethna, Universally sloppy parameter sensitivities in systems biology models, PLoS computational biology, 3 (2007), p. e189.
  * [12] E. L. Ionides, C. Bret , and A. A. King, Inference for nonlinear dynamical systems, Proceedings of the National Academy of Sciences, 103 (2006), pp. 18438–18443.
  * [13] M. J. Keeling and P. Rohani, Modeling Infectious Diseases in Humans and Animals, Princeton University Press, Princeton, NJ, 2007.
  * [14] A. A. King, M. Domenech de Cellès, F. M. G. Magpantay, and P. Rohani, Avoidable errors in the modelling of outbreaks of emerging pathogens, with special reference to ebola, Proceedings of the Royal Society of London B: Biological Sciences, 282 (2015).
  * [15] J. Legrand, R. F. Grais, P.-Y. Boelle, A.-J. Valleron, and A. Flahault, Understanding the dynamics of ebola epidemics, Epidemiology and infection, 135 (2007), pp. 610–621.
  * [16] J. A. Lewnard, M. L. N. Mbah, J. A. Alfaro-Murillo, F. L. Altice, L. Bawo, T. G. Nyenswah, and A. P. Galvani, Dynamics and control of ebola virus transmission in montserrado, liberia: a mathematical modelling analysis, The Lancet Infectious Diseases, 14 (2014), pp. 1189–1195.
  * [17] J. Ma, J. Dushoff, B. M. Bolker, and D. J. Earn, Estimating initial epidemic growth rates, Bull. Math. Biol., 76 (2014), pp. 245–260.
  * [18] M. I. Meltzer, C. Y. Atkins, S. Santibanez, B. Knust, B. W. Petersen, E. D. Ervin, S. T. Nichol, I. K. Damon, and M. L. Washington, Estimating the future number of cases in the Ebola epidemic-Liberia and Sierra Leone, 2014–2015, MMWR Surveill Summ, 63 (2014), pp. 1–14.
  * [19] S. Merler, M. Ajelli, L. Fumanelli, M. F. Gomes, A. P. y Piontti, L. Rossi, D. L. Chao, I. M. Longini, M. E. Halloran, and A. Vespignani, Spatiotemporal spread of the 2014 outbreak of ebola virus disease in liberia and the effectiveness of non-pharmaceutical interventions: a computational modelling analysis, The Lancet Infectious Diseases, 15 (2015), pp. 204–211.
  * [20] C. F. Nielsen, S. Kidd, A. Sillah, E. Davis, J. Mermin, and P. H. Kilmarx, Improving burial practices and cemetery management during an ebola virus disease epidemic-sierra leone, 2014, MMWR Surveill Summ, 64 (2015), pp. 1–8.
  * [21] H. Nishiura and G. Chowell, Early transmission dynamics of ebola virus disease (evd), west africa, march to august 2014, Euro Surveill, 19 (2014), p. 20894.
  * [22] A. Pandey, K. E. Atkins, J. Medlock, N. Wenzel, J. P. Townsend, J. E. Childs, T. G. Nyenswah, M. L. Ndeffo-Mbah, and A. P. Galvani, Strategies for containing ebola in west africa, Science, 346 (2014), pp. 991–995.
  * [23] C. M. Rivers, E. T. Lofgren, M. Marathe, S. Eubank, and B. L. Lewis, Modeling the impact of interventions on an epidemic of ebola in sierra leone and liberia, PLoS Currents Outbreaks, 6 (2014).
  * [24] J. Wallinga and M. Lipsitch, How generation intervals shape the relationship between growth rates and reproductive numbers, Proceedings of the Royal Society B: Biological Sciences, 274 (2007), pp. 599–604.
  * [25] J. S. Weitz and J. Dushoff, Modeling post-death transmission of ebola: Challenges for inference and opportunities for control, Scientific reports, 5 (2015).
  * [26] WHO, http://apps.who.int/gho/data/node.ebola-sitrep, 2014.  [Online; accessed 17-December-2014].
  * [27] WHO Ebola Response Team, Ebola virus disease in West Africa – the first 9 months of the epidemic and forward projections, New England Journal of Medicine, 371 (2014), pp. 1481–1495.

##  Appendix A Appendix

###  A.1 Generating function approach to link characteristic times and
transmission rates

The generalized SEIRD model includes $n_{e}$ number of exposed classes each of
duration $T_{e}/n_{e}$. Given a value of $r_{0}=\tau_{c}^{-1}$ and $\rho_{D}$,
then the generating function approach of Wallinga and Lipsitch [24, 25] can be
used to derive the following relationships:

$\displaystyle{\cal{R}}_{0}$
$\displaystyle=\frac{1+T_{i}r_{0}}{\left(\frac{n_{e}}{n_{e}+T_{e}r_{0}}\right)^{n_{e}}\left(1-\rho_{D}+\frac{\rho_{D}}{1+T_{d}r_{0}}\right)},$
$\displaystyle\beta_{I}$
$\displaystyle=\frac{(1+T_{d}r_{0})(1+T_{i}r_{0})(\rho_{D}-1)}{\left(\frac{n_{e}}{n_{e}+T_{e}r_{0}}\right)^{n_{e}}T_{i}(T_{d}r_{0}(\rho_{D}-1)-1)},$
$\displaystyle\beta_{D}$
$\displaystyle=\frac{(1+T_{d}r_{0})(1+T_{i}r_{0})(\rho_{D})}{\left(\frac{n_{e}}{n_{e}+T_{e}r_{0}}\right)^{n_{e}}fT_{d}(T_{d}r_{0}(\rho_{D}-1)-1)}.$

###  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:

Process | Reaction | rate/probability
---|---|---
Pre-death infection | $S+I\rightarrow E+I$ | $r_{1}=\beta_{I}S\frac{I}{N}$
Post-death infection | $S+D\rightarrow E+D$ | $r_{2}=\beta_{D}I\frac{S}{N}$
Onset of infectiousness | $E\rightarrow I$ | $r_{3}=\sigma=1/T_{E}$
End of infectiousness (survival) | $I\rightarrow R$ | $r_{4}=(1-f)\gamma=(1-f)/T_{I}$
End of infectiousness (death) | $I\rightarrow D$ | $r_{5}=f\gamma=f/T_{I}$
Burial | $D\rightarrow B$ | $r_{6}=\rho=1/T_{D}$
Table 2: Discrete stochastic reactions that transition the state of an
individual

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,
$\mathbf{y}=(N_{0}-1,0,1,0,0)$ at $t=0$. The total rate of outbreak-associated
events is $r_{tot}=\sum_{i=1}^{6}r_{i}$. The time until next event is
determined randomly such that $\delta t\sim\frac{-\log{\chi}}{r_{tot}}$ where
$\chi$ 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.

Figure 4: The distribution of the characteristic time, $\tau_{c}$ increases
with decreasing trigger population. Ensembles of SEIRD stochastic dynamics are
simulated with parameters as in Figure 2. The ensemble includes $10^{4}$
trajectories with an epidemic that persists for 300 days. In each case,
estimates of $\tau_{c}$ are based on fits to the CCC for 42 days after the
trigger population is reached.

###  A.3 Pseudo-Bayesian approach for uncertainty in $\rho_{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
$\rho_{D}$ as uncertain with a uniform prior distribution between 0.1-0.4. We
are interested with models that have a deterministic ${\cal{R}}_{0}$. For each
${\cal{R}}_{0}$, we randomly, uniformly sample $\rho_{D}$ from the prior
distribution. Each sample of ${\cal{R}}_{0}$ and $\rho_{D}$ determine the
growth rate $r_{0}=\tau_{c}^{-1}$ which in turn determines the infection
parameters, $\beta_{I}$ and $\beta_{D}$.

An ensemble of $10^{4}$ trajectories are obtained for each value
${\cal{R}}_{0}$ with $\rho_{D}$ sampled uniformly from [.1 .4] for each
trajectory. The resulting distribution of measured characteristic times
$\hat{\tau}_{c}$ can be interpreted as a marginal distribution across
$\rho_{D}$. Even with a fixed ${\cal{R}}_{0}$, as $\rho_{D}$ varies the
secondary infection distribution changes. Hence, the marginal distribution
across $\rho_{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 $\tau_{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(\tilde{\tau}_{c}|\tau_{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 $\tau_{c}$ and ${\cal R}_{0}$ CIs. The
resulting $\tau_{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].

Figure 5: Variation in regression-based fits of the characteristic time of
case counts between estimates using cumulative or incidence case count data.
The theoretical characteristic time of the underlying model is $\tau_{c}=20$
as shown by the black dashed vertical line. The ICC distribution has a median
of 20.3 days with $95\%$ CIs of 16.6-28.5. The CCC distribution has a median
of 19.9 days with $95\%$ CIs of 17.0-30.1.

We emphasize that our method does not utilize the quality of any individual
fit to generate the CIs for ${\cal{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.

Figure 6: Distributions of standard error of estimated growth rates,
$r_{0}=\hat{\tau}_{c}^{-1}$ for CCCs and ICCs arising from linearly fitting
simulated data assuming deviations are Poisson distributed. The underlying
deterministic model for simulated data has a theoretical characteristic time
of $\tau_{c}=20$ days.

###  A.5 Identifiability problem persists when inferring relative fraction of
post-death transmission from stochastic trajectories

The measured growth rate of an epidemic, $\hat{r}_{0}=\hat{\tau}_{c}^{-1}$ can
be used to infer the basic reproductive number, ${\cal{R}}_{0}$. Using a
generating function approach, it can be shown that ${\cal{R}}_{0}=1/M(-r_{0})$
where $M(z)=\int_{0}^{\infty}e^{za}g(a){\rm d}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 ${\cal{R}}_{0}$ may be
compatible with a single estimate of $\hat{r}_{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<\rho_{D}<1$. Figure 7 shows the measured variation in
the characteristic time, $\langle 1/\hat{r}_{0}\rangle$ 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 $\rho_{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 $\rho_{D}$. Therefore, efforts to constrain estimates of
${\cal{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.

Figure 7: 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
$\rho_{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.

##  Appendix 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.


* Alam et al. (2009) Alam, M-R., Liu, Y. & Yue, D. K. P. 2009 Bragg resonance of waves in a two-layer fluid propagating over bottom ripples. Part I. Perturbation analysis. J. Fluid Mech. 624, 191–224.
  * Baker & Sutherland (2020) Baker, L. E. & Sutherland, B. R. 2020 The evolution of superharmonics excited by internal tides in non-uniform stratification. J. Fluid Mech. 891, R1.
  * Buhler & Holmes-Cerfon (2011) Buhler, O. & Holmes-Cerfon, M. 2011 Decay of an internal tide due to random topography in the ocean. J. Fluid Mech. 678, 271–293.
  * Burns et al. (2020) Burns, K. J., Vasil, G. M., Oishi, J. S., Lecoanet, D. & Brown, B. P. 2020 Dedalus: A flexible framework for numerical simulations with spectral methods. Phys. Rev. Research 2, 023068.
  * Couston et al. (2017) Couston, L-A., Liang, Y. & Alam, M-R. 2017 Oblique internal-wave chain resonance over seabed corrugations. J. Fluid Mech. 833, 538–562.
  * Craik (1971) Craik, A. D. D. 1971 Non-linear resonant instability in boundary layers. J. Fluid Mech. 50 (2), 393–413.
  * Craik et al. (1978) Craik, A. D. D., Adam, J. A. & Stewartson, K. 1978 Evolution in space and time of resonant wave triads - 1. the ’pump-wave approximation’. Proc. Roy. Soc. A 363 (1713), 243–255.
  * Davis & Acrivos (1967) Davis, R. E. & Acrivos, A. 1967 The stability of oscillatory internal waves. J. Fluid Mech. 30 (4), 723–736.
  * Grimshaw (1988) Grimshaw, R. 1988 Resonant wave interactions in a stratified shear flow. J. Fluid Mech. 190, 357–374.
  * Grimshaw (1994) Grimshaw, R. 1994 Resonant wave interactions near a critical level in a stratified shear flow. J. Fluid Mech. 269, 1–22.
  * Grisouard et al. (2011) Grisouard, N., Staquet, C. & Gerkema, T. 2011 Generation of internal solitary waves in a pycnocline by an internal wave beam: a numerical study. J. Fluid Mech. 676, 491–513.
  * Gururaj & Guha (2020) Gururaj, S. & Guha, A. 2020 Energy transfer in resonant and near-resonant internal wave triads for weakly non-uniform stratifications. Part 1. Unbounded domain. J. Fluid Mech. 899, A6.
  * Hall et al. (2013) Hall, R. A., Huthnance, J. M. & Williams, R. G. 2013 Internal wave reflection on shelf slopes with depth-varying stratification. J. Phys. Oceanogr. 43 (2), 248 – 258.
  * Hasselmann (1967) Hasselmann, K. 1967 A criterion for nonlinear wave stability. J. Fluid Mech. 30 (4), 737–739.
  * Ince (1956) Ince, E. L. 1956 Ordinary differential equations. New York, NY: Dover.
  * Kafiabad et al. (2019) Kafiabad, H. A., Savva, M. A. C. & Vanneste, J. 2019 Diffusion of inertia-gravity waves by geostrophic turbulence. J. Fluid Mech. 869, R7.
  * Kirby (1986) Kirby, J. T. 1986 A general wave equation for waves over rippled beds. J. Fluid Mech. 162, 171–186.
  * Klymak et al. (2013) Klymak, J. M., Buijsman, M., Legg, S. & Pinkel, R. 2013 Parameterizing surface and internal tide scattering and breaking on supercritical topography: The one- and two-ridge cases. J. Phys. Oceanogr. 43 (7), 1380 – 1397.
  * Koudella & Staquet (2006) Koudella, C. R. & Staquet, C. 2006 Instability mechanisms of a two-dimensional progressive internal gravity wave. J. Fluid Mech. 548, 165–196.
  * Lahaye & Llewellyn Smith (2020) Lahaye, N. & Llewellyn Smith, S. G. 2020 Modal analysis of internal wave propagation and scattering over large-amplitude topography. J. Phys. Oceanogr. 50 (2), 305–321.
  * Legg (2014) Legg, S. 2014 Scattering of low-mode internal waves at finite isolated topography. J. Phys. Oceanogr. 44 (1), 359 – 383.
  * Legg & Adcroft (2003a) Legg, S. & Adcroft, A. 2003a Internal wave breaking at concave and convex continental slopes. J. Phys. Oceanogr. 33 (11), 2224–2246.
  * Legg & Adcroft (2003b) Legg, S. & Adcroft, A. 2003b Internal wave breaking at concave and convex continental slopes. J. Phys. Oceanogr. 33 (11), 2224 – 2246.
  * Li & Mei (2014) Li, Y. & Mei, C. C. 2014 Scattering of internal tides by irregular bathymetry of large extent. J. Fluid Mech. 747, 481–505.
  * Liang et al. (2017) Liang, Y., Zareei, A. & Alam, M-R 2017 Inherently unstable internal gravity waves due to resonant harmonic generation. J. Fluid Mech. 811, 400–420.
  * MacKinnon et al. (2013) MacKinnon, J. A., Alford, M. H., Sun, O., Pinkel, R., Zhao, Z. & Klymak, J. 2013 Parametric subharmonic instability of the internal tide at 29°N. J. Phys. Oceanogr. 43 (1), 17 – 28.
  * MacKinnon & Winters (2005) MacKinnon, J. A. & Winters, K. B. 2005 Subtropical catastrophe: Significant loss of low-mode tidal energy at 28.9 degrees. Geophys. Res. Lett. 32 (15).
  * Mathur et al. (2014) Mathur, M., Carter, G. S. & Peacock, T. 2014 Topographic scattering of the low-mode internal tide in the deep ocean. J. Geophys. Res. Oceans 119 (4), 2165–2182.
  * Maugé & Gerkema (2008) Maugé, R. & Gerkema, T. 2008 Generation of weakly nonlinear nonhydrostatic internal tides over large topography: a multi-modal approach. Nonlinear Process. Geophys. 15 (2), 233–244.
  * Meyer (1979) Meyer, R. E. 1979 Surface wave reflection by underwater ridges. J. Phys. Oceanogr. 9 (1), 150 – 157.
  * Munk (1966) Munk, W. H. 1966 Abyssal recipes. Deep-Sea Res. Oceanogr. Abstr. 13 (4), 707 – 730.
  * Olbers et al. (2020) Olbers, D., Pollmann, F. & Eden, C. 2020 On PSI interactions in internal gravity wave fields and the decay of baroclinic tides. J. Phys. Oceanogr. 50 (3), 751 – 771.
  * Phillips (1966) Phillips, O. M. 1966 The dynamics of the upper ocean. ambridge University Press.
  * Sutherland (2016) Sutherland, B. R. 2016 Excitation of superharmonics by internal modes in non-uniformly stratified fluid. J. Fluid Mech. 793, 335–352.
  * Varma et al. (2020) Varma, D., Chalamalla, V. K. & Mathur, M. 2020 Spontaneous superharmonic internal wave excitation by modal interactions in uniform and nonuniform stratifications. Dyn. Atmospheres Oceans 91, 101159.
  * Varma & Mathur (2017) Varma, D. & Mathur, M. 2017 Internal wave resonant triads in finite-depth non-uniform stratifications. J. Fluid Mech. 824, 286–311.
  * Vic et al. (2019) Vic, C., Garabato, A.C.N., Green, J.A.M., Waterhouse, A.F., Zhao, Z., Mélet, A., De Lavergne, C., Buijsman, M. C. & Stephenson, G. R. 2019 Deep-ocean mixing driven by small-scale internal tides. Nat. Commun. 10 (1), 1–9.
  * Voelker et al. (2021) Voelker, G. S., Akylas, T. R. & Achatz, U. 2021 An application of WKBJ theory for triad interactions of internal gravity waves in varying background flows. Q. J. R. Meteorol. Soc. 147 (735), 1112–1134.
  * Wunsch (2017) Wunsch, S. 2017 Harmonic generation by nonlinear self-interaction of a single internal wave mode. J. Fluid Mech. 828, 630–647.
  * Young et al. (2008) Young, W. R., Tsang, Y.-K. & Balmforth, N. J. 2008 Near-inertial parametric subharmonic instability. J. Fluid Mech. 607, 25–49.
  * Zhao et al. (2016) Zhao, Z., Alford, M. H., Girton, J. B., Rainville, L. & Simmons, H. L. 2016 Global Observations of Open-Ocean Mode-1 M2 Internal Tides. J. Phys. Oceanogr. 46 (6), 1657–1684.