key: cord-0963605-u6r8itdf authors: Dass, S. C.; Kwok, W. M.; Gibson, G. J.; Gill, B. S.; Sundram, B. M.; Singh, S. title: A Data Driven Change-point Epidemic Model for Assessing the Impact of Large Gathering and Subsequent Movement Control Order on COVID-19 Spread in Malaysia date: 2020-11-23 journal: nan DOI: 10.1101/2020.11.20.20233890 sha: 0a9a3d800b22c3547ea179f10b9c18d85ec21c94 doc_id: 963605 cord_uid: u6r8itdf The second wave of COVID-19 in Malaysia is largely attributed to a mass gathering held in Sri Petaling between February 27, 2020 and March 1, 2020, which contributed to an exponential rise of COVID-19 cases in the country. Starting March 18, 2020, the Malaysian government introduced four consecutive phases of a Movement Control Order (MCO) to stem the spread of COVID-19. The MCO was implemented through various non-pharmaceutical interventions (NPIs). The reported number of cases reached its peak by the first week of April and then started to reduce, hence proving the effectiveness of the MCO. To gain a quantitative understanding of the effect of MCO on the dynamics of COVID-19, this paper develops a class of mathematical models to capture the disease spread before and after MCO implementation in Malaysia. A heterogeneous variant of the Susceptible-Exposed-Infected-Recovered (SEIR) model is developed with additional compartments for asymptomatic transmission. Further, a change-point is incorporated to model the before and after disease dynamics, and is inferred based on data. Related statistical analyses for inference are developed in a Bayesian framework and are able to provide quantitative assessments of (1) the impact of the Sri Petaling gathering, and (2) the extent of decreasing transmission during the MCO period. The analysis here also quantitatively demonstrates how quickly transmission rates fall under effective NPI implemention within a short time period. where S(t), E(t), I(t) and R(t) represents, respectively, the susceptible, exposed, infected and recovered compartments representing the total number of individuals in each compartment at time t (here,ẋ(t) denotes the derivative of x(t) with respect to time t for x ∈ {S, E, I, R}), N is the total population size, and h(S, I) = β S(t) N I(t) is the rate of new infections (or, the number of new cases in the population). The parameters that govern the trajectory of the SEIR model are θ ≡ (β, δ, γ, i 0 , e 0 ) which are, respectively, the transmission rate (i.e., number of individuals in the population an infected person comes in contact with and successfully transmits the disease per unit time), the rate of incubation of the disease, the rate of infectiousness, the initial number of infectives and the initial number of exposed individuals. SinceṠ(t) +Ė(t) +İ(t) +Ṙ(t) = 0, it follows that S(t) + E(t) + I(t) + R(t) = N for all t. In (1)-(4), the incubation period, 1/δ, is inversely proportional to the incubation rate δ, 148 and similarly, the infectious period 1/γ is inversely proportional to the infectious rate γ. The correspondence between rates and exponential sojourn times is only approximately so 150 since it is not so straightforward to establish this correspondence with an individual-based 151 stochastic model given the non-linear nature of the process. Modifications to typical SEIR formulation (1) preprint (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. at least on the timescales over which the epidemic is observed; see, for example, [19, 20] 177 In order to provide a quantitative assessment of the impact of the MCO, we modify (1)- (4) preprint (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this this version posted November 23, 2020. ; https://doi.org/10.1101/2020.11.20.20233890 doi: medRxiv preprint observed infectiousness is given by the following system of ODEs: where the rate of new infections is now given by T * , is chosen so that observed daily cases fall either to the left or right of T * . We denote the observation window to the left of T * by W L which consists of dates from March 1st up to and including T * . The window to the right of T * is denoted by W R which consists of dates after T * up to and including April 28th, 2020. The change-point date T * is inferred from data; it is not taken as March 18th, 2020, the date of the start of MCO. Choosing T * in this data-driven way is justified as the impact of MCO on observed cases may not be realized immediately. The modified SEIR model with change-point T * is governed by the ODE system (6)-(12) for all time points t ∈ W L . For t ∈ W R , the same ODE system (6)-(12) is considered but now with a different set of parameter values in W R compared to W L , and with a new β o expressed as a function of t, The basic reproduction number, R 0 , is defined as the number of secondary infections caused by one primary individual during his/her infectious period. It is the most important quantitative indicator reported to assess whether the disease is in control or not. It is well-known that the threshold value of 1 for R 0 distinguishes between the situations 10 All rights reserved. No reuse allowed without permission. preprint (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this this version posted November 23, 2020. ; where a major epidemic occurs versus the disease dying out eventually in the population. In fact, several studies in the literature report a gradual decrease in R 0 after lockdown [22, 23, 24] . For the SEIR model in (1) -(4), R 0 is given by the well-known formula R 0 = β/γ. Time-varying measures of secondary infections per primary infected, R t , are also available in the literature. However, for the modified SEIR model presented in Section 2.3, R 0 and R t cannot be computed. Therefore, we resort to an alternative quantitative assessement of disease spread -the total number of infections (i.e., generational) caused by the introduction of one additional infectious individual into the infection process at time point t. This procedure is illustrated using the I o compartment. Based on (6)-(12), new values for the ODE system are calculated from time point t onwards with current state values serving as initializations of the ODE system for all except one compartment: For the I o -compartment, the current value I o (t) is replaced by I o (t) + 1/N as the initial value. The new rate of incidence is given by h(S * , I * o , I * c ) over the infectious period of the individual when the ODE system is propagated using (6) The increase in the rate of incidence by the introduction of this individual in the I o compartment is given by Similarly, the increase in the rate of incidence by the introduction of one infectious individual into the I c compartment at time point t can be calculated. This is denoted by ∆ t (I c ). The final increase in the incidence is the mixture where p is the proportion of exposed individuals who enter the I o compartment in (8). The likelihood relating components of the R o -compartment to the total number of daily cases, D t , t ∈ W L ∪ W R , is taken to be the negative binomial probability function, that is, 11 All rights reserved. No reuse allowed without permission. preprint (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this this version posted November 23, 2020. ; https://doi.org/10.1101/2020.11.20.20233890 doi: medRxiv preprint has mean ν ≡ q r/(1 − q) and variance qr/(1 − q) 2 ≡ ν + τ ν 2 ; in (17) instead. More discussion on this aspect is presented later. eters of the ODE system (6)-(12) for t ∈ W L and t ∈ W R , respectively. The collection 251 of all parameters is denoted by The prior on τ is taken to be U (a τ , b τ ) for the entire observation window W L ∪ W R . The prior on β o is elicited indirectly via a reparametrization: For W L , we take β o = 273 All rights reserved. No reuse allowed without permission. preprint (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this this version posted November 23, 2020. ; https://doi.org/10.1101/2020.11.20.20233890 doi: medRxiv preprint 13 All rights reserved. No reuse allowed without permission. preprint (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this this version posted November 23, 2020. ; https://doi.org/10.1101/2020.11.20.20233890 doi: medRxiv preprint preprint (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this this version posted November 23, 2020. ; https://doi.org/10.1101/2020.11.20.20233890 doi: medRxiv preprint respectively, the posterior of Θ can be derived using Bayes theorem as 15 All rights reserved. No reuse allowed without permission. preprint (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this this version posted November 23, 2020. ; https://doi.org/10.1101/2020.11.20.20233890 doi: medRxiv preprint where D = {D t : t ∈ W L ∪ W R } is the collection of observed daily cases from March 1st until April 28, 2020, L(Θ) is the entire likelihood for D and π(Θ) is the prior on Θ as described in Section 2.5.2. Bayesian inference is carried out using Monte Carlo importance sampling. A total of M samples, Θ i for i = 1, 2, · · · , M , are generated from the prior specification π(Θ). The likelihood, L(Θ i ), is computed for each Θ i and normalized to obtain weights for i = 1, 2, · · · , M . To compute L(Θ i ), one requires to numerically evaluate the solution of the ODE system (6)-(12). This is achieved using the deSolve package in R. The Bayesian computational algorithm described here is developed using R and the RStudio R user interface. Through this importance sampling step, an approximation to the Maximum-a-Posteriori (MAP) estimator of Θ is given bŷ with the approximation becoming more accurate as M → ∞. Resampling Θ i s with preprint (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. 17 All rights reserved. No reuse allowed without permission. preprint (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. preprint (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this this version posted November 23, 2020. ; https://doi.org/10.1101/2020.11.20.20233890 doi: medRxiv preprint that they first increase for time points t ≤ T * followed by a significant drop for t > T * . Hence, we conclude that an exponential rise in cases occurred right after the completion 22 All rights reserved. No reuse allowed without permission. preprint (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. preprint (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this this version posted November 23, 2020. ; https://doi.org/10.1101/2020.11.20.20233890 doi: medRxiv preprint authors would also like to thank the Director General of Health Malaysia for his permission 457 to publish this article. Appendix: Prior Elicitation on Remaining Parameters 459 We describe the prior elicitation for the initial number of infectious and exposed individ- preprint (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this this version posted November 23, 2020. preprint (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this this version posted November 23, 2020. ; https://doi.org/10.1101/2020.11.20.20233890 doi: medRxiv preprint -Detection 488 of A New Case Infected by The 2019 Novel Coronavirus (2019-nCoV) in Malaysia Latest Situation of the Coronavirus Disease 2019 (COVID-19) Infection in Malaysia Latest Situation of the Coronavirus Disease 2019 (COVID-19) Infection in Malaysia Sri Petaling Tabligh gathering remains Msia's largest Covid-19 cluster On the spread of epidemics in a closed heterogeneous population Epidemiological Models with Parametric Heterogeneity An introduction to compartmental modeling for the 583 budding infectious disease modeler A contribution to the mathematical theory of 587 epidemics Robust qualitative estimation of time-591 varying contact rates in uncertain epidemics Using Mathematical Models to Assess Responses to an Outbreak of an 595 Emerged Viral Respiratory Disease, Final Report to the Australian Government De-596 partment of Health and Ageing Content/mathematical-models Basic and effective reproduction 601 numbers of COVID-19 cases in South Korea excluding Sincheonji cases Report of the WHO-China Joint Mission on Coronavirus 606 Disease 2019 (COVID-19) Bayesian inference in an extended SEIR model with nonpara-610 metric disease transmission rate: an application to the ebola epidemic in sierra leone Bayesian non-parametric inference for stochastic 614 epidemic models using Gaussian Processes World Health Organization, Coronavirus disease 2019 (COVID-19) Situation Comparing SARS-CoV-2 with SARS-CoV 622 and influenza pandemics Malaysia to boost virus testing with S Korean test kits Incubation period of 2019 novel coron-630 avirus (2019-nCoV) infections among travellers from Wuhan, China The authors would like to acknowledge the Global Challenges Research Fund (GCRF)