Hindsight is 2020 vision: a characterisation of the global response to the COVID-19 pandemic

Abstract

Background

The global impact of COVID-19 and the country-specific responses to the pandemic provide an unparalleled opportunity to learn about different patterns of the outbreak and interventions. We model the global pattern of reported COVID-19 cases during the primary response period, with the aim of learning from the past to prepare for the future.

Methods

Using Bayesian methods, we analyse the response to the COVID-19 outbreak for 158 countries for the period 22 January to 9 June 2020. This encompasses the period in which many countries imposed a variety of response measures and initial relaxation strategies. Instead of modelling specific intervention types and timings for each country explicitly, we adopt a stochastic epidemiological model including a feedback mechanism on virus transmission to capture complex nonlinear dynamics arising from continuous changes in community behaviour in response to rising case numbers. We analyse the overall effect of interventions and community responses across diverse regions. This approach mitigates explicit consideration of issues such as period of infectivity and public adherence to government restrictions.

Results

Countries with the largest cumulative case tallies are characterised by a delayed response, whereas countries that avoid substantial community transmission during the period of study responded quickly. Countries that recovered rapidly also have a higher case identification rate and small numbers of undocumented community transmission at the early stages of the outbreak. We also demonstrate that uncertainty in numbers of undocumented infections dramatically impacts the risk of multiple waves. Our approach is also effective at pre-empting potential flare-ups.

Conclusions

We demonstrate the utility of modelling to interpret community behaviour in the early epidemic stages. Two lessons learnt that are important for the future are: i) countries that imposed strict containment measures early in the epidemic fared better with respect to numbers of reported cases; and ii) broader testing is required early in the epidemic to understand the magnitude of undocumented infections and recover rapidly. We conclude that clear patterns of containment are essential prior to relaxation of restrictions and show that modelling can provide insights to this end.

Background

A unique feature of the coronavirus disease 2019 (COVID-19) pandemic has been the rapid and widespread availability of data through online platforms [1–3]. These data enable the analysis of various patterns of outbreak containment, and provide an unparalleled opportunity to learn how to respond to a new surge of COVID-19 or future pandemics. It is important to learn from the past to prepare for the future.

Travel restrictions, increased hygiene education, social distancing, school and business closures, and complete lockdowns [4–6] are examples of non-pharmaceutical intervention (NPI) strategies that many countries have introduced to slow transmission rates and relieve pressure on healthcare systems in the absence of a vaccine or treatment for COVID-19 [7]. Modelling is at the forefront of determining the efficacy of these measures in reducing severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) transmission and quantifying risk of future outbreaks along with their potential severity [8–11]. This understanding is crucial given the deleterious sociological and economic impacts of many NPIs [4, 12–14].

The global modelling community has provided insight into the transmissibility of SARS-CoV-2 [15–17], global risk of spread through transport networks [18, 19], forecasting and prediction [20–22], and evaluation of interventions [10, 23]. Techniques include: empirical approaches such as phenomenological growth curves [22]; data-driven, statistical approaches using non-linear autoregressive models [24]; and mechanistic models based on epidemiological theory [25] with various extensions [26, 27].

We aim to learn about the global pattern of behaviour among countries based on the trajectories of reported cases, recoveries and fatalities as provided by Johns Hopkins University (JHU) [18, 28]. Notwithstanding acknowledged drawbacks in relying on reported cases [2], we argue that such data will be the main source of information for government and health managers in future scenarios. We avoid imposing the specific complex history of intervention measures for each country by including a novel regulatory mechanism that captures the changes in community behaviour in response to rising confirmed cases. Through model calibration, we infer the country specific response timing and strength. As a result, we analyse the overall effect of interventions and community responses across diverse regions.

We characterise the response of 158 countries to the COVID-19 outbreak for the period 22 January to 9 June 2020. This time frame encompasses the period during which initial measures were imposed by most countries and the period in which some countries started to relax restrictions. Our analysis is a broad assessment of the global response to the COVID-19 pandemic, and reflects how countries in early phases of outbreak may have adjusted their strategies to reduce the time to recovery. We find that very large outbreaks are characterised by a delayed response, whereas countries that observed a decrease in active cases during the early period of study are characterised by high case identification rates. For many countries, the transmission rates were declining in the later period of the study. However, large unobserved infected population counts were estimated. Our analysis confirms that multifaceted approaches that include NPIs, increased testing, contact tracing, isolation and quarantine measures are effective in reducing the severity of COVID-19 outbreaks world-wide. We also demonstrate, that the magnitude of undocumented cases substantially impacts uncertainty in risks of subsequent flare-ups after restrictions are relaxed. We conclude that wider testing is essential to reduce this uncertainty to reliably evaluate risk of future waves.

Methods

Ethics

Ethics approval was not required for this study because data are publicly available from the JHU coronavirus data repository [1] (https://github.com/CSSEGISandData/COVID-19). These data are daily reported confirmed case, recovery, and death numbers that do not contain any confidential or identifiable patient data. Theses publicly available data have been widely utilised by governments, news outlets and other epidemiological studies.

Data summary

Daily counts of reported confirmed COVID-19 cases, recoveries and deaths for each country are obtained from the JHU coronavirus resource center [1, 18] (publicly available at https://github.com/CSSEGISandData/COVID-19). We refer the reader to the “Discussion” section for comments on this data source. Population data for 2020 were obtained from United Nations Population Division estimates [29].

We analyse three different time periods: i) 22 January–30 March; and ii) 22 January–13 April; and iii) 22 January–9 June. These periods are selected as they broadly represent the time period of the initial outbreak of COVID-19; covering the initial exponential growth period, the first epidemic peak, and subsequent initial recovery period of many countries. We use these time points to look at the changes in key model parameters relating to a countries responses over time.

Countries are included in the analysis for a give time period provided the cumulative number of confirmed COVID-19 cases exceeded 100 at least one day prior to the end of the particular analysis period. While the specific threshold is empirically chosen as 10², the lower bound on the initial case numbers is needed to ensure that the infectious population is large enough for sufficient mixing to occur so that a compartmental model is a reasonable approximation (i.e., as threshold of 10¹ initial cases will lead to poor approximations and 10³ cases will miss the early stages of the pandemic in many countries). The initial case number threshold should be exceeded at least one day before the end of the analysis period to ensure there are at least two observations in the time series for all countries. Using these inclusion criteria, we obtain N=98 countries for the period of 22 January–30 March, N=121 countries for period of 22 January to 13 April, and N=158 countries for the period 22 January–9 June. Countries included in one period are not removed from subsequent periods and lower values of N in earlier periods reflect the fact that fewer countries had experienced outbreaks by that time.

Analysis summary

For each country, i=1,2,…,N, the JHU maintains a time-series, ${\mathcal {D}}_{i} = \left [\left \{C_{t,i},R_{t,i},D_{t,i}\right \}_{T \geq t \geq 0}\right ]$, where C_t,i,R_t,i, and D_t,i are, respectively, the cumulative confirmed cases, case recoveries and case deaths on day t for country i, t=0 is the first day such that C_t,i≥100 and t=T is the end of the study period. Since there are variations in reporting protocols across countries and time as well as data curation challenges [2], caution is necessary in the interpretation of our analysis across all countries over time.

Bayesian parameter inference is applied over three time periods. The first period, 22 January to 30 March, is used to assess the community response to the initial outbreak of COVID-19. The second period, including data up to 13 April, encompasses the time period in which the efficacy of the community response starts to become evident. Finally, the third period includes data up to the 9 June, in which many countries had started to relax restrictions. We also consider in this analysis the prevention of future waves, and highlight the sensitivity of system dynamics to the uncertainty in unobserved infectious individuals.

Mathematical model

A stochastic epidemiological compartmental model is used to describe the spread of COVID-19 within a single country over the time period t∈(0,T]. The assumed well-mixed population of size P is comprised of six compartments: susceptible, S_t; infectious, I_t; confirmed active cases, A_t; case recoveries, R_t; case fatalities, D_t; and unconfirmed recoveries, $R_{t}^{u}$. The population that is susceptible to the SARS-CoV-2 infection, S_t, can be infected by individuals from the unobserved infectious population, I_t, including both symptomatic and asymptomatic infections. The active confirmed cases, A_t, are those who have tested positive for COVID-19 but have not yet recovered or died. We assume individuals in A_t are isolated from the susceptible community (e.g., self-isolated, quarantined, or hospitalised) and no longer contribute to new infections. Importantly, A_t need not be symptomatic, but may have been identified from contact tracing protocols or community wide testing. R_t and D_t are, respectively, the population of confirmed cases that recover or die. These correspond to the recoveries and deaths reported in the JHU data. Lastly, $R_{t}^{u}$ is the population of infected individuals that recover or die without being tested for COVID-19; these individuals no longer spread the infection but do not contribute to the reported recovery and fatality counts. The cumulative confirmed cases, as reported by the JHU, can be obtained by C_t=A_t+R_t+D_t.

The populations S_t,I_t, and $R_{t}^{u}$ are not observable and are latent variables in our model. Therefore, strategies for managing the spread of the virus, such as NPIs, are informed by the observables, A_t,R_t, and D_t. Media coverage, official government information, and health authority reports based on these observables may also affect the behaviour of individuals. For example, frequent reports on growing case numbers may increase voluntary self-isolation; conversely, media coverage that downplays the risk of infection or seriousness of the disease may lead to widespread non-compliance with health advice or government regulations. We model this dynamic introducing a feedback loop within the transmission process.

A schematic of this system that highlights the state transitions and the feedback loop is given in Fig. 1. The dynamics can be described by the differential equations,

$$ \begin{aligned} {\dot{S}_{t}} &= -g\left(A_{t},R_{t},D_{t}\right)S_{t} I_{t}/P, \quad {\dot{I}_{t}} = -(\gamma + \eta\beta)I_{t} + g(A_{t},R_{t},D_{t})S_{t} I_{t}/P, \\ {\dot{R}^{u}_{t}} &= \eta \beta I_{t}, \quad {\dot{A}_{t}} = \gamma I_{t} - (\beta + \delta)A_{t}, \quad {\dot{R}_{t}} = \beta A_{t}, \quad\text{and}\quad {\dot{D}_{t}} = \delta A_{t}. \\ \end{aligned} $$

(1)

Fig. 1

Schematic of epidemic model including a regulatory mechanism inducing a feedback loop. State transitions are marked by arrows with superscripts indicating respective rate parameters. Here, observable quantities can inform individual behaviour to inhibit transmission in the latent SIR model

Here, g(·)>0 is the transmission rate function, γ>0 is the identification rate, β>0 is the case recovery rate, δ>0 is the case fatality rate, and η>0 is the latent removal rate relative to the case recovery rate. Initial conditions for the observables, A₀,R₀,D₀, are obtained from the JHU data. To capture uncertainty in community spread at early time we set I₀=κA₀, where κ>0 is the relative number of unobserved cases. We assume initially $R_{0}^{u} = 0$ and S₀=P−C₀−I₀. Although Eq. (1) shows a deterministic system for ease of interpretation, we apply a stochastic equivalent for the analysis in this work (See Supplementary Material).

The novel feedback mechanism provides a general framework to describe how communities change their behaviour as case numbers rise. This is similar to the influence of media reports that have been the subject of study for influenza and HIV [30, 31]. However, our approach includes a response strength parameter.

We define a so-called reporting function,

$$ U(A_{t},R_{t},D_{t}) = w_{A} A_{t} + w_{R} R_{t} + w_{D} D_{t}, $$

(2)

where the weights w_A,w_R,w_D≥0, represent the relative weighting of observables in contributing to information that influences individual behaviour, introduction of NPIs, and subsequent compliance with government regulation or health advice. In the context of this work, the weights w_A,w_R, and w_D have a very important interpretation, but we first need to present more details of the feedback mechanism.

We consider a nonlinear transmission rate of the form,

$$ g(A_{t},R_{t},D_{t}) = \alpha_{0} + \alpha f(U(A_{t},R_{t},D_{t})), $$

(3)

where the response function, f(·)∈[0,1], is a decreasing function with respect to U(·),α is the controllable transmission rate such that αf(·) is a transmission rate that decreases as the reporting function increases and α₀ is the residual transmission rate as f(·)→0. The strength of the response s=(1−f(·))×100% is the percentage reduction in community transmission, excluding residual transmission α₀.

For the response function we assume the form

$$ f(U(\cdot)) = \frac{1}{1+U(\cdot)^{n}}, $$

(4)

where the parameter n≥0 controls the rate of decrease with respect to the reporting function. This form is selected for two reasons. Firstly, it generalises techniques that capture the influence of media reports during epidemics [30]. Secondly, the weights from Eq. (2) have an important interpretation. This can been seen by noting that values for A_t,R_t and D_t that satisfy the condition U(A_t,R_t,D_t)=1, indicate the threshold case numbers that leads to a response strength of 50%, that is, f(U(·))=1/2 leading to g(·)=α₀+α/2. The effect of the slope parameter, n, and the weights are shown in Fig. 2. If U(·)=0, that is no cases are reported, or w_A=w_R=w_D=0, indicating no perceived risk, then the model reduces to an SIR model in the unobserved population with transmission rate α₀+α.

Fig. 2

Effect of parameters on the response function. a The effect of the slope parameter n. Note, as n increases, the faster f(·)→0. For any n we have f(·)=1/2 (dashed black), at the point U(·)=1 (dotted black). b The effect of weights on the response function for the special case w_A=w_R=w_D=w>0 for constant n=5. Note the point at which f(·)=1/2 corresponds to A_t+R_t+D_t=1/w. That is, as w increases the lower the number of cases are required to influence the community to reduce the spread

For n≤1 the shape of f(U(·)) starts to decline rapidly levelling out (Fig. 2a). Increasing n>1 results in a decreasing sigmoid curve with an inflection point at U(·)=1 in which a population response strength reaches 50%. Small n describes a population that does not significantly reduce the transmission rate until the U(·) is large. Conversely, larger n describes a population that acts decisively as a response that rapidly reduces transmission around U(·)=1. Large values of weights w_A,w_R,w_D correspond to lower acceptable thresholds of cases, including active, recovery and death counts. Lower weights lead to delayed responses. That is, the parameter n relates to the rate of intervention introduction and the weights relate to decision thresholds and subsequent compliance. Importantly, our approach does not distinguish between different NPIs and voluntary population behaviour, but rather models the net effect that reporting has on transmission rates.

We focus on the reporting function with w_R=w_D=0, that is, U(A_t,R_t,D_t)=w_AA_t (see “Discussion” section and Supplementary Material for alternatives). Since A_t will increase and decline over to course of an outbreak, the model can exhibit oscillations that are essential for understanding the potential for flare-ups and multiple waves.

Bayesian analysis

Parameter inference

Our model has up to 11 parameters: two transmission rate parameters, α₀ and α; case recovery rate β; case identification rate γ; case death rate δ; relative latent recovery rate η; response slope parameter n; the initial infected scale factor κ; and the weights of the reporting function w_A,w_R, and w_D. We assume w_R=w_D=0 and infer the parameters θ=[α₀,α,β,γ,δ,η,n,κ,w_A] (see SupplementaryMaterial for sensitivity analysis for the general case).

Using the daily case data, ${\mathcal {D}}_{i}$, for each country i∈[1,2,…,N], we infer model parameters within a Bayesian framework by sampling the joint posterior distribution,

$$ {p({{\theta}}\mid{{\mathcal{D}}_{i}})} \propto {p({{\mathcal{D}}_{i}}\mid{{\theta}})} {p({\theta})}, $$

(5)

where ${p({{\mathcal {D}}_{i}}\mid {\theta })}$ is the likelihood and p(θ) is the prior distribution.

We rely on adaptive sequential Monte Carlo for approximate Bayesian computation (SMC-ABC) [32–35] to obtain approximate posterior samples since the likelihood is intractable (Supplementary Material). We use independent uniform priors, $\alpha _{0},\alpha,\beta,\gamma,\delta,\eta \sim \mathcal {U}(0,1), n \sim \mathcal {U}(0,20), \kappa \sim \mathcal {U}(0,100)$, and $\log _{10} w_{A} \sim \mathcal {U}({-6},{-2})$.

Assessment of model fit and prediction

The highly variable nature of the COVID-19 pandemic makes it notoriously difficult to predict [20, 21]. Our purpose is not to provide forecasts, but rather to capture the dynamic effects of changes in community behaviour during the outbreak. As a result, our model needs to be able to capture the overall trends in daily cases.

Model fit is assessed through sampling the posterior predictive distribution

$$ {p({{\mathcal{D}_{s}}}\mid{{\mathcal{D}}_{i}})} = \int {p({{\mathcal{D}_{s}}}\mid{{\theta}})}{p({{\theta}}\mid{{\mathcal{D}}_{i}})}\, \mathrm{d}{\theta}, $$

(6)

where ${\mathcal {D}_{s}}$ is simulated data as generated by the model. We compute the 50% and 95% credible interval (CrI) of ${p({{\mathcal {D}_{s}}}\mid {{\mathcal {D}}})}$ by: 1) generating simulated data for each posterior sample generated by the SMC-ABC sampler; and 2) computing quantiles of the simulation state distribution at each observation time.

The posterior predictive distribution is also used to asses flare-up risk. After model calibration using data up to 9 June, we continue simulations to 24 June and obtain CrIs for oscillatory behaviour relating to localised flare-ups and additional waves.

Parameter point estimation and uncertainty quantification

Parameter point estimates are also obtained from the approximate posterior sample with the lowest average discrepancy with the observed data (See SupplementaryMaterial). Parameter uncertainty is reported using 95% CrI for the marginal approximate Bayesian posterior distributions (See Supplementary Material). This uncertainty quantification encompasses all plausible parameter combinations within the achieved acceptance threshold of the ABC-SMC method.

Correlation analysis

We use the point estimates to evaluate which factors have had the greatest impact on the COVID-19 outbreak evolution across countries. For each country, we manually classify the state of the outbreak for each country on 13 April based on the trend in the daily reported cases and active case numbers. These stages are: the growth stage–characterised as an increasing trend in daily reported cases numbers; the post-peak stage–characterised by declines in daily case numbers, indicating the curve is flattening; the recovery stage–characterised by declines in active case numbers. Spearman’s rank-order correlation coefficients are computed between each parameters and observed data at T= 13 April (i.e., cumulative case numbers, C_T, recoveries, R, and deaths, D_T).

Results

Assessment of model fit

For most countries, the 95% CrI contains the daily case, recovery and death data (See Supplementary Material), and the 50% CrI overlaps the main trend in the data. Exceptions to this are largely consistent with reporting delay effects, such as weekly seasonality as evident in the daily cases for Germany (Figure S13). For some countries, particularly those that responded rapidly (e.g., Australia Figure S10m), the early portion of the time series sits in the tails of the 95% CrI and does not overlap with the 50% CrI. In a few cases, sudden spikes in daily numbers (e.g., Recoveries in Germany Figure S13h) are not completely captured within the 95% CrI, however, the model does match overall trend well with many daily numbers remaining within the 95% CrI. For example, the possible decline in daily cases numbers for the United States (Figure S19g) is captured by the lower bound of the 95% CrI, however, the uncertainty of this trend on 13 April is indicated by the increasing upper bound. In some extreme cases, such as changes in reporting methodology from China on 13 February [36] (Figure S12), subsequent inaccuracies occur (Figure S12). We discuss potential model improvements to account for this in the “Discussion” section. Notwithstanding this, our model appears to capture the overall trends in the trajectories to facilitate a broad comparative analysis of global responses.

Example model outputs using the parameter point estimate, based on the lowest expected data discrepancy (See “Methods” section), are provided in Fig. 3 (See Sup-plementary Material for other examples). In particular, compare Fig. 3a–b with Figure S19g–i. Despite further increases in daily cases being highly plausible on the 13 April, the point estimates appropriately exclude these trajectories from the response analysis.

Fig. 3

Examples of model fit using parameter point estimates: a–b United States, c–d Germany, e–f Australia, g–h United Kingdom, i–j South Korea, and k–l New Zealand. Vertical bars indicate daily reported cases (yellow) and deaths (red). The 50% (dark shaded region) and 95% credible intervals (light shaded region) of the posterior predictive distributions are plotted against the observational data. Credible intervals are computed using n=100 stochastic simulations for the given point estimate. Full posterior predictive distributions are presented in the Supplementary Material

Characterisation of responses

Based on our correlation analysis (See “Methods” section), these key parameters are: the case identification rate, γ, the relative initial undocumented infections, κ, and the response weight w_A. Therefore, we evaluate our parameter inferences and point estimates for the three analysis periods. These results reveal that as the pandemic evolves, we can learn more about the possible response histories and latent infected populations. Figure 4 shows pairwise scatter plots for the key parameters w_A,γ, and κ for each of the three analysis periods.

Fig. 4

Pairwise scatter plots of point estimates of each assessed country (grey points) for the key parameters related to the management of an COVID-19 outbreak up to: a–c 30 March; d–f 13 April; and g–i 9 June. a,d,g w_A versus γ; b,e,h w_A versus γ; and c,f,i κ versus γ. For each time period, countries with the ten largest confirmed case counts are highlighted (red points) along with representative countries that were recovering or managed to control the outbreak (green). Labels identify the country by ISO-3166 alpha-3 code

The interactions between parameters are complex, and the differences between time periods deserves some interpretation. However, we first highlight some overarching trends across all time periods and then discuss specific details. One clear trend is countries with that largest numbers of cumulative cases (Fig. 4, top ten countries for cumulative cases numbers indicated in red) tend to have lower response weights, typically w_A≈10⁻⁴. Some more extreme cases are as low as w_A≤10⁻⁵, for example, the United States (USA; Fig 4d,g), Russia (RUS; Fig 4g) and Brazil (BRA; Fig 4g). Small w_A indicates a delayed response in which the transmission rate did not decline significantly until active cases, A_t, increased to larger numbers. This is consistent with reported delays in response across Europe and the United States [4]. Low case identification rates, γ<0.01, and higher relative initial unobserved infections κ>10 are also characteristics of countries will large cumulative confirmed case counts. Countries that controlled the outbreak during the period of study, such as Australia (AUS), New Zealand (NZL), South Korea (KOR), and Taiwan (TWN) are characterised by either rapid responses, w_A>10⁻³, lower relative initial cases numbers, κ, or higher case identification rates γ.

On 30 March (Fig. 4a–f), the top ten countries having the largest cumulative cases (highlighted red) are: the United States (USA), Italy (ITA), Spain (ESP), China (CHN), Germany (DEU), France (FRA); Iran (IRN); the United Kingdom (GBR); Switzerland (CHE); and Belgium (BEL). Of these, only China was recovering (highlighted green). We also highlight South Korea (KOR) as the only other country in recovery during this period, and Taiwan (TWN) as substantial community outbreak was avoided altogether. China and South Korea are characterised by a higher identification rate (China γ=0.96; 95% CrI [0.01,0.90]; South Korea γ=0.28; 95% CrI [0.17,0.93]) and lower relative initial undocumented cases (China κ=0.13; 95% CrI [0.08,35.2]; South Korea κ=2.56; 95% CrI [0.50,5.93]) which is indicative of their strict testing, isolation, and tracing regimes [37, 38]. Taiwan (TWN), with a high response weight (w_A=10^−1.9; 95% CrI [10^−5.9,10^−1.1]) responded very rapidly, having established strong public health response mechanisms after the 2003 severe acute respiratory syndrome (SARS) outbreak. This is also reflected in a low level of initial undocumented cases κ=6.48 (95% CrI [0.18,21.95]) for Taiwan [39]. This is in stark contrast to Iran (IRN; κ=33.79; 95% CrI [4.7,91.41]) that experienced substantial community transmission ahead of the first reported cases. The large response weight for Iran is results in almost no effective reduction in community transmission since α₀ is much larger than α (See Supplementary Material), which could reflect the effects of large gatherings [7, 40, 41].

By 13 April, the situation changes as COVID-19 spreads to more countries and response strategies are altered. The top ten countries having the largest cumulative case numbers had changed to: The United States (USA); Spain (ESP); Italy (ITA); France (FRA); Germany (DEU); The United Kingdom (GBR); China (CHN); Iran (IRN); Turkey (TUR); and Belgium (BEL). By this time, Australia (AUS) and New Zealand (NZL) were included in the ranks of countries that were starting to recover. In Fig. 3d–e, there is a substantial decrease in the response weight for many of the worst affected countries (USA, w_A=10^−5.2 95% CrI [10^−5.7,10^−3.2]; FRA, w_A=10^−4.7 95% CrI [ 10^−5.4,10^−1.3; ESP, w_A=10^−4.4 95% CrI [ 10^−4.7,10^−3.9; DEU, w_A=10^−4.4 95% CrI [ 10^−4.8,10^−4.1; IRN, w_A=10^−4.4 95% CrI [ 10^−5.2,10^−1.8; ITA w_A=10^−4.2 95% CrI [ 10^−5.8,10^−1.2; GBR w_A=10^−4.0 95% CrI [ 10^−5.1,10^−1.1;). This indicates that, in light of data between 31 March–13 April, the community response to the outbreak was even more delayed than earlier data indicated. The very high values of w_A for New Zealand (w_A=10^−2.56; 95% CrI [10^−5.43,10^−1.25]) demonstrates that a rapid response has been a key factor in keeping cumulative case number low (C_T<1,500). Many countries in the top ten cumulative case numbers still had a larger point estimates for κ>10 (Fig. 4e,f), with the United Kingdom having the largest point estimate of κ=42.73 (95% CrI [0.46,91.52]), which could be the result of early unobserved transmission prior to abandonment of “herd immunity” targets in favour of social distancing, and closing of non-essential business and schools [10, 42]. Australia has a very similar characterisation to the United Kingdom in terms of response weight (w_A=10^−3.36; 95% CrI [10^−3.59,10^−3.22]) but with much lower κ=2.72 (95% CrI [0.24,2.78]) and higher γ=0.16 (95% CrI [0.16.0.99]) that likely reflects fact that many of confirmed cases within Australia during this period were imported cases and local community transmission was low [17]. Small increases in γ for Germany (γ=0.15; 95% CrI [0.01,0.87]), Italy (γ=0.10; 95% CrI [0.0,0.92]) and France (γ=0.07; 95% CrI [0.01,0.84]), and a large increase in γ for the United States (γ=0.21; 95% CrI [0.0,0.89]); this is possibly a reflection of increased testing capabilities within these countries between 31 March to 13 April [43]. There is also decrease in κ overall, indicating that the number of early undocumented infections could be lower than previously thought. Especially for Germany with κ=3.27 (95% CrI [0.44,38.63]).

In the period up to 9 June (Fig. 4g–h), Brazil (BRA), Peru (PER), Russia (RUS) and India (IND) have replaced Belgium (BEL), China (CHN), Iran (IRN) and Turkey (TUR) in top ten countries for cumulative case numbers as the epicentre of the COVID-19 pandemic shifts away from Europe. These four new countries in the top ten cases list are characterised by low response weights (BRA, w_A=10^−5.2 95% CrI [10^−5.7,10^−1.1]; IND, w_A=10^−4.1 95% CrI [10^−5.8,10^−1.0]; PER, w_A=10^−3.8 95% CrI [10^−5.6,10^−1.1]; RUS, w_A=10^−4.9 95% CrI [ 10^−5.6,10^−4.7]), low identification rates (BRA, γ=0.03 95% CrI [0.01,0.90]; IND, γ=0.01 95% CrI [0.01,0.88]; PER, γ=0.004 95% CrI [0.002,0.83]; RUS, γ=0.05 95% CrI [0.02,0.92]) and high relative numbers of initial undocumented infections (BRA, γ=79.34 95% CrI [1.19,80.59]; IND, γ=47.62 95% CrI [0.95,92.97]; PER, κ=15.47 95% CrI [0.80,89.93]; RUS, κ=19.39 95% CrI [0.44,61.41]). Many of the mainland European countries were recovering and relaxing restrictions imposed by intervention strategies, consequently the parameters for Germany (DEU), Spain (ESP), France (FRA), and Italy (ITA) have not changed much from the previous analysis other than a further reduction in the estimated relative initial undocumented case numbers for Germany (κ=0.11; 95% CrI [0.01,43.80]). While the United Kingdom (GBR) now has higher response weight, this is unfortunately offset by a larger initial relative number of undocumented infections. For the United States (USA), there is a decline in the identification rate (γ=0.14 95% CrI [0.01,0.96]).

Overall assessment of the global response to the COVID-19 outbreak

Our correlation analysis (See “Methods” section) provides an overview of the relationship between model parameter estimates and the magnitude of the COVID-19 outbreak. Here we use the point estimates computed for 121 countries for time period 22 January to 13 April. Figure 5 highlights the output from this analysis. The lower diagonal section of Fig. 5 show the distribution of point estimates and outbreak stage classifications, whereas the upper diagonal show the correlation coefficients.

Fig. 5

Distributions of model parameter point estimates along with observered cumulative confirmed cases C_T, recoveries R_T and deaths D_T at T= 13 April. Pairwise scatter plots on the lower diagonal indicate the stage of the COVID-19 outbreak for that country: growth stage (red circles), post-peak stage (purple triangles), or recovery stage (green squares). Histograms on the diagonal show the distribution of parameters across all countries within each outbreak stage. Spearman correlation coefficients between each point estimate and observed case numbers with the sign and strength of the correlation indicated by the colour-map (positive correlations in red and negative correlations in blue)

The two parameters with the strongest correlation with large case numbers, C_T, are the response weight w_A (Spearman’s ρ=−0.4357), and the relative initial undocumented cases numbers, κ, (Spearman’s ρ=0.4439). These parameters relate to the response behaviour of communities at the early stages of the outbreak. Recall A_t=1/w_A is the critical number of active cases to invoke a response strength of 50%, corresponding to a reduction in transmission rate of $\frac {\alpha _{0} +\alpha /2}{\alpha _{0} + \alpha } \times 100\%$. A smaller value of w_A indicates a delay in community response, since larger numbers of active cases are required to invoke a response strength of 50%. Large values for κ are indicative of community transmission occurring ahead of the earliest reported cases. There is a strong negative correlation between the identification rate γ and κ (Spearman’s ρ=−0.7896), indicating that countries with strong testing and contact tracing regimes were able to minimise the amount of undocumented community transmission at early time. The identification rate also has low negative correlation (Spearman’s ρ=−0.2284) with the residual transmission rate, α₀, meaning that countries with stronger testing regimes also improved maximum efficacy of other interventions. The response slope n had only weak correlations with any other parameters. This parameter relates to the rate of change in community behaviour before and after the critical A_t=1/w_A point (See Fig. 2). The weak positive correlation between κ and n could mean that countries adopting a more gradual introduction of interventions (lower n) also tended to have less undocumented community transmission initially.

While the three outbreak classification stages (See “Methods” section) are not well separated in all parameters, there are a few trends to note. Countries in the recovery stage tend to have lower residual transmission, α₀, larger regulated transmission, α, larger case recovery and identification rates, β and γ, lower death rates δ, and lower relative undocumented initial infections κ. For countries in the growth stage of the outbreak, the converse is true. The post-peak stage have, on average, parameter values that sit between the recovery and growth stages. Some countries that experienced large numbers of cases in this analysis period are also in the post-peak or recovery stage, whereas others in the growth stage had only small numbers of cases at this time. In some case, such as countries of South America, this analysis could have helped highlight the importance of interventions early in the pandemic.

Avoiding multiple waves

Reliably estimating the number of undocumented infections, including the extent of asymptomatic by infectious case, is crucial to avoid flare-ups that potentially can lead to multiple waves of COVID-19 spread [44]. Recent evidence suggests asymptomatic individuals having a substantial role in the spread of COVID-19 [45]. Quantification of uncertainty in the unobserved infectious population is crucial for planning the timing of easing of restrictions.

Due to the form of our response function (See “Methods” section), we model both how communities introduce and subsequently relax interventions in response to active case numbers. As A_t declines then f(·) increases to simulate increased mixing of the population. The posterior predictive distribution can demonstrate oscillatory behaviours, the magnitude of which depends on the evolution of the undocumented infectious population I_t. Figure 6 demonstrates this. Uncertainty in the daily case numbers, driven by parameter uncertainty and undocumented infectious population uncertainty, pre-empts the small flare-ups in numbers. To obtain this behaviour, all possible evolutions of I_t that are consistent with the observed daily cases must be taken into account.

Fig. 6

Example of small secondary oscillations in model behaviour using the model fit against daily case data (yellow bars) for a Australia and b South Korea up to 9 June (dashed line). The posterior predictive simulations are continues up to 24 June to demonstrate the uncertainty in potential case increases after relaxation of restrictions. Actual daily case numbers for the period 10–24 June (red bars) also demonstrate increases within the credible intervals (dark blue 50% CrI; light blue 95% CrI)

This highlights the importance of conservative uncertainty estimation for forecasting recovery from the pandemic. While it looks encouraging to see predictions that indicate potential future declines in active cases counts, the uncertainties in the number of undocumented cases and model parameters indicate caution in predicting the timing of consistent declines in active case numbers. Therefore, it is essential that communities remain vigilant in fast-evolving situations such as this.

Discussion

We have applied a novel stochastic epidemiological model to characterise the response to the first wave of the COVID-19 pandemic. We find that the worst affected countries (in terms of confirmed case numbers), are characterised by a delayed response (small w_A), allowing case numbers to rise before interventions became effective. However, increased testing and isolation protocols (large γ) have demonstrably reduced the longer term impact, as demonstrated by China and South Korea. Many countries seem to be learning from these collective experiences, with more rapid responses (large w_A). Unfortunately, we also identify that the number of undocumented cases likely substantially exceeded the confirmed cases for many countries (large κ). It is important to emphasise that we do not make any specific guidelines for particular countries in improving specific non-pharmaceutical intervention (NPI) or testing strategies. However, in light of our analysis, we advise that intervention mechanisms be mobilised rapidly without waiting for large numbers of cases to be confirmed. This has been a key characteristic of countries that have successfully managed the initial outbreak, such as Australia, New Zealand and Taiwan.

The data on reported daily cases from the Johns Hopkins University coronavirus resource center have some limitations [2]. Firstly, when aggregated at a country level, these data do not take into account high levels of spatial heterogeneity. To account for potential bias, future work could consider a sensitivity analysis on the level of individual cities or provinces where available. The delay between onset dates and reported dates for new case also potentially introduces bias and data spikes (as noted in “Methods” and “Results” sections). While we do apply any corrections for this reporting delay, techniques to resolve this problem could be considered from the studies of Influenza and Dengue epidemics [46]. Lastly, the Johns Hopkins data does is not curated to distinguish between cases acquired through local transmission as opposed to imported cases. Further modelling extensions that account for details captured in alternative data sources, such as the European Center for Disease Control and Prevention (ECDC) [3], should also be considered.

We focus on tracking model parameter estimates in time. Here, countries are included as outbreaks occur. Partitioning countries based on geographical location is another possibility for future exploration. Such separation may highlight overall differences in response patterns between different global regions. For example, how do the response trends differ as the pandemic moved from Asia to Europe and then to the Americas?

Our model, like any model, has fundamental assumptions that are necessarily introduced. We treat each country as a single well-mixed population. While our approach does include important features such as undocumented infections and a variable transmission rate, more advanced analysis could be performed by considering disease spread through a network [9, 47–49] of well-mixed populations, such as provinces, states or cities. This would assist is capturing social factors that could also influence COVID-19 transmission, such as spatial variation in population density and large population movement such as those that occur during times of festival [5, 40]. We also treat each country as a closed system, whereas realistic sources and sinks through inclusion of an international travel network could enable us to track the impact of decisions of one country on connected countries.

Other details of the model could also be extended. We treat active confirmed cases as non-infectious due to quarantine and isolation. In reality, active confirmed cases can still spread the virus to medical staff. We also apply the reasonable approximation that there is no re-detection or re-infection, however, new evidence is questioning the validity of this assumption [42, 45, 50]. It is also possible that seasonal effects related to climate could also cause transmission rates to change, although the evidence suggests that this is not a substantial effect in the pandemic stage [51]. A finer granularity of classes of susceptible individuals (i.e., at risk), incubation periods, severity of symptoms, and climate effects would also enable more detailed analysis for an individual country [17, 21]. However, our reduced set of classes and interactions represents a trade-off between realism and broad applicability to worldwide data. A Bayesian hierarchical modelling approach could also be applied to better capture heterogeneity across countries.

Future work should also consider methods that account for heterogeneity in terms of population distribution, age-structure, socio-economic factors and cultural aspects within individual countries. This is particularly true for continental scale countries with many diverse cultures like Brazil and India. While the net effects of many of these factors will be captured in the rate parameters α and α₀, bias in parameter inferences may occur without explicit consideration of non-uniform transmission patterns. As a result bias correction factors may need to be derived. Interactions with other epidemics may also introduce bias in countries that are experiencing other outbreaks beyond COVID-19, such as Brazil [52].

Our model framework is flexible through the inclusion of a response function (See “Methods” section) in the virus transmission mechanism and may be extended to other scenarios. In this manuscript, we have only considered the case of a response dependent on the number of active confirmed cases, leaving a sensitivity analysis for the more general from for the Supplementary Material. This response function could be further extended to include economical factors or be modified to be a function of state and time. This would enable a wide range of behaviours to be explored, such as specific timings of enforced NPIs, and subsequent lifting of restrictions when active case go below a threshold.

Conclusions

Our work confirms that a multi-pronged approach to combat COVID-19 is essential. Firstly, early introduction of testing and effective contact tracing protocols and quarantine effectively reduces the uncertainty in the unobserved infected population (i.e., low κ and high γ). Intervention strategies are also essential and are most effective when introduced early (high w_A). These results demonstrate the utility of modelling combined with high quality, immediately available data for providing insight into the early stages of the pandemic.

It is hoped that our work might be used to inform future responses to outbreaks of COVID-19 or other pandemics. The message is clear: to avoid multiple waves we must not be complacent in response to an outbreak as the earliest confirmed cases arise. We also highlight the importance of wider testing to effectively reduce uncertainty in predictions of case numbers, recoveries and deaths.