Estimation of end-of-outbreak probabilities in the presence of delayed and incomplete case reporting

We use a hidden state model in which the instantaneous reproduction number $R_{t}$ at time $t$, the daily incidence of new local infections $I_{t}$ and related quantities are treated as hidden states. The number of new case notifications $C_{t}$ each day is treated as an observed state.

The instantaneous reproduction number $R_{t}$ is assumed to follow a Gaussian random walk in logarithmic space:

where $\delta_{t}$ and $\sigma_{t}$ are the mean and standard deviation of the random walk step on day $t$. We set $\delta_{t}=0$ and $\sigma_{t}$ to be a constant for most of the simulated time period. However, in the two case studies that we consider, in a short time window around the start of the interventions, we set $\delta_{t}<0$. This allows $R_{t}$ to decrease rapidly in response to the interventions.

To allow for individual heterogeneity in transmission rates, we assume that each infected individual is assigned a “transmission multiplier”, representing variation in infectiousness and/or contact rates between individuals. We assume the transmission multipliers are independent, identically distributed random deviates from a Gamma distribution with shape $k_{r}$ and scale $1/k_{r}$ (i.e. mean $1$). Thus the total of the transmission multipliers for people infected on day $t$ is

where $n_{t}$ is the number of imported infections arising on day $t$ and $\varepsilon$ is the average infectivity of imported infections relative to local infections. The variable $Y_{t}$ represents the aggregate infectivity of people who were infected on day $t$.

The number of new local infections on day $t$ follows a standard renewal equation [19] but with the aggregate infectivity of infections $s$ days previously ($Y_{t-s}$) used in place of the number of infections $s$ days previously ($I_{t-s}$) to drive the number of local infections on day $t$:

where $g_{s}$ is the probability mass function for the generation interval distribution. Under this formulation, if the reproduction number was a constant $R$, then the total number of people infected by a randomly selected infected individual would follow a negative binomial distribution with mean $R$ and dispersion parameter $k_{r}$ [16]. In the limit $k_{r}\to\infty$, $Y_{t}$ is deterministically equal to $I_{t}+\varepsilon n_{t}$, meaning that Eq. (3) reduces to the standard Poisson renewal equation. Smaller values of $k_{r}$ correspond to larger variance of $Y_{t}$, representing greater individual heterogeneity in transmission rates.

Note that the formulation of Eq. (3) relies on the assumption that the number of secondary infections arising from each infected individual is independent of other infected individuals. This means that, although the model includes individual heterogeneity, it does not include network-type effects such as a higher probability of highly connected individuals infecting other highly connected individuals. Furthermore, Eq. (3) assumes that the generation interval distribution is known. This is a standard modelling assumption, but we caution that estimating the generation interval distribution from epidemiological data is subject to various biases, particularly during the exponential growth phase of a newly emerging outbreak [22, 23].

To account for delays between infection and case notification, we assign each infected individual a notification time (regardless of whether or not they are actually notified as cases) according to an infection-to-notification distribution with probability mass function $u$. Among individuals infected on a given day $s$, the number of these individuals $Z_{st}$ with a notification time on day $t$ is drawn from a multinomial distribution:

Note the infection-to-notification time encompasses the incubation period plus any additional delay from symptom onset to notification. This could be explicitly modelled using two separate distributions if required, but we here we consider a single distribution for the overall time from infection to notification.

The total number of individuals with a notification time on day $t$ is $\sum_{s=1}^{t}Z_{st}$. We assume that infections have a fixed probability $\alpha\in(0,1]$ of being notified as cases. Hence, the expected number of case notifications on day $t$ is $\mu=\alpha\sum_{s=1}^{t}Z_{st}$. To allow for noise in daily case notifications, we use a negative binomial distribution for the observed number of case notifications $C_{t}$ on day $t$ [18]:

where $k_{c}$ is the dispersion parameter for the observed data (note that unlike $k_{r}$, $k_{c}$ does not impact the transmission process, only the variance of the observed data). Whilst we retain the negative binomial distribution for model flexibility, in practice we set $k_{c}=\infty$, meaning that the number of daily observed cases is Poisson distributed.

We fit the model to data on daily local case notifications using a bootstrap particle filter [24], using data on imported cases (where available) as seed infections. We simulate a set of $m=10^{5}$ particles according to Eqs. (1)–(4). At each daily time step $t$, the likelihood $w^{(j)}$ for each particle $j=1,\ldots,m$ is calculated using the data $C_{t}^{\mathrm{obs}}$ for the number of case notifications on day $t$. It follows from Eq. (5) that the likelihood is given by

where $\mu^{(j)}$ is the expected number of case notifications on day $t$ for particle $j$, and $f_{NB}$ is the probability mass function for a negative binomial distribution with mean $\mu^{(j)}$ and dispersion parameter $k_{c}$.

Particles are then resampled, with replacement, with probability $w^{(j)}$. After resampling, the particles represent samples from the posterior distribution over the hidden states conditional on the observed data up to day $t$ [25]. To avoid sample degeneracy in the time series of the particles, we used fixed-lag resampling, which means that instead of resampling the entire history of each particle at each time step, we only resample the most recent $t_{\mathrm{lag}}$ days [20].

As in [10], the probability that there are are no new infections on or after day $t$, assuming a constant reproduction number $R$, is given by

where

and $F$ is the cumulative distribution function for the generation time distribution, defined by $F_{s}=\sum_{s^{\prime}=1}^{s}g_{s^{\prime}}$. Note that in Eq. (8), $I_{t-s}$ is replaced by $Y_{t-s}$ to account for the individual heterogeneity in transmission rates and the effect of imported infections.

Because our model explicitly includes delays from infection to notification, it is possible that there could still be cases notified on or after day $t$, even if there are no new infections. To account for this, we also defined the probability $P_{00}$ that there will be no new infections and no notified cases on or after day $t$. As shown in Supplementary Material Section S1, $P_{00}$ may be expressed as

We also calculated the probability of ultimate extinction (i.e. probability that the number of active infections eventually becomes zero, even with an unlimited susceptible population) under constant reproduction number $R$:

where $q\in[0,1]$ is the probability of ultimate extinction for an outbreak that starts with a fully infectious seed case (see Supplementary Material Section S1 for full derivation). PUE is the probability that the transmission chains stemming from any current infections self-extinguish, as opposed to generating a large outbreak that only ends due to depletion of the susceptible population (or reintroduction of control measures). This is equivalent to the “probability of stochastic extinction” defined by [16], but adapted to consider starting with all currently active infections at time $t$. In general, PUE will be larger than $P_{0}$ as there is a possibility that there will be one or more future infections but the outbreak nevertheless self-extinguishes due to stochastic effects.

We calculated the values of $\gamma_{t}$ and $\phi_{t}$ for each of the $m$ particles according to Eqs. (8) and (9) on each day $t$, before resampling the particles according to the likelihood function given the data on the number of case notifications on day $t$. This ensured that the end-of-outbreak probabilities were real-time estimates, i.e. are based only on data that was available prior to day $t$ and not subsequent data.

The pre-intervention reproduction number, denoted $R_{t}^{\mathrm{pre}}$, for each of the $m$ particles was defined to be the average estimated value of $R_{t}$ in the 14 days prior to the start of the intervention. This value was calculated conditional on data up to the start of the intervention but not subsequent data. This provided a sample of $m$ values from the distributional estimate for $R_{t}^{\mathrm{pre}}$.

We defined $P_{0}$, $P_{00}$ and PUE as the probabilities of the respective outcome, under the assumption that the intervention was lifted on day $t$ and the reproduction number consequently reverted to its pre-intervention value. We estimated $P_{0}$, $P_{00}$ and PUE for each particle at each time $t$ by choosing a random sample from the distribution of $R_{t}^{\mathrm{pre}}$ and using this as the value of $R$ in Eqs. (7), (9) and (10). The overall estimates for $P_{0}$, $P_{00}$ and PUE were found by averaging over all $m$ particles. These estimates therefore include the effect of uncertainty in the estimated pre-intervention reproduction number. Note all end-of-outbreak probability estimates assumed that there would be no further imported infections on or after day $t$.

For the New Zealand Covid-19 outbreak, we used Ministry of Health data on the daily number of Covid-19 case notifications in New Zealand between February and June 2020 [26]. We assumed that the generation time (for wildtype SARS-CoV-2) was a gamma distribution with mean 5.05 days and standard deviation 1.94 days [27]. We assumed that the time from infection to notification was a gamma distribution with mean $\pm$ standard deviation of either 7.7 $\pm$ 3.2 days or 11.2 $\pm$ 4.7 days. These correspond to an incubation period of 5.5 $\pm$ 2.3 days, plus values reported by Hendy et al. [28] for the self-reported delay from symptom onset to isolation and from symptom onset to notification, respectively. We use these as indicative values representing a short and a long notification delay scenario respectively, and note that exact values will be outbreak-specific. We truncated the generation time and reporting time distributions at a maximum of 15 days and 25 days respectively, and set $t_{\mathrm{lag}}=30$ days to ensure the resampling lag was greater than both of these maxima [20].

Cases were classified by the Ministry of Health as either local or imported. We assumed that the date of infection for all imported cases was $5$ days prior to date of notification (which means that the first imported case had an infection date of 21 February 2020). We also assumed that imported infections had the same notification probability $\alpha$ as local infections. To model this, we generated a set of unreported imported infections with infection dates resampled with replacement from the infection dates for the notified imported cases. We assumed that imported infections had an average infectivity of $\varepsilon=0.5$ relative to local cases, representing the effect of self-isolation measures for arriving travellers. This value for $\varepsilon$ is consistent with the findings of [29] based on contact tracing data. For simplicity we assumed that imported and local cases had the same generation time distribution. We assumed that any imported cases with an infection date after the introduction of mandatory government-managed quarantine (10 April 2020) had a negligible risk of causing local infections and could therefore be ignored.

Because the case ascertainment rate was highly uncertain, we tested different values for the notification probability $\alpha$ of $0.4$ and $0.7$. We also ran the model for $\alpha=1$ to provide a baseline for assessing the impact of including under-reporting in the model. We also investigated different values for the transmission overdispersion parameter of $k_{r}=\infty$ (representing a Poisson offspring distribution, i.e. no heterogeneity in transmission multipliers), $k_{r}=1$ and $k_{r}=0.2$ (representing increasing levels of heterogeneity). These values span the estimated value of $k_{r}=0.41$ for SARS-CoV-2 from a recent meta-analysis [30].

To initialise the model, we assumed that there were zero infections prior to the first imported infection. We drew the value of $R_{t}$ at the start of the simulation from a prior distribution, assumed to be a Gamma distribution with mean $2$ and standard deviation $1$. We set the time window for the rapid change in $R_{t}$ to be 7 days starting on 23 March 2020, the day on which it was announced that the country would move to “Alert Level 4” two days subsequently. During this $7$-day time window, we set the daily random walk step to have mean $\delta_{t}=-0.1$ and standard deviation $\sigma_{t}=0.2$. This corresponds to a broad prior for the effect of the intervention on $R_{t}$, with an expected aggregate reduction in $R_{t}$ by a factor of $\exp(-7\times 0.1)\approx 0.5$, which is approximately consistent with previous analyses [28]. The chosen values were also found visually to give a good qualitative fit to the data and we provide a sensitivity analysis in Supplementary Material Section S2). All model parameter values are shown in Table 1.

We also analysed data from an EVD outbreak with 54 cases that occurred in Équateur Province of DRC in 2018. Following laboratory confirmation of the first two cases on 8 May, the start of the outbreak was declared and a response team (the Ebola Response Team; ERT) was deployed by the DRC government and international partners [31, 32]. The ERT implemented a range of measures, including ring vaccination, active case finding, contact tracing, case isolation and treatment, laboratory testing, and community engagement. On 24 July, 42 days after the final case recovered, the outbreak was declared over and the ERT was withdrawn [32].

Symptom onset dates were available for all cases so, consistent with previous analyses of this outbreak [10], we analysed the data by date of symptom onset rather than date of notification. This is unlikely to have a substantial effect on model results because, in the latter stages of the outbreak when end-of-outbreak probabilities become relevant, the notification date was the same as the symptom onset date for most cases.

We assumed the generation time was a gamma distribution with mean 15.3 days and standard deviation 9.3 days. This corresponds to the serial interval estimate of [33]. However, since the generation time and serial interval distributions for EVD are thought be very similar, we contend that it is reasonable to use this as the generation time distribution [34]. We tested two different values for the mean and standard deviation of the infection-to-notification time distribution (see Table 1). These correspond to two different estimates for the mean and standard deviation of the incubation period for EVD [35, 36]. We truncated the generation time and reporting time distributions at a maximum of 50 days and 30 days respectively, and set $t_{\mathrm{lag}}=50$ days. We ran the model for three different values for the notification probability $\alpha$ of 0.8, 0.9 and 1 (we considered these relatively high values due to the extensive case finding activities undertaken by the ERT). We also investigated the same three values for the overdispersion parameter $k_{r}$ as for the Covid-19 outbreak. These values span the value of $k_{r}=0.18$ estimated for EVD by Althaus [37]

To initialise the model, we assumed that the first case resulted from a single zoonotic spillover event and all subsequent cases resulted from local human-to-human transmission. We drew the value of $R_{t}$ at the start of the simulation from a prior distribution, assumed to be a Gamma distribution with mean $2.5$ and standard deviation $1$. We set the time window for the rapid change in $R_{t}$ to be 7 days starting on 8 May 2018, the date on which the ERT arrived. During this time window, we set the daily random walk step to have mean $\delta_{t}=-0.1$ and standard deviation $\sigma_{t}=0.2$, as for the Covid-19 outbreak model.

The daily number of locally acquired Covid-19 case notifications in New Zealand peaked at 61 on 4 April 2020 and the last notified case of the outbreak was on 22 May 2020. For the default model parameters (see Table 1), the median reproduction number varied between $1.7$ and $2.2$ in the period prior the start of the intervention, falling to between $0.4$ and $0.5$ after the intervention on 25 March 2020 (Figure 1a). With the assumed case ascertainment rate of $\alpha=0.4$, new local infections peaked at approximately 150–200 per day (Figure 1b) and the model provided a good qualitative fit to the data on local case notifications (Figure 1c).

The estimated probability $P_{0}$ that there would be no future infections if the reproduction number reverted to its pre-intervention value first became non-zero around the beginning of May 2020 (Figure 1c, blue curve). The value of $P_{0}$ tended to increase gradually during intervals with no new case notifications, and to drop sharply when sporadic new cases were notified, for example on 14-15 May and 22 May. After the last case notification on 22 May, $P_{0}$ increased steadily from around 30%, reaching 95% on 6 June 2020 (hereafter, when describing all of our analyses in the text, we refer to relevant probabilities as percentage values). The actual date on which the outbreak was declared to have been eliminated and control measures were relaxed was 8 June 2020. Note a threshold value of 95% for the end-of-outbreak probability is chosen here for illustrative purposes only. Other choices of threshold are equally possible and in reality the appropriate choice of threshold for relaxation of control measures will depend on the epidemiological and socioeconomic context (see Discussion for more detail). See Supplementary Material Section S2 for additional discussion regarding uncertainty in the estimated value of $P_{0}$.

Increasing the assumed notification probability $\alpha$ led to higher estimated values of $P_{0}$ (Figure 2, solid red and yellow curves). For example, the value of $P_{0}$ reached 95% on 3 June 2020 for $\alpha=0.7$ and on 1 June 2020 for $\alpha=1$ (Table 2). This was because a higher assumed case ascertainment rate meant that, for a given number of case notifications, the inferred number of infections was lower and therefore reached zero sooner. However, it was notable that the difference between the estimated time at which $P_{0}$ reached 95% for $\alpha=0.4$ and $\alpha=1$ was only 5 days. Since end-of-outbreak declarations based on this type of analysis would only be expected to be made when the estimated end-of-outbreak probability is close to 1, this result suggests that end-of-outbreak declaration dates may in some situations be relatively insensitive to the assumed case ascertainment rate.

Increasing the mean time $t_{n}$ from infection to notification from $7.7$ days to $11.2$ days meant that the estimated value for $P_{0}$ started to increase earlier and generally was higher at any given time than with a smaller $t_{n}$ (Figure 2, dashed curves). This counterintuitive result occurred because, during the intervention period, the estimated reproduction number was well under the critical value of $1$, meaning that incidence of new infections was on an exponentially decaying trajectory. Observed data at time $t$ gives information about infections that happened around approximately $t-t_{n}$. Hence, the larger $t_{n}$, the longer incidence has been exponentially decaying in the period between infection and observation. Notably however, the time at which $P_{0}$ reached 95% was only $2$–$4$ days earlier with $t_{n}=11.2$ days than with $t_{n}=7.7$ days (Table 2). Additional results showing the hidden states $R_{t}$ and $I_{t}$ for different values of $\alpha$ and $t_{n}$ are provided in Supplementary Figure S1.

We also investigated model results for different levels of individual heterogeneity in transmission and using different definitions for the end-of-outbreak probability. The probability $P_{00}$ of no future infections or notifications (Figure 3a, dotted curves) was lower than the probability of no future infections ($P_{0}$). This is because the delay from infection to notification means that it is possible for cases to be reported on or after day $t$ even if there are no new infections. In contrast, the probability of ultimate extinction (PUE) (Figure 3a, solid curves) was higher than $P_{0}$. This is because transmission chains have a non-zero probability of self-extinguishing, even when the reproduction number is greater than $1$.

Increasing the level of individual heterogeneity in transmission (i.e. reducing the value of the dispersion parameter $k_{r}$) led to higher estimates for the end-of-outbreak probability, whichever of three definitions was used (Figure 3b,c). This was expected as it is well known that increasing individual heterogeneity in the number of secondary infections (decreasing $k_{r}$) increases the likelihood of stochastic extinction [16]. For example, with $k_{r}=0.2$ (Figure 3c), the probability $P_{0}$ of no future infections reached 95% on 30 May 2020, which was 7 days earlier than in the scenario without heterogeneity (Figure 3a). Additional results showing the hidden states $R_{t}$ and $I_{t}$ for different values of $k_{r}$ are provided in Supplementary Figure S2.

A sensitivity analysis on the random walk parameters for the time-varying reproduction number is provided in Supplementary Figure S3 and shows that the results are not highly sensitive to these parameters.

We went on to analyse the 2018 Équateur Province EVD outbreak dataset. The highest number of daily cases was six on 4 May 2018, just before the arrival of the ERT on 8 May. For the default model parameters (see Table 1), the median reproduction number varied between $2.0$ and $2.9$ in the period prior to the start of the intervention, falling to approximately $0.3$–$0.4$ after the intervention (Figure 4a).

Because of the longer generation time of EVD compared to Covid-19, the estimated probability $P_{0}$ of no future infections did not start to increase substantially above zero until around 2 weeks after the last reported case on 2 June 2018, and subsequently increased smoothly towards $1$ (Figure 4c). This contrasted with the results for Covid-19 where the estimated end-of-outbreak probabilities started to increase after a few days with no newly reported cases but dropped down after a new case notification, leading to sawtooth-shaped curves.

The effects of the notification probability $\alpha$ and mean time to notification $t_{n}$ on the estimated probability $P_{0}$ of no future infections were broadly similar for EVD as for Covid-19. Increasing either $\alpha$ or $t_{n}$ tended to increase $P_{0}$ (Figure 5). For the parameter values investigated, the estimated probability $P_{0}$ of no future infections reached 95% in the range 22 July to 30 July 2018 (Table 2). The actual date on which the ERT was withdrawn was 24 July 2018. Additional results showing the hidden states $R_{t}$ and $I_{t}$ for different values of $\alpha$ and $t_{n}$ are provided in Supplementary Figure S4.

The relative behaviour of the alternative definitions of the end-of-outbreak probability and the effect of heterogeneous transmission patterns were also similar for EVD as for Covid-19. Reducing the dispersion parameter $k_{r}$ increased the estimated end-of-outbreak probabilities (Figure 6). This was particularly true for the PUE (Figure 6, solid curves), which, for the most highly overdispersed scenario (Figure 6c) was approximately $0.1$ even during the middle of the outbreak. This reflects the fact that small outbreaks with highly overdispersed transmission patterns have a significant probability of stochastic extinction, even if the reproduction number is well above $1$.

The other two end-of-outbreak metrics, $P_{0}$ and $P_{00}$ (Figure 6, dashed and dotted curves respectively), were very similar to each other for EVD. This was a consequence of the fact that the generation time was long relative to the time from infection to notification, meaning the likelihood of there being cases yet to be notified but no longer infectious was small. Additional results showing the hidden states $R_{t}$ and $I_{t}$ for different values of $k_{r}$ are provided in Supplementary Figure S5. A sensitivity analysis on the random walk parameters for the time-varying reproduction number is provided in Supplementary Figure S6 and shows that the results are not highly sensitive to these parameters.

The authors would like to thank the Isaac Newton Institute for Mathematical Sciences, Cambridge, for support and hospitality during the programme “Modelling and inference for pandemic preparedness”, where initial work on this research article was undertaken. This research was supported by EPSRC grant EP/R014604/1 (MJP, WSH and RNT). MJP acknowledges funding from Te Niwha Infectious Diseases Research Platform, Institute of Environmental Science and Research, grant number TN/P/24/UoC/MP. WSH and RNT would like to thank members of the Infectious Disease Modelling group in the Wolfson Centre for Mathematical Biology at the University of Oxford for useful discussions about this research. The authors are grateful to two anonymous reviewers for helpful comments.