A Bayesian Updating Scheme for Pandemics: Estimating the Infection Dynamics of COVID-19
Shuo Wang, Xian Yang, Ling Li, Philip Nadler, Rossella Arcucci, Yuan Huang, Zhongzhao Teng, Yike Guo
1 Abstract βEpidemic models play a key role in understanding and responding to the emerging COVID-19 pandemic. Widely used compartmental models are static and are of limited use to evaluate intervention strategies with the emerging pandemic. Applying the technology of data assimilation, we propose a Bayesian updating approach for estimating epidemiological parameters using observable information for the purpose of assessing the impacts of different intervention strategies. We adopt a concise renewal model and propose new parameters by disentangling the reduction of instantaneous reproduction number πΉ π into mitigation and suppression factors for quantifying intervention impacts at a finer granularity. Then we developed a data assimilation framework for estimating these parameters including constructing an observation function and developing a Bayesian updating scheme. A statistical analysis framework is then built to quantify the impact of intervention strategies by monitoring the evolution of these estimated parameters. By Investigating the impacts of intervention measures of European countries, the United States and Wuhan with the framework, we reveal the effects of interventions in these countries and the resurgence risk in the USA. Index Terms βCOVID-19, Data assimilation, Bayesian updating, Renewal process, Epidemiology, Non-pharmaceutical intervention. I. I NTRODUCTION n response to the COVID-19 pandemic, governments have taken non-pharmaceutical intervention measures. Common measures include travel restriction, school and non-essential business closure and social distancing, as well as early isolation of confirmed patients. Recently, as the first-wave epidemic peak has faded away in many countries, the accumulated observations of epidemic growth [1] and corresponding intervention policies [2] shed more insights on how the interventions worked. Meanwhile, many governments have switched into the phase to reopen economic and social activities, with attention on tamping down possible resurgences. However, the recent second-wave outbreak in some countries and regions (e.g. the United States, Hong Kong) alerts us to monitor the epidemic evolution carefully while intervention measures are being relaxed. * Corresponding author: Yike Guo (email: [email protected]).
Mathematical models play a key role in understanding and responding to the emerging COVID-19 pandemic [3]β[5]. Compartmental models (e.g. SIR, SIER) and time-since-infection models (i.e. renewal process-based models) are the two well-known approaches describing the underlying transmission dynamics [6], [7]. The compartmental models describe the transmission among sub-populations while the renewal process-based approach starts from the inter-individual transmission. Despite different nomenclatures and applications, each model contains parameters characterizing the epidemic dynamics. One of the most well-known parameters is the reproduction number π , which represents the average number of secondary cases that would be induced by an infected primary case [8]. This key parameter is related to the final epidemic size of infectious disease [9]. Intervention measures aim to maintain the reproduction number under one so that the epidemic can be contained along with time. Thus, estimation of time-varying π will reflect the impacts of intervention. The basic reproduction number π ! is the reproduction number at the beginning of the epidemic outbreak, when the susceptible population is approximately infinite and without intervention measures. When various intervention measures are being introduced, the instantaneous reproduction number π " (also called effective reproduction number) is of greater interest. To gain insights into epidemic evolution, most existing studies such as [3], [10] focus on estimating time-varying instantaneous reproduction number π " . π " is defined as the average number of secondary cases that would be generated by an infected primary case at a time π‘ when conditions remained the same thereafter [8], reflecting the real-time transmission dynamics. This could help governments to monitor the evolution of COVID-19 and update intervention policies accordingly [11]. However, the nowcasting of π " from reported data is not an easy task. Several approaches have been proposed to estimate π " with different advantages [12]β[14], but the timeliness and accuracy are still of concern. Nowcasting results are affected by different factors, such as assumptions of the epidemic models, statistical inference methods and uncertainty of data resources. Inappropriate interpretation or imprecise estimation of π " are Ling Li is with the School of Computing, University of Kent, Kent, UK. Yuan Huang and Zhongzhao Teng are with the Department of Pure Mathematics and Mathematical Statistics, University of Cambridge, Cambridge, UK. Yike Guo is with Hong Kong Baptist University, Hong Kong, China.
A Bayesian Updating Scheme for Pandemics:
Estimating the Infection Dynamics of COVID-19
Shuo Wang I π " , we introduce two complementary parameters, the mitigation factor ( π " ) captures the effect of shielding susceptible population (e.g. through social distancing), and the suppression factor ( π· " ) captures the effect of isolating the infected population (e.g. through quarantine) to stop virus transmission. We propose a novel method to estimate these parameters by taking the data assimilation approach of using Bayesian updating methods. We use daily reports of confirmed cases as the observation. A deconvolution method is used to build an observation function to estimate the infection cases by adjusting the incubation time and report delay. The evolution of the time-varying infectiousness profile (i.e. π " and π· " ) is estimated from the adjusted epidemic curve through a Bayesian approach of assimilation. Such a fine-grained infectiousness profile enables us to quantify the impacts of various intervention measures in a comprehensive way. The paper is structured as follows: We introduce the related work in Section II. In section III, we present the overview of a time-varying renewal process model where the two parameters π " and π· " are proposed. In section IV, we present in detail the Bayesian updating scheme for estimating the dynamic parameters. In section V, we develop a statistical analysis method of assessing the intervention impacts based on the estimated results and the report of intervention policies. In section VI, as applications of our approach, we investigate the impacts of intervention measures of European countries, the United States and Wuhan to illustrate the importance of this development. II. R ELATED W ORK
At the beginning of COVID-19 outbreak in Wuhan, China, compartmental models (e.g. SIR, SEIR model) have been used to investigate the epidemic dynamics [16]β[18], where the basic reproductive number was estimated from the models with static parameters. With the spread of COVID-19 worldwide, renewal process-based models (i.e. time-since-infection model) are also being widely used in the study of COVID-19. The R package βEpiEstimβ [12], [13] is the most widely used in estimating the time-varying π " with a sliding window. In [10], βEpiEstimβ was applied to infer π " via the discrete renewal process for policy impact assessment. Similar work has been done in [3] to infer π " using βEpiEstimβ from laboratory-confirmed cases in Wuhan and hence evaluated the impact of non-pharmaceutical public health interventions. The work in [11] has pointed out that the infection data is usually not available and death data was used as observation for π " updating. Instead of simply applying βEpiEstimβ, they estimated π " by employing the renewal equation as a latent process to model infections and connecting the infections to death data via a generative mechanism. However, the estimated π " is in a piecewise form and the number of changing points was assumed to be determined by the imposed interventions. [19] estimates π " from the death data as well while linking the disease transmissibility to mobility using the renewal equation. In general, [11] and [19] explicitly formulated the π " βs updating function by introducing external factors (e.g. interventions and mobility). Thus, the estimated π " curve is largely constrained by the factors that are considered in the model. Data Assimilation [20] lends itself naturally to this problem since it provides a framework to enable dynamically updating the model states and parameters when new observations become available while also taking into account model and observation uncertainty. Data assimilation technologies, such as Kalman filter and variational method [21], have been widely used in signal tracking, oceanology, environment monitoring and weather forecasting where physical models and observation data are assimilated to produce accurate prediction. Data assimilation for epidemiological modelling was first proposed in [22] where compartment models were used as the underlying model for assimilation. In [25] and [26], estimating time-varying parameters in the compartment models was further investigated. To the authorsβ best knowledge, our work is the first study of applying data assimilation to the renewal process-based model. III. E PIDEMIC MODELLING OF
COVID-19
TRANSMISSION
In this section, we propose a time-varying renewal process with two complementary parameters π " and π· " to model the evolving infectiousness profile. We adopted a time-varying renewal process for epidemic modeling. The renewal process [8] of infectious disease transmission is: πΌ(π‘) = + πΌ(π‘ β π) π½(π)ππ (1) where
πΌ(π‘) is the incident infection on time π‘ and π½(π) is the infectiousness profile. The infectiousness profile means a primary case who was infected π time ago (i.e. with the infection-age π ) can now generate new secondary cases at a rate of π½(π) , describing a homogenous mixing process. The infectiousness profile π½(π) is related to biological, behavioral and environmental factors. We can calculate the reproduction number π as the area under curve of π½(π) , which is the overall number of secondary cases infected by a primary case. Further, π½(π) can be rewritten as: π½(π) = π β π€(π) (2) where the unit-normalized transmission rate π€(π) is the probability density function of generation time, i.e. the interval between the primary infection and the secondary infection. In 3 the early stage without intervention, the infectiousness profile remains time-independent as the baseline π½ ! (π) which describes the transmission dynamics when the susceptible population is infinite. The corresponding π is the well-known basic reproduction number π ! . In reality, the infectiousness profile π½(π) will evolve with time π‘ , therefore we introduce π½ " (π) to address the change in its distribution caused by intervention measures. To quantify the impacts of intervention measures to the evolution of π " , we propose two factors: suppression and mitigation to disentangle the intervention effects. Here we use two complementary metrics π " and π· " modelling the suppression and mitigation factors respectively, as illustrated in Figure 1. The suppression effects mainly shorten the infectious period of the infected population, corresponding to the truncation of π½(π) along the horizontal axis. We use a time-varying parameter π· " to denote the effective infectious window induced by suppression. The mitigation effects attenuate the overall infectiousness by shielding the susceptible population, corresponding to the scaling on the vertical direction. We introduce another time-varying parameter π " to describe this attenuation effect induced by mitigation. Formally, we parameterize the evolution of the infectiousness profile as: π½ " (π) = 4π½ ! (π) β π " π < π· " " (3) Accordingly, the instantaneous reproductive number π " can be derived: π " = π " β + π½ ! (π)ππ $ ! ! (4) Therefore, the impact of intervention measures on π " reduction is disentangled: mitigation factor π " attenuates the overall infectiousness through shielding the susceptible population and suppression factor π· " shortens the infectious period through isolating the infected population. It is noted that the π " can be derived from π " and π· " which provide more mechanistic details about the evolution of the infectiousness profile. IV. A DAPTIVE P ARAMETER E STIMATION
We aim to develop a comprehensive framework to estimate parameters of renewal process models using Bayesian updating approach of data assimilation, especially the three key parameters: < π " , π " , π· " >. The estimation is essential for quantify the impacts of different interventions through monitoring the evolution of < π " , π " , π· " >. This framework contains building an observation function to map observations to model state, modelling and Bayesian updating as shown in Figure 2 and 3. By applying the observation function, we reconstruct the number of daily infections from reports of confirmed cases, taking into account the incubation time and report delay with a deconvolution algorithm. Then < π " , π " , π· " > is estimated through a Bayesian approach of data assimilation. A. Reconstruction of daily infection from reported cases
In data assimilation, model states and parameters can be updated using new observation data. It is important for parameter estimation that proper observation is chosen, and an observation function can be built which maps observations to a state variable (usually regarded as the output of the model). In this study, the observations we have chosen are from the reported number of confirmed cases. The model output is daily infection incidence through the renewal process. However, such observations experience an inevitable time delay between the actual infection time and the reporting date (Figure 2). This includes an incubation time (i.e. the period between infection and onset of symptoms) and confirmation period (i.e. the period between onset and officially reported after being tested). The confirmed cases reported on time π‘ were actually infected within a past period and the reported number is the convolution result of the historical daily infection. Here, we define an observation function to reconstruct the Fig. 2. Reconstruction of daily infection from the confirmed cases using deconvolution algorithms. The time delay between the infection and onset and report is demonstrated (top). The estimated distribution between infection and report is presented which is used for deconvolution (bottom). Fig. 1. Disentangling the reduction of reproduction number into mitigation and suppression factors. π (π) for π β {0, π} from infection to report (Figure 2). Denoting the epidemic curve of reported infection cases πΌ? %:β = {πΌ? % , πΌ? ( , β¦ , πΌ? " } and the epidemic curve of confirmed cases πΆ %:β = {πΆ % , πΆ ( , β¦ , πΆ " } , the reported infection with an observation process of past infections can be modelled as a Poisson process: πΆ " ~ ππππ π ππ(ππππ = J π (π)πΌ? ")**+" ) (5) Estimate the daily reported infection curve πΌ? %:β given the daily confirmed cases curve πΆ %:β and infection-to-confirmed time distribution π %:, is an ill-posed deconvolution problem and can be solved using Richardson-Lucy (RL) iteration method [25]. The initial guess πΌ? %:"! is the confirmed cases curve πΆ %:β shifted back by the mode of the infection-to-confirmed time distribution. Let πΆ? -. = β π (π)πΌ? ")*.*+" be the expected number of confirmed cases on day π of iteration π , and π " be the probability that a reported case resulting from infection on day π‘ will be observed as defined in [25]. Then the iteration of πΌ? " is computed by an expectation-maximization (EM) algorithm as: πΌ? "./% = πΌ? ". π " J π (π β π‘)πΆ " πΆ? ".-0" (6) A normalized π ( statistics is used as the stop criterion of the iteration: π ( = 1π J (πΆ? -. β πΆ " )πΆ? -." < 1 (7) where π is the total number of data points. It is of note that the reported number of confirmed cases constitute the lower bound of the real infection due to the lack of mass test and the existence of asymptomatic cases. However, as long as the detection rate remains consistent, the scaling of reconstructed data does not affect the following inference of transmission dynamics. B. Bayesian Updating for Parameter Estimation
Following the Bayesian updating approach of data assimilation, we propose an instantaneous estimation method. For the defined epidemiology renewal process, the daily incident infection πΌ " is the state variable and can be assimilated from the reconstructed infection data from observation. The evolution of the state πΌ " is governed by the renewal process with the time-varying infectiousness profile π½ " (π) , parameterized with π " and π· " . Here we present a Bayesian framework to monitor the evolution of π " and π· " using the daily reports of confirmed cases (Figure 3). Our updating scheme employs a two-level hierarchical model for the inference of time-varying parameters [26]. Let us denote the observed daily incidence of infection till time step π‘ as πΌ? %:β = {πΌ? % , πΌ? ( , β¦ , πΌ? " } . Suppose pTπ β)% |πΌ? %:β)% W is the estimated distribution of π = [π, π·] at time step π‘ β 1 . Under the assumption of consistent detection rates, the observed daily incidence πΌ? " also satisfies the renewal process. The low-level model predicts the observation (i.e. reconstructed daily infection) given a parameter set through the renewal process: pTπΌ? β Zπ β , πΌ? %:β)% W ~ ππππ π ππ(ππππ = J π½ " (π; π β )πΌ? ")*")%*2% ) (8) where a Poisson process of observing the infected cases is assumed. This describes the likelihood of observing the new incidence data given history observations and parameter value π β . The high-level model describes the evolution of the model parameters π " and π· " through transforming the joint distribution: pTπ β ZπΌ? %:β)%
W = T β pTπ β)%
ZπΌ? %:β)%
W (9) where
T(. ) is a transformation function defining the temporal variations of the π . The prior knowledge of parameter distribution is transferred to the next time step π‘ by the high-level model T . Under the scenario without interventions, the parameters π " and π· " fluctuate around the baseline values. Therefore, we can assume a random walk of π in the parameter space as the high-level model. The update of joint parameter distribution is by convoluting with a Gaussian kernel with variance π % . When the intervention is introduced on time π , the random walk of π is altered where the variance of the Gaussian kernel will become π ( . The transformation T(. ) is defined as:
Fig. 3. Illustration of the Bayesian updating framework for estimating suppression and mitigation factors. We employ a two-level hierarchical model: For each time step, the low-level model (i.e. renewal process) provides the likelihood of π ! , π· ! (green). The posterior (orange) is calculated through the element product of the likelihood and the prior (blue) from the previous time step. To generate the prior for next time step, we use the high-level model (i.e. the transformation T ) to induce the evolution of parameters. The high-level model is a piecewise gaussian random walk process where the fluctuations of π ! and π· ! differ before and after an intervention time. The instantaneous reproduction number π ! can be derived from the posterior distribution of π ! and π· ! . T β p(π) = ap(π) β K " (π) π‘ < πp(π) β K (π) π‘ β₯ π (10) where K " (π) and K (π) are the Gaussian kernels before and after the deployment of intervention at time π . This high-level model includes three hyperparameters: variances before and after intervention: π % and π ( , and the change-point time π . Let us denote the hyperparameters πΌ = [π % , π ( , π] . After seen the latest observation πΌ? β , the posterior estimation of π is update by the Bayes rule: pTπ β ZπΌ? %:β
W = T β pTπ β)%
ZπΌ? %:β)%
W β pTπΌ? β Zπ β , πΌ? %:β)% WpTπΌ? β ZπΌ? %:")%
W (11)
This step reflects the Bayesian principle in the key updating step in Kalman filtering [21]. Unlike the Kalman filtering method where uncertainty is explicitly modelled through a covariance matrix under the Gaussian assumption, we directly use posterior probability to capture the uncertainty of estimation. The posterior is usually intractable but can be approximated through grid-based methods. Given a set of hyperparameters πΌ - , the hybrid model evidence can be calculated as [26]: pTπΌ? %:β ZπΌ - W = + pTπΌ? %:β , π β ZπΌ - Wππ β (12) Finally, the posterior estimation pTπ β ZπΌ? %:β W can be averaged across the hyperparameter grids weighted by the hybrid model evidence. The posterior mean and confidence intervals of π " and π· " as well as the corresponding π " are obtained in a dynamic manner. The prior of π ! at the first timestep is set uninformative as a uniform distribution with the pre-set lower and upper limits (e.g., the upper limit for the European countries is set to 8 in the experiment). The shape of π½ ! (π) is adapted from the distribution of generation time interval π€(π) reported by Ferretti et al.[5] We applied the above framework to infer the epidemic evolution in 14 European countries, states in the US and Wuhan city, China in Section VI. The codes of the our framework is released as an open-source package (https://github.com/whfairy2007/COVID19_Bayesian). V. E VALUATION OF I NTERVENTION M EASURES
With the estimated results from the above Bayesian updating scheme, now we can perform statistical analysis between the evolution of the transmission dynamics and the implementation of intervention measures. The whole framework containing data reconstruction, dynamic modelling, Bayesian updating, statistical analysis is presented in Figure 4. In this section, we introduce the quantification of intervention measures and the statistical method. A. Data Source
For the observations, we use the aggregated data of publicly available daily confirmed cases of 14 Europe countries (Austria, Belgium, Denmark, France, Germany, Ireland, Italy, Netherlands, Norway, Portugal, Spain, Sweden, Switzerland and the United Kingdom) and 52 states of the United States from John Hopkins University database [1]. The data include the time series of confirmed cases from January 22nd to June 8 th th π " of intervention measures during the analysis period (accessed on June 9 th π " β€ π " β€ π " β€ π " β€
80% and Level 4: 80%< π " β€ π " . B. Calculation of intervention policy indices
We categorize the dates within our analysis period in European countries into five different response levels, based on the overall stringency index π " . To identify the representative measures of each response level, we calculate the quantification indices of the eight intervention measures. Descriptions of the eight intervention measures and the quantification methods are provided in [2]. For each intervention measure, the Oxford report provides an ordinal scale quantification π£ of the strength of j-th policy implementation and a binary flag π Fig. 4. Components of the quantification framework. The evolution of mitigation and suppression factors are estimated using the infection data reconstructed from the daily reported confirmed cases. Given the history of government responses, the impacts of intervention measures are quantified by correlating the inferred epidemic parameters to response levels. π‘ . Following similar practice use in the Oxford report, we normalize the implementation of each intervention measure as π = max (0, π£ + 0.5π β 0.5)π Γ 100% (13) where π is the maximum value of the indicator π . To assign a label of response level to each measure, we calculate the change of mean policy indices across different response levels. The response level with the largest increase is considered as the level that the measure belongs to (i.e. the measure is a representative measure of this response level). For example, the mean index of school closure showed the largest increase from Level 0 to Level 1, so we consider this is a representative measure of Level 1. The representative measures of each response level are listed in Table 1. C. Regression analysis of the intervention impacts
We performed a retrospective analysis of the time-varying transmission dynamics during different response levels in Europe countries. First, the evolution history of π " and the overall stringency index π " are obtained using the above framework. The stringency index π " is categorized into five response levels. We fit a log-linear mixed-effect model, where the logarithm of π " is the outcome variable and categorical stringency index is the predictor. The logarithm is used to obtain the intervention impacts on the relative change of π " [27]. We performed a partial-pool analysis by assuming the impacts of intervention measure (slopes) shared across all selected European countries while the basic reproduction number π ! (intercept) varies due to environmental and social factors. The regression formula is written as: ln π = π ! + J π * β π· + πΎ + π π = 1,2, β¦ ,14 (14) where π is the estimated reproduction number of j-th country, π ! is the fixed effect term of ln π ! and π * is the fixed effects of interventions in response level π . π· is the dummy variable that takes the value 1 if and only if the response status is at Level k. πΎ is the random effect term following zero-mean Gaussian which explains the difference of ln π ! across countries and π is the Gaussian error term. Equation 14 associates the relative changes in π to the fixed effects of response levels, and can be rewritten into its marginal form as: ln(1 + π β π ! π ! ) = J π * β π· *6*2% (15) Therefore, the relative change of π due to the intervention measures in k -th response level can be derived from π * (i.e. βπ /π ! = ππ₯π(π * ) β 1 ). Country-specific ln π ! can be estimated as π ! + πΎ at the Level 0. The statistical analysis is performed using the R package βlme4β. The fixed effect is considered significant with P value<0.05. The 95% confidence intervals (CI) are estimated using bootstrap method. The assumption of normality is checked by inspecting the quantile-quantile plot of the residuals. The same procedure is also applied to the analysis of π· " and π " to quantify the suppression and mitigation factors, respectively. The results are demonstrated in Table 1. VI. R ESULTS A. Validation on simulated data
We simulated an artificial epidemic outbreak with a time-varying infectiousness profile using renewal process. The generation time intervals were adapted from Ferretti et al.[5]. The simulation period includes 50 days and an intensive intervention measure is induced on day 35 altering the transmission dynamics. Before the intervention, the ground-truth π " followed Gaussian random walk with a mean of 2.5. After the intervention (50% π " reduction and 67% π· " reduction), the mean of π " was reduced to 0.5 (black line). We validate the effectiveness of our approach in capturing the sudden change of π " evolution induced by interventions, which is hard to be detected by traditional sliding window-based methods (Figure 5). We compared the results using our approach (red line with 95% confidence intervals) to the results computed by the R package βEpiEstim v2.2β [12] (blue) which is a sliding window-based method widely used for π " estimation. We observed that the ground-truth π " is well estimated within our confidence interval. In particular, the sharp change of π " caused by the intervention is captured immediately by our approach while there is a lag using the sliding window-based method. B. Evaluation of Intervention measures in Europe Countries
In this part, we applied the proposed framework to analyze the epidemic evolution in the 14 European Countries and also Wuhan. With the inferred < π " , π " , π· " >, we can then assess the impacts of intervention measures. Figure 6 demonstrates the reconstruction of daily infections in the UK from the reported confirmed cases. The infected-to-report delay between report and infected time is composed of the incubation period (a lognormal distribution with a mean of 5.5 days and a standard deviation of 2.1 days [5]) and the onset-to-report period (a gamma distribution with a mean of 4.9 days and a standard deviation of 3.3 days [10]). The blue bars in Figure 6 indicate the number of confirmed cases. After deconvolving the confirmed numbers using infected-to-report delay, we got the infected curve, which is colored in red in Fig. 5. Validation of the proposed Bayesian updating scheme. π " of the UK from the infected curve. The missing values in the infected curve are replaced by the average mean of the neighbouring numbers. green bar is the posterior mean of estimated π " . To quantitatively show the impacts of different strength levels of interventions, Table 1 summarizes the statistical analysis results of 14 European countries. It shows different reduction rates of < π " , π " , π· " > for different response levels. The relative reduction of < π " , π " , π· " > compared to the minimal response (Level 0 where π " is set to π ! ) was estimated for each response level. With soft response (Level 1), the corresponding intervention measures (e.g. school closure, quarantine of international arrivals from high-risk regions) are correlated with a relative reduction of π " by 35% showing both strong suppression effect ( π· " shortening 22%) and mitigation effect ( π " reduction 29%). With strong response (Level 2), the relative reduction of π " increases to 60% with a strong mitigation effect ( π " reduction 56%). But the suppression effect ( π· " shortening 26%) is similar to that of Level 1, indicating marginal incremental suppression effect. This observation shows a consistency with the aim of representative intervention measures on this level (e.g. cancelling public events, restrictions on gathering and internal movements) to reduce the contact rates among the population. The emergent response (Level 3) shows substantial relative reduction of reproductive number ( π " reduction 71%) with suppression ( π· " shortening 37%) and mitigation ( π " reduction 67%) effects, correlated to the intensive measures (e.g. workplace closure and stay-at-home requirements). A similar degree of reductions is found for Level 4 ( π " reduction 74%; π· " shortening 40%; π " reduction 70%) while the stringency of intervention measures is higher. We find that our estimated evolving patterns of π " and π· " correspond well to the serial strategies taken by some European countries, such as the βcontain-delay-lockdownβ route taken in the UK. Apart from the results of 14 European Countries, Figure 8 also shows the results of applying our method to the data from Wuhan, where the greens bars indicate the posterior mean of π " during the outbreak of COVID-19. We can see that at the early stage of the pandemic, the π " levels are above 1. After the lockdown intervention has taken effect, π " has experienced a sharp decrease from 23 rd Jan. When the centralized quarantine policy has been enforced from the beginning of February, the π " values then largely remain below zero (the spike around 14 th Feb is due to misreporting). Figure 9 compares the reductions in < π " , π " , π· " > for different response levels between European Countries and Wuhan. From the analysis of Wuhan data, the strong impact of lockdown is clearly demonstrated with the immediate relative reduction of π " by 58%. We also observed that the combination of lockdown, centralized quarantine and immediate admission of confirmed patients starting from Feb 2nd in Wuhan was associated with a more substantial relative reduction of π " with strong suppression and mitigation effects. Fig. 6. Reconstruction of daily infections from the report of confirmed cases in UK. The forward convolution on reconstructed data (black line) matches well with actual reported data (blue bars), validating the correctness of the deconvolution method. Fig. 7. Estimated evolution of transmission dynamics in UK. The black line represents the reconstructed daily infection number and the green bar is the posterior mean of estimated π ! . Fig. 8. Estimated evolution of transmission dynamics in Wuhan. The black line represents the reconstructed daily infection number and the green bar is the posterior mean of estimated π ! . Two major events (city lockdown measure from Jan 23 rd and centralized quarantine from Feb 2 nd ) are annotated with red arrows. C. Resurgence risks in United States
We also used the proposed framework to estimate the epidemic evolution in different states of the United States. We observed that, as of the week ending May 31 st , the averaged reproduction number π " in 30 states exceeds 1 (Figure 10). These could be related to the recent lift of government restrictions and alert us to take a close monitoring on the epidemic evolution. At the time of preparing this paper (June 18 th th June 2020 have experienced an increased number of daily confirmed cases compared to that of May 31 st , and 14 states have recorded all-time high after May 31 st . When we prepare the final version in early August, this alarming prediction of a second wave outbreak is unfortunately proven true for all the states listed. So far, the application of the framework to many countries and the retrospective impact analysis of intervention measures in European countries indicate the effectiveness of our approach in monitoring π " . This can be further validated by predicting the evolution of π " , π· " and π " and projected infections in future study. Our current study has several limitations. Firstly, the reporting protocols and standards of confirmed cases, as well as the detection rates, vary among countries. However, as long as the reporting bias is consistent over time, the inference results of π " , π· " and π " should not be affected. We also note that the implementation of multiple intervention measures within a short interval makes it challenging to quantify the impact of a single measure which needs further statistical analysis. VII. C ONCLUSIONS
In conclusion, we propose a comprehensive Bayesian updating approach to timely estimate parameters of COVID-19 epidemic models. The disease transmission dynamics is modelled by renewal equations with time-varying parameters. Instead of purely focusing on estimating instantaneous reproduction number π " , we introduce two complementary parameters, the mitigation factor ( π " ) and the suppression factor ( π· " ), to quantify intervention impacts at a finer granularity. A Bayesian updating scheme is adopted to dynamically infer model parameters. By monitoring and analyzing the evolution of the estimated parameters, impacts of intervention measures in different response levels can be quantitatively assessed. We have applied our method to European countries, the United States and Wuhan, and reveal the effects of interventions in these countries and the resurgence risk in the USA. Our work opens a promising venue to inform policy for better decision-making in response to a possible second-wave outbreak. A CKNOWLEDGMENT
We express our sincere thanks to all members of the joint analysis team between Imperial College London, University of Cambridge and University of Kent and Hong Kong Baptist University. We thank Yuting Xing for helping collect epidemic data in Wuhan and the United States. We thank Siyao Wang and Liqun Wu for their efforts on developing a digital tracing app for validation and visualization.
Fig. 9. The relative reduction of mitigation factor π ! and suppression factor π· ! under different response levels compared to minimal response level. TABLE I. THE RELATIVE REDUCTION OF MITIGATION FACTOR AND SUPPRESSION FACTOR
UNDER DIFFERENT RESPONSE LEVELS OF E UROPEAN COUNTRIES
Response Representative Measures Impact of measures π ! relative reduction Suppression effect π· ! relative reduction Mitigation effect π ! relative reduction Level 0 Minimal response
No mandatory restrictions 0 0 0
Level 1 Soft response
Closing schools, International travel controls. 35% CI: [25%, 45%] 22% CI: [17%, 27%] 29% CI: [18%, 38%]
Level 2 Strong response
Cancel public events, Restrictions on gathering, Restrictions on internal movement. 60% CI: [54%, 65%] 26% CI: [21%, 30%] 56% CI: [50%, 61%]
Level 3
Close workplace, Close public transport, Stay-at-home requirements. 71% CI: [68%, 74%] 37% CI: [35%, 40%] 67% CI: [64%, 70%]
Level 4 Emergent response
74% CI: [71%, 77%] 40% CI: [37%, 42%] 70% CI: [66%, 73%]
9 R
EFERENCES [1] E. Dong, H. Du, and L. Gardner, βAn interactive web-based dashboard to track COVID-19 in real time,β
Lancet Infect. Dis. , vol. 20, no. 5, pp. 533β534, May 2020. [2] T. Hale, A. Petherick, T. Phillips, and S. Webster, βVariation in government responses to COVID-19,β 2020. [3] A. Pan et al. , βAssociation of Public Health Interventions With the Epidemiology of the COVID-19 Outbreak in Wuhan, China,β
JAMA , vol. 323, no. 19, p. 1915, May 2020. [4] R. Li et al. , βSubstantial undocumented infection facilitates the rapid dissemination of novel coronavirus (SARS-CoV2),β
Science (80-. ). , vol. 3221, no. March, p. eabb3221, 2020. [5] L. Ferretti et al. , βQuantifying SARS-CoV-2 transmission suggests epidemic control with digital contact tracing.,β
Science , vol. 6936, no. March, pp. 1β13, 2020. [6] E. Vynnycky and R. White,
An introduction to infectious disease modelling . OUP oxford, 2010. [7] N. C. Grassly and C. Fraser, βMathematical models of infectious disease transmission,β
Nat. Rev. Microbiol. , vol. 6, no. 6, pp. 477β487, 2008. [8] C. Fraser, βEstimating individual and household reproduction numbers in an emerging epidemic,β
PLoS One , vol. 2, no. 8, 2007. [9] J. Ma and D. J. D. Earn, βGenerality of the final size formula for an epidemic of a newly invading infectious disease,β
Bull. Math. Biol. , vol. 68, no. 3, pp. 679β702, 2006. [10] K. Leung, J. T. Wu, D. Liu, and G. M. Leung, βFirst-wave COVID-19 transmissibility and severity in China outside Hubei after control measures, and second-wave scenario planning: a modelling impact assessment,β
Lancet , vol. 395, no. 10233, pp. 1382β1393, Apr. 2020. [11] S. Flaxman et al. , βEstimating the effects of non-pharmaceutical interventions on COVID-19 in Europe,β
Nature , pp. 1β5, 2020. [12] R. N. Thompson et al. , βImproved inference of time-varying reproduction numbers during infectious disease outbreaks,β
Epidemics , vol. 29, no. August, 2019. [13] A. Cori, N. M. Ferguson, C. Fraser, and S. Cauchemez, βA new framework and software to estimate time-varying reproduction numbers during epidemics,β
Am. J. Epidemiol. , vol. 178, no. 9, pp. 1505β1512, 2013. [14] J. Wallinga and P. Teunis, βDifferent epidemic curves for severe acute respiratory syndrome reveal similar impacts of control measures,β
Am. J. Epidemiol. , vol. 160, no. 6, pp. 509β516, 2004. [15] D. Adam, βA guide to R-the pandemicβs misunderstood metric.,β
Nature , vol. 583, no. 7816, pp. 346β348, 2020. [16] N. Imai, I. Dorigatti, A. Cori, C. Donnelly, S. Riley, and N. Ferguson, βReport 2: Estimating the potential total number of novel Coronavirus cases in Wuhan City, China,β 2020. [17] Q. Li et al. , βEarly transmission dynamics in Wuhan, China, of novel coronavirus-infected pneumonia,β
N. Engl. J. Med. , vol. 382, no. 13, pp. 1199β1207, 2020. [18] J. T. Wu, K. Leung, and G. M. Leung, βNowcasting and forecasting the potential domestic and international spread of the 2019-nCoV outbreak originating in Wuhan, China: a modelling study,β
Lancet , vol. 395, no. 10225, pp. 689β697, 2020. [19] P. Nouvellet et al. , βReport 26: Reduction in mobility and COVID-19 transmission.β [20] M. Asch, M. Bocquet, and M. Nodet,
Data assimilation: methods, algorithms, and applications . 2016. [21] Z. Chen, βBayesian filtering: From Kalman filters to particle filters, and beyond,β
Statistics (Ber). , vol. 182, no. 1, pp. 1β69, 2003. [22] C. J. Rhodes and T. D. Hollingsworth, βVariational data assimilation with epidemic models,β
J. Theor. Biol. , vol. 258, no. 4, pp. 591β602, 2009. [23] L. M. A. Bettencourt and R. M. Ribeiro, βReal time bayesian estimation of the epidemic potential of emerging infectious diseases,β
PLoS One , vol. 3, no. 5, p. e2185, 2008. [24] L. Cobb, A. Krishnamurthy, J. Mandel, and J. D. Beezley, βBayesian tracking of emerging epidemics using ensemble optimal statistical interpolation,β
Spat. Spatiotemporal. Epidemiol. , vol. 10, pp. 39β48, 2014. [25] E. Goldstein, J. Dushoff, M. Junling, J. B. Plotkin, D. J. D. Earn, and M. Lipsitch, βReconstructing influenza incidence by deconvolution of daily mortality time series,β
Proc. Natl. Acad. Sci. U. S. A. , vol. 106, no. 51, pp. 21825β21829, 2009. [26] C. Mark, C. Metzner, L. Lautscham, P. L. Strissel, R. Strick, and B. Fabry, βBayesian model selection for complex dynamic systems,β
Nat. Commun. , vol. 9, no. 1, 2018. [27] A. Agresti,
An introduction to categorical data analysis . John Wiley & Sons, 2018. Fig. 10. The averaged π ! values in different states of the United States. We report the result of averaged π ! in the US during the week ending May 31st 2020, which is ranked by the averaged π !!