key: cord-0580334-eualiq2i authors: Yang, Xian; Wang, Shuo; Xing, Yuting; Li, Ling; Xu, Richard Yi Da; Friston, Karl J.; Guo, Yike title: Bayesian data assimilation for estimating epidemic evolution: a COVID-19 study date: 2020-12-22 journal: nan DOI: nan sha: 735521a54957cee8ac479b5b7aa364806cd6c37a doc_id: 580334 cord_uid: eualiq2i The evolution of epidemiological parameters, such as instantaneous reproduction number Rt, is important for understanding the transmission dynamics of infectious diseases. Current estimates of time-varying epidemiological parameters often face problems such as lagging observations, averaging inference, and improper quantification of uncertainties. To address these problems, we propose a Bayesian data assimilation framework for time-varying parameter estimation. Specifically, this framework is applied to Rt estimation, resulting in the state-of-the-art DARt system. With DARt, time misalignment caused by lagging observations is tackled by incorporating observation delays into the joint inference of infections and Rt; the drawback of averaging is overcome by instantaneously updating upon new observations and developing a model selection mechanism that captures abrupt changes; the uncertainty is quantified and reduced by employing Bayesian smoothing. We validate the performance of DARt and demonstrate its power in revealing the transmission dynamics of COVID-19. The proposed approach provides a promising solution for accurate and timely estimating transmission dynamics from reported data. retrospective analysis of real-world scenarios. Our system provides the insights of impacts of different intervention polices and highlights the effectiveness of undergoing mass vaccination. Epidemic modelling is important for understanding the transmission dynamics and responding to the emerging COVID-19 pandemic [1] - [8] . Since the pilot work by Kermack and McKendrick [9] , various epidemic models with different governing equations have been developed to describe the transmission dynamics of infectious diseases [10] . For common infection diseases such as influenza, the epidemiological parameters are related to the nature of the virus and treated as constants during the epidemic outbreak. These models are not applicable to the emerging COVID-19 pandemic where extensive government control measures have been implemented and revised. Due to the impacts of control measures, the epidemiological parameters (e.g., infection rates) linked to human behaviours could change substantially. In particular, the instantaneous reproduction number ! has drawn extensive attention, defined as the expected number of secondary cases occurring at time , divided by the number of infected cases scaled by their relative infectiousness [11] . Estimating such timevarying parameters from epidemiological observations (e.g., daily report of confirmed cases) is useful for nowcasting transmission [12] , for retrospectively assessing intervention impacts and for developing vaccine strategies [6] , [13] , [14] . All the applications depend on a reliable system estimating the time-varying parameters with accuracy and timeliness. Imprecise estimation or inappropriate interpretation could feed misinformation. Several systems [10] , [15] - [18] have been proposed to estimate the time-varying epidemiological parameters in practice; however, this remains a challenging task due to the following issues [19] : a. Lagging observations. Given a mathematic model of transmission dynamics, to infer the time-varying parameters, the number of infections should be the ideal observable data. However, the actual infection number is unknown and can only be inferred from other epidemiological observations (e.g., the daily confirmed cases). Such observations are lagging behind the infection events due to inevitable time delays between the infection of individual patients and the detection of the cases (e.g., days for symptom onset [20] ). Direct parameter estimation from lagging observations without adjusting for the time delay results in the temporal inaccuracy of estimates [19] , [21] . To address this problem, a two-step strategy, first estimating infections from epidemiological observations with a temporal transformation followed by parameter estimation, has been commonly used in practice [21] . The simple temporal shift of observations by the mean observation delay turns out not sufficient for the relatively long observation delay or the rapidly changing transmission dynamics, which are seen in the COVID-19 pandemic [21] . Backward convolution method (i.e., subtracting time delay, with a given distribution, from each observation time) leads to an over-smooth reconstruction of the infection number and bias for parameter estimation [10] . Deconvolution methods [22] through inversing the observation process are mathematically more accurate but sensitive to the optimisation procedures (e.g., stopping criterion) of the ill-posed inverse problem. In addition, the estimated result of infection number is often calculated as a point estimate so that the uncertainty from the observation process is neglected [23] . Taking an alternative method to the two-step strategy, we are investigating a new Bayesian approach that could jointly estimate both infection number and epidemiological parameters with uncertainty by explicitly parameterising the observation delay. b. Averaging inference. There are two general paradigms to deal with the challenge of estimating time-varying parameters: 1) reformulating the problem into an inference of static or quasi-static parameters, so that various methods for static parameter estimation can be used; 2) developing inference methods for explicit time-varying parameter estimation. For the first approach, the time-varying parameter is usually parameterised with several static parameters (e.g., the initial value and the exponential decay rate [17] ). The quasi-static method is to assume such a slow evolution of the parameter that could be treated as static within a short period. For example, Cori et al. [16] proposed a slidingwindow method 'EpiEstim' using a segment of observations for the averaging inference of ! , assuming ! remains the same within the sliding window. But this assumption does not apply to the rapidly changing transmission dynamics, and the window size affects the accuracy. Best practices of selecting the sliding window are still under investigation [21] . Instead of adopting a local sliding window, Flaxman et al. [13] defined several periods according to the dates of intervention measures, assuming a constant ! within each period. This approach requires additional information about the intervention timeline, which could be inaccurate, and does not capture the abrupt change of ! . In contrast to these windowbased methods, data assimilation [24] is a window-free alternative approach that has been less explored for parameter estimation in computational epidemiology. Applying sequential Bayesian inference [25] , [26] , data assimilation supports instantaneous updating of model states upon the arrival of new observations. The Bayesian model selection mechanism [27] can also be used for modelling the switching transmission dynamics under interventions, avoiding the drawback of averaging inference. Different from the common compartment models [28] - [30] used in concurrent data assimilation studies of COVID-19 modelling, we use the renewal process taking account for the changing infectiousness of the virus during the infection period. Moreover, we propose the Bayesian smoothing scheme that allows the correction of historical estimate based on subsequent observations. c. Quantification of uncertainty. The credibility of parameter estimation is of equal importance compared to the estimate itself, especially for policymaking. The uncertainty come from different sources, including the intrinsic uncertainties of epidemic modelling, data observation and inference processes. Firstly, the uncertainty of epidemiological models affects the final estimates. For example, ! estimation is found to be sensitive to the assumed distribution of generation time intervals [21] . Secondly, the uncertainty, resulting from systematic errors (e.g., weekend misreporting) and random errors (e.g., spike noise) in the observation processes should be properly quantified. During the COVID-19 pandemic, for example, we have seen different reporting standards and time delay across countries and regions, with different levels of uncertainty. Thirdly, the uncertainty could be enlarged or smoothed in the inference processes. For example, the use of a sliding window could smooth the parameter estimation but may simultaneously miscalculate the uncertainty, due to the overfitting within the sliding-window. To provide reliable credible intervals (CrI) of parameter estimates, the three aforementioned types of uncertainty should all be considered and reported as part of the final estimates. As a state-of-the-art package for ! estimation, EpiEstim (Version 2) [31] allows users to account for the uncertainty from epidemiological parameters by resampling over a range of plausible values. However, the uncertainty from imperfect observations and the side effects associated with the sliding window cannot be processed by this tool. Recently, 'EpiNow' [18] was proposed to integrate the uncertainty of observation process, but the inference is still based on the sliding window. In this work, we deal with model and data uncertainty in the data assimilation framework [24] with a Bayesian smoothing mechanism to enable both the latest and historical observations to be continuously integrated into inference flow, thereby alleviating spurious variability of estimations. In order to tackle these practical issues, we propose a comprehensive Bayesian data assimilation system, for estimating time-varying epidemiological parameters together with their uncertainty. In particular, we focus on the joint estimation of infection numbers and ! as a real-world application. Compared to the Bayesian approach for estimating the basic reproduction number " at the beginning stage of an epidemic break [32] , the sequential updating scheme is developed in our system. The evolution of the transmission dynamics is described by a hierarchical transition process, which is informed by newly data formulated with explicit observation delay. A model selection mechanism is built in the transition process to detect abrupt changes under interventions. We propose a Bayesian data assimilation approach, as illustrated in Figure 1 , to estimate the time-varying parameters based on epidemiological observations. This framework is applicable to various epidemic models when the governing equations and observation functions are available. Given an epidemic model (e.g., renewal process), we can construct a latent state ! at time which consists of the time-varying variables and parameters of the governing equations. The epidemiological observations #:% up to the latest observation time are made during the observation process of the latent state ! . The problem of estimating the timevarying parameters can be formulated as a Bayesian inference problem of ( ! | #:% ) for each time step . In contrast to inferring the 'pseudo' dynamics (i.e., reformulating into a static/quasi-static problem), our method directly estimates the 'real' dynamics by assimilating information from the observations for the epidemic model forecast. Please refer to the Method section for more detailed explanations. To apply the proposed Bayesian data assimilation approach in a real-world problem, we developed the 'DARt' (Data Assimilation for ! estimation) system for the ! estimation. The The performance of the DARt system is validated and compared to that of the state-of-the-art EpiEstim and EpiNow2 systems through simulations and real-world applications. The results confirm its power of estimation and adequacy for practical use. We have made the system available online for broad use in ! estimation for both research and policy assessment. § Validation through simulation Due to the lack of ground-truth ! in real-world epidemics, we conduct a set of simulation experiments by using synthetic data for validation. Figure 2 illustrates the design of simulation experiments where a synthesised ! , is adopted as the ground truth to validate its estimated 2 ! . We also estimated ! using the state-of-the-art ! estimation package EpiEstim [31] and EpiNow2 [18] to compare the effectiveness in overcoming the three aforementioned issues (i.e., lagging, averaging and uncertainty). Validation experiment of the DARt system on simulated data. First, the groundtruth ! sequence is synthesised using piecewise Gaussian random walk split by several abrupt change points. The sequence of incident infection ! is simulated based on a renewal process parameterised by the synthesised ! . The observation process includes applying a convolution kernel that represents the probabilistic observation delay to obtain the expected observation ̅ ! and adding Gaussian noise that represents the reporting error to obtain the noisy 'real' observation ! . The inputs (in grey) to the DARt system are the distributions of generation time, observation kernel and simulated noisy observation ! . The system outputs are the estimated 2 ! , estimated ̂! and change indicator 6 ! . These outputs are compared with the synthesised ! , ! and the time of abrupt changes. Also, the observation function is applied to the estimated ! to compute the estimated observation 7 ! with uncertainty, which is compared to the 'real' observation. In the simulation experiments, we compare the performance of DARt with that of two comparative methods: EpiEstim and EpiNow2. These two models are applied under their default settings with a 7-day smoothing window. EpiEstim takes observations ! as the infection without adjusting for the observation delay, which is a 'plug-and-play' use of EpiEstim as reported in [1] . As the current implementation of EpiNow2 1 only supports Gamma and Lognormal distribution as time delays, we set the generation time distribution to be a Gamma distribution with shape and scale equal to 4.44 and 1.89 (obtained by fitting the Weibull distribution reported by Ferretti et al. [3] using the Gamma distribution). With the simulated ! curve and the generation time distribution, we follow the renewal process to simulate the infected curve ! (initialized to be 1). Then, the observation curve of onset cases ̅ ! is generated using the incubation time distribution [3] (i.e., the lognormal distribution with log mean and standard deviation of 1.644 and 0.363 days respectively) as the observation time delay. Similar to the experiments in other related work [12] , [13] , all comparative models start estimation when the daily observation exceeds a threshold number, which is set to be 10 in our experiments. All simulation results are shown in Figure 3 and our observations are summarized as follows. • Correctness of estimation: Figure 3 (B) compares the synthesised " with the estimated " from DARt, EpiEstim and Epinow2. We can see that " from DARt matches the synthesised " • Correctness of estimation: Figure 3 (D) shows the simulated " , DARt estimated " , and EpiNow2 estimated " . We can find that the DARt estimated " with 95% CrI match well the simulated " . In contrast, the " curve from Epinow2 does not align with the simulated " very well. In particular, the peak value of " from EpiNow2 deviates greatly from the simulated value. • Accuracy in recovering observations : Figure 3 (E) compares the distributions of reconstructed " from DARt and that of EpiNow2. We can find that compared with " from EpiNow2, " from DARt with 95% CrI can generally match well with the simulated " . Epidemic dynamics in these four regions. The inference results for ! and reconstructed observations for these four countries are shown in Figure 4 . For Wuhan, the observation data are the number of onset cases compiled retrospectively from epidemic surveys, while for Hong Kong, UK and Sweden, the observation data are the number of reported confirmed cases. Notably, we use the onset-to-confirmed delay distribution from [33] together with the distribution of incubation time proposed in [3] to approximate the observation delay. As the ground-truth ! is not available, we validate the results by checking whether the estimated distributions of observation well cover the observation curve of ! . As shown in the top panel of each subplot in Figure 4 , the CrIs of estimated ! distributions (in blue) covered most parts of the original observations (in yellow), confirming the reliability of our ! estimation. Figure 4 (A) shows the results of testing using Wuhan's onset data [1] . We observe that there was a sharp decrease in Wuhan's ! after 21 st of Jan 2020, which is also illustrated by ! as the probability of abrupt changes peaked at this time (in green bars). A strict lockdown intervention has been enforced in Wuhan since 23 rd of January 2020. This sharp decrease in Wuhan's ! is likely to be the result of this intervention, indicating its impact. The small offset between the exact lockdown date and the time of sharp decrease might be due to noisy onset observations and approximated incubation time distribution. After lockdown, ! decreased smoothly, indicating that people's awareness of the disease and the precaution measures taken had made an impact. Since the beginning of February 2020, the value of ! remained below 1 for most of the time with the enforcement of quarantine policy and increases in hospital beds to accept all diagnosed patients. It is noted that the onset curve has a peak on 1 st February 2020, due to a major correction in the reporting standard. Neither ! nor ! curve from our model were severely suffered from this fluctuation, showing the robustness of our model thanks to the smoothing mechanism. The results from Wuhan suggest that our switching mechanism can overcome the issue of averaging and automatically detect sharp changes in epidemiological dynamics. The results from DARt are also compared with results from EpiEstim and EpiNow2. EpiEstim generates results with significant local variations and delays, while EpiNow2 can derive a smooth ! curve with no obvious delays. However, the immediate impact of lockdown cannot be well detected by EpiNow2. Figure 4 (B) shows the inferred results from Hong Kong that reported confirmed cases [34] during the most recent outbreak from November 2020 to March 2021. In Hong Kong, the number of infections remains low for the most time and the government has continuously imposed soft interventions. In the middle of November 2020, a newly imported case has triggered a new outbreak, resulting in a large ! value. However, the ! level returned to be around 1 very soon as the government has further tightened social distancing measures at that time. From late January 2021, Hong Kong started to implement mandatory lockdown in the restricted areas. Since then, the number of daily cases remains at low level. Compared with the results from EpiEstim and EpiNow2, the results from DARt have similar trend with the others. The delays in the results of EpiEstim are still significant. We also investigate the performance of DARt using different types of observations and present the results of ! estimation using onsets and confirmed cases in the supplementary document. Figure 4 (C) shows the inference results from the United Kingdom's reported confirmed cases [35] . It is noted that the United Kingdom was one of the first countries in the world to authorise the emergency use of COVID-19 vaccines. Since early December 2020, the United Kingdom rolled out its COVID-19 mass vaccination programme. By mid-February 2021, the United Kingdom had successfully hit its target of 15 million first-dose COVID-19 vaccinations, encompassing the top four priority groups for vaccination. As of 22nd April 2021, the UK had reached its target of 33 million (63% adults) first-dose COVID-19 vaccinations and 11 million (21% adults) second-dose. After 3 months since the mass vaccination programme started, the number of infection cases continuously followed a downward trend. However, since further easements of COVID-19 restrictions in mid-May 2021, ! gradually increased. During Euro 2020 football match (from June 11 to July 11), the ! value remained above 1. An immediate decrease in ! occurred when Euro 2020 was finished and since then ! remained around 1. The results from DARt, EpiEstim and EpiNow2 are generally consistent. However, DARt has accurately detected and responded to the impact of the completion of Euro 2020 in mid-July. In addition to studying the whole country, we applied DARt to typical cities in England to investigate the local epidemic dynamics as well (please refer to section 4 of the supplementary document). To summarise, DARt has been applied to four different regions for investigating the transmission dynamics of COVID-19 to demonstrate its real-world applicability and effectiveness. Consistent with the findings in the simulation study, DARt has shown its advantages in the following aspects: 1) Instantaneity --DARt adopts a window-free sequential Bayesian inference approach to detect and indicate abrupt epidemic changes; 2) Robustness --with Bayesian smoothing, the ! curve from DARt is stable at the presence of observation noise; 3) Temporal accuracy --DARt performs a joint estimation of ! and ! by explicitly encoding the lag into observation kernels. In this paper, we have proposed a Bayesian data assimilation scheme for estimating the timevarying epidemiological parameters based on observations. To study a real-world application scenario, we focus on estimating ! and provide a state-of-the-art ! estimation tool, DARt, supporting a wide range of observations. In the DARt system, epidemic states can therefore be updated using newly observed data, following a data assimilation process in the framework of sequential Bayesian updating. For the model inference, a particle filtering/smoothing method is used to approximate the ! distribution in both forward and backward directions of time, ensuring the ! at each time step assimilates information from all time points. By taking the Bayesian approach, we have emphasised the uncertainty in ! estimation by accommodating observation uncertainty in likelihood mapping and introduced Bayesian smoothing to incorporate sufficient information from observations. Our method provides a smooth ! curve together with its posterior distribution. We have demonstrated that inferred ! curves can describe different observations accurately. Our work is not only important in revealing the epidemical dynamics but also useful in assessing the impact of interventions. The sequential inference mechanism of ! estimation takes into account the accuracy of time alignment and provides an abrupt change indicator. Different from approaches of directly incorporating interventions as co-factors into epidemic model [13] , [37] , our method offers a promising method for intervention assessment. The third approximation is implicit in the use of a particle filter to approximate the posterior distributions over model state variables -including ! -with a limited number of samples (i.e., particles). Particle filtering makes no assumption about the form of posterior distributions. On the contrary, the variational equivalent of particle filter, namely variational filtering [38] provides an analytical approximation to the posterior probability and can be regarded as limiting solutions to an idealised particle filter, with an infinite number of particles [39] . Considering the importance of both the mean value of ! and its estimation uncertainty is important for advising governments on policymaking, an analytical approximation is desirable to help properly quantify uncertainty is desirable. Finally, change detection is approximated by the change indicator ! , which is included as part of the latent state and inferred during particle filtering. This work opens an avenue to explore variational Bayesian inference for switching state models [40] . Crucially, variational procedures enable us to assess model evidence (a.k.a. marginal likelihood) and hence allow automatic model selection. Examples of Variational Bayes and model comparison to optimise the parameters and structure of epidemic models can be found in previous studies [41] . These variational procedures can be effectively applied to change detection. In conclusion, our work provides a practical scheme for accurate and robust estimation of timevarying epidemiological parameters. It opens a new avenue to study epidemic dynamics within the Bayesian data assimilation framework. We provide an open-source ! estimation package as well as an associated Web service that may facilitate other people's research in computational epidemiology and the practical use for policy development and impact assessment. The proposed Bayesian data assimilation framework for estimating epidemiological parameters include three main components: 1) a state transition model -describing the evolution of the latent state; 2) an observation function -defining an observation process and describing the relationship between the latent state and observations; 3) a sequential Bayesian engine: estimating statistical reason time-varying model parameters with uncertainty by assimilating prior state information provided by the transition model and the newly available observation. In this section, we introduce a real-world application of the proposed data assimilation framework to estimate one of the key epidemiological parameters, ! . The modelling epidemic dynamics is characterised by the renewal process, which is the foundation of our state transition model. We then describe the observation function, linking a sequence of infection numbers with the observation data. Next, we present a detailed state transition model and propose the sequential Bayesian update module. Common ! estimation methods include compartment model-based methods (e.g., SIR and SEIR [42] ) and time-since-infection models based on renewal process [31] . Their relationships are discussed in supplementary document Section 1.1. Comparative studies have been conducted in [21] to show that EpiEstim, one of the renewal process-based methods, outperforms other methods in terms of accuracy and timeliness. Given the renewal process, the key transition equation derived from the process is: where ! is the number of incident infection cases on day , . is the time span of the set { -}, and individualis the probability that the secondary infection case occurs days after the primary infection, describing the distribution of generation time [10] . The profile ofis related to the biological characteristics of the virus and is generally assumed to be timeindependent during the epidemic. Considering the simplicity and superior performance of applying the renewal process to model epidemic dynamics, our work adopts Equation (1) as the basic transition function for joint estimation of ! and ! . In epidemiology, the daily infection number ! cannot be measured directly but is reflected in observations such as the case reports of onset, confirmed infections and deaths. There is an inevitable time delay between the real date of infection and the date reporting, due to the incubation time, report delay, etc. Taking account of this time delay, we model the observation process as a convolution function between kernel , and the infection number in / most recent days. where ! is the observation data, andis the probability that an individual infected is detected on day . / is the maximum dependency window. It is assumed that the past daily infections before this window do not affect the current observation ! . Since there is a delay between observation and infection, we suppose the most recent infection that can be observed by ! is at the time − , where is a constant determined by the distribution of observation delay. To accommodate various observation types (e.g., the number of daily reported cases, onsets, deaths and infected cases), DARt will choose the appropriate time delay distributions accordingly. For example, for the input of onsets, the infected-to-onsets time distribution is chosen to be the kernel in the observation function. For the input of daily reported cases, the infected-to-onset and the onset-to-report delays are used together as the kernel in the observation function. These delay distributions can be either directly obtained from literature or inferred from case reports that contain individual observation delays [3] . Detailed descriptions of the observation functions for different epidemic curves can be found in the supplementary document (Section 1.2). In Figure 5 is a vectorised form of infection numbers ! , * indicates the most recent infection that can be detected at time is from the time * due to observation delay, and 2 is the length of the vector ! * such that ! is only relevant to ! * and ! * )# only depends on ! * via the renewal process. § State transition model In our model, indirectly observable variables ! and ! are included in the latent state. The state transition function for ! is commonly assumed to follow a Gaussian random walk [18] or constant within a sliding window as implemented in EpiEstim. Such a simplification cannot capture an abrupt change in ! under stringent intervention measures. To capture such abrupt changes, we introduce an auxiliary binary latent variable ! to indicate the switching dynamics of the epidemiological parameters under interventions without assuming a pre-defined evolution pattern (e.g., constant or exponential decay). ! = 0 indicates a smooth evolution corresponding to minimal or consistent interventions; ! = 1 indicates an abrupt change of corresponding to new interventions or outbreak. The smooth evolution is modelled as a Gaussian random walk while the abrupt change is captured through resetting the parameter memory by assuming a uniform probability distribution for the next time step of estimation. Doing so provides an automatic way of framing a new epidemic period that was manually done in [13] . The transition of ! is modelled as a discrete Markovian process with fixed transition probabilities controlling the sensitivity of change detection: where ( !&# , 3 * ) is a Gaussian distribution with the mean value of !&# and variance of 3 * , describing the random walk with the randomness controlled by 3 . U[0, !&# + ∆] is a uniform distribution between 0 and !&# + ∆ allowing sharp decrease while limiting the amount of increase. This is because we assume that ! can have a significant decrease when intervention is introduced but it is unlikely to increase dramatically as the characteristics of disease would not change instantly. The transition of the change indicator ! , is modelled as a discrete Markovian process with fixed transition probabilities: where is a value close to and lower than 1. The above function means that the value of ! is independent of !&# , while the probability of Mode II (i.e., ! = 1) is quite small. This is because it is unlikely to have frequent abrupt changes in ! . For the incident infection ! , the state transition can be modelled based on Equation (1) Therefore, all the historical information needed to infer ! * is available from ! * &# , i.e., ! * only depends on ! * &# (i.e., being Markovian). The state transition process and observation process are illustrated in Figure 6 . The latent state in our model is then defined as ! =< ! * , ! * , ! * >, which contribute to ! at time . The state transition function of ! * is therefore Markovian: The most recent infection that can be observed by ! is at the time * = − where is a constant determined by the distribution of observation delay. Suppose 2 is the length of the vector ! * = [ ! * &% $ )# , ! * &% $ )* , … , ! * ] such that ! is only relevant to ! * and ! * only depends on ! * &# via the renewal process. Therefore, 2 ≥ max( . , / − + 1). The case that 2 = / − + 1 is depicted in this figure. § Forward Filtering We formulate the inference of the latent state ! =< ! * , ! * , ! * > with the observations ! as within a data assimilation framework. A sequential Bayesian filtering approach is adopted to infer the time-varying latent state, which updates the posterior estimation using the latest observations following the Bayes rule. This approach differs from the fixed prior in the Bayesian inference of static parameters. This filtering mechanism computes the posterior distribution of the latent state by assimilating the forecast from the forward transition model with the information from the new epidemiological observations. For the implementation of this Bayesian updating process, we adopt a particle filter method [26] to efficiently approximate the posterior distribution through Sequential Monto Carlo (SMC) sampling. This eschews any fixed-form assumptions for the posterior -of the sort used in variational filtering and dynamic causal modelling [38] . Let us denote the observation history between time 1 and as #:! = [ # , * , … , ! ]. Given that previous estimation ( !&# | #:!&# ) and new observation ! , we would like to update the estimation of ! , i.e., ( ! | #:! ) following the Bayes rule with the assumption that #:! is conditionally independent of #:!&# given ! : where ( ! | #:!&# ) is prior and ( ! | ! ) is the likelihood. The prior can be written in the marginalised format: where ! is assumed to be conditionally independent of #:!&# given !&# , and the transition ( ! | !&# ) is defined in Equation (6) based on the underlying renewal process. The likelihood ( ! | ! ) can be calculated assuming the observation uncertainty follows a Gaussian distribution: where is the observation function with a kernel chosen according to the types of observations By substituting Equation (8) into Equation (7), we obtain the iterative update of ( ! | #:! ) given the transition ( ! | !&# ) and likelihood ( ! | ! ): . Compared with other parameter estimation methods [13] , [16] , Bayesian data assimilation takes the advantage of additional information to smooth inference results with reduced uncertainty caused by incomplete observations. More specifically, the smoothing mechanism can be described as: given a sequence of observations #:% up to time and filtering results ( ! | #:! ), for all time < , the state estimates are smoothed as: Comparing with the sliding-window (i.e., averaging inference) approaches, our sequential Bayesian updating with backward smoothing mechanism features an instantaneous epidemiological parameter estimation and smoothing uncertainty through the utilisation of all available observations. More details can be found in the supplementary document. We obtained daily onset or confirmed cases of four different regions (Wuhan, Hong Kong, Sweden, UK) from publicly available resources [1] , [34] - [36] . For Wuhan, we adopted the daily number of onset patients from the retrospective study [1] ............................................................................................................ ............................................................................................................... 14 CITIES ........................................................................................................ 15 5 EXPERIMENTAL SETTING ..................................................................................................................... 17 6 SENSITIVITY ANALYSIS ..................................................................................................................... .......................................................................... 20 REFERENCES ................................................................................................................................................. 21 1 Mathematical Models Origin 0(" − ))/) (S10) Then we have the infectiousness profile ?(", )), representing the effectiveness rate at which an infectious individual with infection-age ) produces secondary cases at time ": ?(", )) = #(") ?()) C()) (S11) The corresponding instantaneous reproduction number ((") is derived from the integral of infectiousness profile ?(", )): ((") = < ?(", ))/) $ " (S12) It is the average number of people that someone infected at time " is expected to infect, if conditions remain unchanged (i.e. susceptible population, infectiousness rate, recovery rate). From the above derivation, we observe that the infectiousness profile ?(", )) and corresponding instantaneous reproduction number ( ! are composed of three factors: #("), ?()) and C()). #(") represents the depletion of susceptible individuals: the decline of #(") will reduce the susceptible population size leading to possible herd immunity. ?()) represents the infectiousness of individuals with infection-age ). This is related to biological (e.g. viral shedding) and behavioural (e.g. contact rates) factors. C()) represents the recovery rate of individuals with infection-age ): faster recovery will result in shorter infectiousness period and smaller reproduction number. All three factors, #("), C()) and ?()) can be altered by the implementation of control measures along with time. Decomposition of the Infectiousness Profile. ( ! is determined by the evolution of the infectiousness profile ?(", )) according to the Equation (S12). Further, the infectiousness profile ?(", )) can be rewritten as: ?(", )) = ((") M(", )) (S13) where M(", )) = ?(", ))/ ∫ ?(", ))/) is called the distribution of generation time, representing the probability distribution of infection events as a function of infection-age ). That is, the distribution of time interval between the primary infection and subsequent secondary infection. In principle, the distribution of generation time is time-varying due to the three factors in Equation (S11), which increases the complexity of parametric modelling. Most existing studies assume a time-invariant generation time distribution (i.e. M(", )) = M()) ) while the introduction of control measures results in the change of (("). Under this assumption, Equation (S10) can be rewritten as: This is the core formula for ( ! estimation from the infection data. That is, Or, we can have the corresponding discretised version: where P * is the time span of the set {M % }. This decomposition of infectiousness profile into ( ! and time-invariant generation time distribution M % is one of the fundamental formulae for ( ! estimation in the existing literature (e.g. the well-known package 'EpiEstim' [3] and this paper). Formulation of the Observation Function. The infection number 0 ! is the ideal data source for ( ! estimation according to Equation (S16). However, it is impossible to obtain the exact number of real-time infections through intensive screening. Instead, the infections are usually observed from the statistical reports of related events (i.e. epidemic curves) such as the daily report of confirmed cases, onset cases and death number. There is an inevitable delay between the occurrence of infecting events and the events being reported. That means these epidemic curves do not reflect the current incidence of infection 0 ! . We clarify this by formulating the observation function of the transmission dynamics. In the framework of data assimilation, 0 ! is the state variable of the dynamic epidemic system and its update is driven by the parameter ( ! as described by Equation (S16). The aggregated reports T ! (e.g., daily confirmed cases, deaths) are the observing results of the state variable through an observation function U: where U is the observation function and T ! is the observation result. Supplementary Figure 1) . The reports of onset cases are usually compiled retrospectively from epidemic surveys of confirmed cases, which can be represented as: where W % ., is the probability that the symptom onset occurs X days after the initial infection date for a reported case, and / , indicates the T ! , can only cover information of infections at least / , days before " . The value of / , is determined by W ., . The distribution W ., is determined by the biological factors of the virus and has been investigated in the previous reports [7] , which is considered time-independent. We use the symbol ⨂ to denote the convolution operation. The epidemic curve of daily confirmed cases T ! + is observed from where W % .+ is the probability that a confirmed case is reported X days after the initial infection date, and / + has a similar definition with / , . The distribution W .+ includes two parts: the time between infection to symptom onset W ., , and the time between symptom onset to reported confirmation W ,+ : The former part W ., is usually similar across regions while the latter time delay W ,+ varies a lot due to test policies and screening capabilities. C. Death Reports. The epidemic curve of death T ! is observed from where Yis the observed mortality rate of infected cases, W % .is the probability that a confirmed case is reported dead X days after the initial infection date, and / -has a similar definition with / , . Yand W .vary across different countries and periods due to capacities of treatment [5] . The time-varying renewal process can be formulated through the framework of state space hidden Markov models. The instantaneous reproduction number ( ! and daily incident infections 0 ! are the two latent variables of the state space models, whose dependence is described by Equation (S16). Consider two evolution modes of ( ! : emerging smooth changes when interventions are being steadily introduced/relaxed, and undergoing an abrupt change due to intensive interventions (e.g., lockdown). We introduce another latent variable Z ! to automate the switch between these two modes, which will be discussed in detail in the next section. The observations T ! are the observed results of 0 ! through Equation (S18), (S19) and (S21). We are interested in inferring the evolution of ( ! (along with 0 ! and Z ! ) upon the realtime update of observations T ! . As revealed in Equation (S14) and Equation (S17), the observations experience time delay with respect to the update of the latent state, due to the lagging and averaging effects of convolution in Equation (S18)-(S21). Thus, the changes of ( ! cannot be reflected in time, due to the incubation time and observation delay. In other words, accurate estimation of ( ! at time " relies on future observations, which imposes the challenges of timely estimation. Therefore, we focus on two inference aims: 1. Given the latest observation, how to give a near real-time estimate of ( ! and -equally if not more important -how to assess the uncertainty of the results? 2. Upon update of the real-time observations, how to modify estimations at all previous time steps and assess the uncertainties to make them more accurate taking into account the new information? These two aims correspond to the two fundamental problems in Bayesian updating, namely the filtering and smoothing problems to be discussed in the next section. [ ! =< ( ! * , ] ! * , Z ! * > is defined as the latent state observed by T ! at time ". Since there is a delay between observation and infection, we suppose the most recent infection that can be observed by T ! is at the time " * = " − /, where / is a constant determined by the distribution of observation delay. Suppose P 1 is the length of the vector such that T ! is only relevant to ] ! * via Equation (S18)-(S21) and 0 ! * only depends on ] ! * &) via the renewal process. We formulate the estimate of the latent state [ ! from the observed reports T ! as a data assimilation problem. A sequential Bayesian approach is adopted to infer timevarying latent state, which are composed of two phases: forward filtering and backward smoothing. The likelihood can be calculated assuming an observation is with Gaussian variance: where U is chosen accordingly to the types of reports and I + 3 is the variance of observation error that can be approximated empirically (detailed settings can be found in Supplementary Section 4). As the likelihood function has explicitly considered observation noise, the function is more robust to noise compared to Poisson likelihood (used in EpiEstim). The results of using the Poisson likelihood for the same synthesised data as depicted in Figure 3 For example, ( ! * may experience an abrupt decrease due to the lockdown policy on time " * . Under this circumstance, the epidemic history does not provide much information about the latest ( ! * , where we name it as Mode II corresponding to Z ! * = 1. Formally, the evolution of ( ! * is described by the switching dynamics conditioned on Z ! * : where e(( ! * &) , I 5 3 ) is a Gaussian distribution with the mean value of ( ! * &) and variance of I 5 3 , describing the random walk with the randomness controlled by I 5 . U[0, ( ! * &) + ∆] is a uniform distribution between 0 and ( ! * &) + ∆ allowing abrupt decrease while limiting the amount of increase. This is because we assume that ( ! * can undergo a big decrease when intervention is introduced but it is unlikely to have a dramatic increase in one day as the characteristics of disease would not change instantly. In our model, we assume a discrete Markovian chain process for Z ! * with the transition probabilities listed in Supplementary Table 1, meaning that the probability of having an abrupt change in ( ! * is low. This assumption is realistic as most of the time the ( ! * curve is undergoing smooth change. Supplementary Table 1 . Transition probabilities of Z ! * . Finally, we use the renewal process to provide transition of ] ! * 2) : , ] ! * &) Then by integrating out which provides the iterative calculation of G([ ! |T ):' ) from time P backwards to time ". Our model includes both 0 ! and ( ! into the latent state and jointly estimates them together. The Supplementary Figure 3 compares 0 ! from our method with 0 ! by using the back calculation function from EpiNow2 1 (for the synthesised data mentioned in the main manuscript). We can see that our method returns a 0 ! curve with CrI that is closer to the simulated infection curve. infection by DARt is drawn in black with 95% CrI. The ground-truth simulated infection is in red and the back calculated infection is in yellow. The integrals in the filtering problem (Equation (S24)) and the smoothing problem (Equation (S30)) are intractable, thus we introduce a Sequential Monte Carlo (SMC) method called 'particle filter' to infer the latent state [8] . In Monte Carlo method, the continuous distribution of a random variable z~{(F) is approximated by | independent samples with importance weights: where {(F) is an arbitrary probability distribution and | independent samples z ;~{ (F) are drawn from the distribution with the normalized importance weight ~;. u < ' (F) denotes the Dirac delta mass located at the i-th sample z ; . These discrete samples are also called 'particles' in particle method, whose locations and weights are used to approximate the intractable integral. If F is a time-dependent state variable, we can update the samples to approximate the distribution through the Sequential Importance Sampling (SIS) technique [9] . Therefore We applied DARt to monitor seven cities in England as shown in Supplementary Figure In our experiments, the generation time and observation delay distributions are adopted from the previous reports [6] , [10] . To truncate these distributions into a fixed length, we discard the time points with the kernel values smaller than 0.1 resulting in the length of ! ! as 7. The initial guess of " " at # = 0 is set to be uniformly distributed from 1 to 5. We set & = 0.1 for getting smooth " " change in Model I. In Model II, we set ∆= 0.5 for all regions. To implement the particle filter, the number of particles is set to 200 for approximating distributions. The variance of observation error & # $ is estimated empirically. We first calculate the 7-day moving average observations. By subtracting the moving average from the observation, we obtain a difference curve, approximating random observation fluctuations. The next step is to perform the 7-day moving average calculation again on the squared value of the difference curve, where the resulted curve is regarded as the observation error variance. Finally, we use a Gaussian distribution as the likelihood function (Equation (S25)), where its variance is approximated by the observation error variance curve. To investigate the robustness of different comparative methods, Figure 6 shows the " " Here, we would like to conduct experiments to investigate how different settings of the truncation threshold would influence the " " estimation results. The supplementary Figure 7 In this paper, we focus on dealing with uncertainties brought by noisy observations. However, we are also quite interested in how our model would behave when uncertainties in the distributions of generation time and observation delay exist. The simulation results in Figure 3 from the main manuscript have illustrated the performance of DARt. As indicated in the Supplementary Table 2 , the simulation results further quantitively confirms our findings by calculating the mean and standard deviation of estimation differences over time. We can easily observe that: 1) compared with EpiNow2 and EpiEstim, DARt achieves much smaller estimation errors; 2) compared with DARt without smoothing, the full implementation of DARt with smoothing can reduce estimation errors, showing that smoothing has greatly contributed to uncertainty reduction. The Supplementary Table 2 . Simulation results using data synthesized in the main manuscript. " " -mean/L " -mean and Rt-sd/L " -sd are the mean and standard deviation of the differences between synthesized " " /L " and estimated " " /L " . Association of Public Health Interventions With the Epidemiology of the COVID-19 Outbreak in Wuhan Substantial undocumented infection facilitates the rapid dissemination of novel coronavirus (SARS-CoV-2) Quantifying SARS-CoV-2 transmission suggests epidemic control with digital contact tracing Impact of non-pharmaceutical interventions (NPIs) to reduce COVID-19 mortality and healthcare demand', Imp. Coll. COVID-19 Response Team Impact assessment of non-pharmaceutical interventions against coronavirus disease 2019 and influenza in Hong Kong: an observational study Effect of non-pharmaceutical interventions to contain COVID-19 in China Modeling COVID-19 scenarios for the United States Modelling transmission and control of the COVID-19 pandemic in Australia A contribution to the mathematical theory of epidemics Estimating individual and household reproduction numbers in an emerging epidemic Practical considerations for measuring the effective reproductive number, Rt Nowcasting and forecasting the potential domestic and international spread of the 2019-nCoV outbreak originating in Wuhan, China: a modelling study Estimating the effects of non-pharmaceutical interventions on COVID-19 in Europe Model-informed COVID-19 vaccine prioritization strategies by age and serostatus Different epidemic curves for severe acute respiratory syndrome reveal similar impacts of control measures A new framework and software to estimate time-varying reproduction numbers during epidemics Real Time Bayesian Estimation of the Epidemic Potential of Emerging Infectious Diseases Estimating the time-varying reproduction number of SARS-CoV-2 using national and subnational case counts A guide to R-the pandemic's misunderstood metric Temporal dynamics in viral shedding and transmissibility of COVID-19 Practical considerations for measuring the effective reproductive number, Rt Reconstructing influenza incidence by deconvolution of daily mortality time series A Bayesian Updating Scheme for Pandemics: Estimating the Infection Dynamics of COVID-19 Data assimilation: methods, algorithms, and applications Bayesian filtering: From Kalman filters to particle filters, and beyond A tutorial on particle filtering and smoothing: Fifteen years later Bayesian model selection for complex dynamic systems Bayesian sequential data assimilation for COVID-19 forecasting Sequential Data Assimilation of the Stochastic SEIR Epidemic Model for Regional COVID-19 Dynamics Bayesian dynamical estimation of the parameters of an SE(A)IR COVID-19 spread model Improved inference of time-varying reproduction numbers during infectious disease outbreaks Bayesian estimation of the basic reproduction number in stochastic epidemic models First-wave COVID-19 transmissibility and severity in China outside Hubei after control measures, and second-wave scenario planning: a modelling impact assessment Latest local situation of COVID-19 Department of Health and Social Care, 'Coronavirus (COVID-19) in the UK: Cases Daily number of new reported cases of COVID-19 by country worldwide Estimating clinical severity of COVID-19 from the transmission dynamics in Wuhan Variational filtering Generalised filtering Variational learning for switching state-space models Testing and tracking in the UK: A dynamic causal modelling study The estimation of the effective reproductive number from disease outbreak data A contribution to the mathematical theory of epidemics A new framework and software to estimate time-varying reproduction numbers during epidemics Association of Public Health Interventions With the Epidemiology of the COVID-19 Outbreak in Wuhan Estimating the effects of non-pharmaceutical interventions on COVID-19 in Europe First-wave COVID-19 transmissibility and severity in China outside Hubei after control measures, and second-wave scenario planning: a modelling impact assessment Temporal dynamics in viral shedding and transmissibility of COVID-19 A tutorial on particle filtering and smoothing: Fifteen years later Sequential importance sampling for nonparametric Bayes models: The next generation Quantifying SARS-CoV-2 transmission suggests epidemic control with digital contact tracing Dr Simon Wang, at the Language Centre, HKBU, has helped improve the linguistic presentation of this manuscript.