key: cord-0628429-9ae26r2z
authors: Zhang, Yong; Yu, Xiangnan; Sun, HongGuang; Tick, Geoffrey R.; Wei, Wei; Jin, Bin
title: COVID-19 infection and recovery in various countries: Modeling the dynamics and evaluating the non-pharmaceutical mitigation scenarios
date: 2020-03-31
journal: nan
DOI: nan
sha: af9fa74295a96130db158a3e963623671a079286
doc_id: 628429
cord_uid: 9ae26r2z

The coronavirus disease 2019 (COVID-19) pandemic radically impacts our lives, while the transmission/infection and recovery dynamics of COVID-19 remain obscure. A time-dependent Susceptible, Exposed, Infectious, and Recovered (SEIR) model was proposed and applied to fit and then predict the time series of COVID-19 evolution observed in the last three months (till 3/22/2020) in various provinces and metropolises in China. The model results revealed the space dependent transmission/infection rate and the significant spatiotemporal variation in the recovery rate, likely due to the continuous improvement of screening techniques and public hospital systems, as well as full city lockdowns in China. The validated SEIR model was then applied to predict COVID-19 evolution in United States, Italy, Japan, and South Korea which have responded differently to monitoring and mitigating COVID-19 so far, although these predictions contain high uncertainty due to the intrinsic change of the maximum infected population and the infection/recovery rates within the different countries. In addition, a stochastic model based on the random walk particle tracking scheme, analogous to a mixing-limited bimolecular reaction model, was developed to evaluate non-pharmaceutical strategies to mitigate COVID-19 spread. Preliminary tests using the stochastic model showed that self-quarantine may not be as efficient as strict social distancing in slowing COVID-19 spread, if not all of the infected people can be promptly diagnosed and quarantined.

The novel coronavirus disease 2019 (COVID-19) outbreak, a respiratory illness that started (first detected) in late December 2019, is a pandemic infecting >336,000 people in more than 140 countries with the average fatality rate of 4.4% globally (data up to 3/22/2020) [1] . The COVID-19 pandemic is infiltrating almost every aspect of life, damaging global economy and altering both man-made and natural environments. Urgent actions have been taken but further effective and efficient strategies are promptly needed to confront this global challenge. To address this challenge and promptly guide the next efforts, it is critical to model the transmission, infection, and recovery dynamics of COVID-19 pandemic. Mathematical models are among the necessary tools to quantify the COVID-19 dynamics and are the primary objective motivating this study.

This study aims to model the COVID-19 evolution for representative countries with apparent coronavirus cases, including China, United States (U.S.), Italy, Japan, and South Korea.

COVID-19 spread in these countries under different starting (initiation and detection) times. For example, China had passed the peak of the coronavirus outbreak and hit a milestone with no new local infections on 3/19/2020 (79 days from the onset, 12/31/2019, in Wuhan, China), while the U.S. coronavirus cases soared past 10,000 on the same day. This study will apply the core characteristics of COVID-19 outbreak obtained in China to estimate the COVD-19 spread in the U.S., as well as other countries where the number of affected people has not yet reached its peak.

In addition, this new pandemic may last for a relatively longer time than expected [2] . No vaccine against SARS-CoV-2 (severe acute respiratory syndrome coronavirus 2) is currently available [3] . Indeed, a vaccine for prevention and infection control may not be ready before March 2021 for COVID-19, considering the minimum 4~5 weeks for trial and at least 1 year for safety evaluation and final deployment. Efficient strategies are therefore needed to mitigate the 4 COVID-19 outbreak. Possible non-pharmaceutical scenarios, such as isolation of cases and contact tracing, can be evaluated using mathematical models [4] , to identify the most efficient strategy going forward. This is another major motivation and the secondary task of this study.

To address the questions mentioned above, this study is organized as follows. Section 2 proposes an updated SEIR model for COVID-19, where "S", "E", "I", and "R" stand for Susceptible, Exposed, Infectious, and Recovered people, respectively [5] . This model is then applied to fit and predict the COVID-19 spread in various provinces and major cities in China, resulting in abundant datasets to derive the core characteristics of the COVID-19 dynamics of transmission/infection and recovery. Section 3 predicts the spread of COVID-19 in the other countries using the knowledge gained from China. Section 4 proposes a fully Lagrangian approach to model the spatiotemporal evolution of COVID-19, and then applies it to evaluate non-pharmaceutical scenarios to mitigate the virus spread. Section 5 reports the main conclusions.

Several mathematical models have been applied for epidemic analysis of COVID-19. The most widely used one, so far, is the well-known SEIR model. For example, Peng et al. [6] (using data up to 2/16/2020) proposed a generalized SEIR model to successfully estimate the key epidemic parameters of COVID-19 in China, and predicted the inflection point and ending time of confirmed COVID-19 cases. The SEIR model was also applied by Li et al. [7] (using the observed data up to 2/6/2020) to compare the effect of city lockdowns on the transmission dynamics in different cities in China. The SIR model with time-dependent transmission and recovering rates was used by Chen et al. [8] (using data up to 2/20/2020) to analyze and predict the number of confirmed cases of COVID-19 in China. The SIR model was extended by Wang et al. [9] (using data up to 2/12/2020) to incorporate various time-varying quarantine protocols for assessing 5 interventions on the COVID-19 epidemic in China. The SEIR model and its modifications were also successfully applied by others [10] [11] [12] [13] , mostly for assessing the early spreading of in China. Previous applications of the popular SEIR model, however, may contain high uncertainty since they had limited data access for only a short period of the COVID-19 outbreak.

As will be shown below, the COVID-19 dynamics have changed dramatically in the last three months, likely due to the improvement/adjustment of screening/testing techniques, public hospital system capabilities, and the government's control policies for contagious diseases. An updated, much better version of the SEIR model is therefore needed and can be reliably built now for China, since the detailed datasets are now available and the coronavirus outbreak (i.e., # of people infected) has passed its peak in China.

Other models have also been used for specific purposes related to COVID-19, including the global metapopulation disease transmission model to project the impact of travel limitations on the epidemic spread [14] , the transmission model for risk assessment [15] , and the synthetic contact matrix model for reproductive ratios for COVID-19 [16] . To the best of our knowledge, the classical stochastic epidemic models, such as the discrete or continuous time Markov chain model and the stochastic differential equation model, have not yet been applied for COVID-19 spread scenarios. A preliminary stochastic model for evaluating COVID-19 spread in a city will be developed and applied in Section 4.

This section focuses on the deterministic model. The classical SEIR model considers constant parameters, which may not well capture time-dependent dynamics. Hence, the SEIR model is updated in Section 2.1., and then tested in Section 2.2.

The classical SEIR model containing four populations (S, E, I, and R) takes the form [17] :

where is the stock of susceptible population; is the number of persons exposed to or in the latent period of the disease; is the stock of infected; is the stock of recovered population; r denotes the number of susceptible people whom the infected people contact daily; N is the sum of all the four groups of people ( ( ) + ( ) + ( ) + ( ) = =constant representing the constancy of population N); p is the constant rate of infection (i.e., representing the probability for the infected people to transform the susceptible people into infected ones); is the constant rate for the exposed person transformed into the infected one; and is the constant recovery rate (defining the speed for the infected person to be cured or expired).

To allow for possible time-sensitive rates for COVID-19 evolution, we revise model (1):

where represents the number of deaths (which is one component in I); is the rate for the healthy, susceptible person to be transferred to the infected one by the exposed people (note that COVID-19 patients in the incubation period might be contagious too); and is the number of 7 healthy susceptible people that are contacted by the exposed people daily. Now, the infection rates can change with time, and the infected persons are removed from the risk of infection also with a time-dependent rate of ( ). If the recovered individuals can return to the susceptible status due to for example loss of immunity, then the partial differential equation (PDE) (2d) for the time rate of change of needs one more (sink) term: / = ( ) − , where is the rate of the recovered individuals returning to the susceptible status.

We add the fractional-order PDE (2e) containing the death probability of (while the other

, which is the Caputo fractional derivative [18, 19] with order (0 < ≤ 1). When the order =1, model (2e) reduces to the classical integer-order PDE for the death evolution. The fractional PDE (2e) is used here for two reasons.

First, the evolution of deaths and cures may be characterized by a random process, since the exact time for the (recovered) person to be initially infected is unknown (i.e., the patient that died or was cured today may have been diagnosed yesterday or last week). Second, some patients may not be treated in time after being infected, making the death toll to evolve with a time memory. Therefore, we extend the classical mass-balance equation of death cases to the fractional PDE to characterize the random property and memory impact embedded in the temporal evolution of mortality. The fractional PDE (2e) and its classical version will be compared below using real data.

We apply the SEIR model (2) to fit the infection and recovery of coronavirus in China up to 3/7/2020, and then predict the future evolution ( Figure 1 ). [20] (with solutions shown by the dotted lines in Figure 1a ) accurately fits the observed data at the early stage, but then overpredicts the spread of COVID-19 observed after 2/12/2020. Dynamics of transmission for COVID-19, especially the recovery rate, therefore, changed in time in China, likely due to the time-dependent conditions in for example medical care mentioned above. A dynamic SEIR model, therefore, may be preferred for modeling COVID-19 spread in China.

In addition, compared to the fractional PDE (2e), the best-fit solution using the classical PDE / = for death evolution (see the black, dotted line in Figure 1a ) slightly overestimates the late-time growth of mortality. The actual death toll in Hubei province grew slower than that estimated by a constant rate model, indicating that the memory impact may affect the late-time dynamics of death that can be better captured by the fractional FDE (2e).

The best-fit solutions using model (2) fit the evolution of the infected and recovered populations well for the data recorded from Hubei province and three large cities closely related to Wuhan, China ( Figure 1 ). The model can also predict well the observed time series of COVID-19 spread from 3/8/2020 to 3/22/2020 for most places, except for Shanghai City (Figure 1d ). This is due to the oversea cases imported to Shanghai, whose number was increasing quickly after 3/3/2020, causing inconsistency of population and failure of the model. Shanghai Pudong

International Airport, one of the two airports located in Shanghai City, is the eighth-busiest airport in the world and the busiest international gateway of mainland China. If deleting the coronavirus 9 cases imported from overseas, model (2) can predict well the COVID-19 data in Shanghai ( Figure   2 ). Therefore, model (2) works well for various places in China, while external sources can easily break the internal evolution, especially the asymptotic status, of COVID-19 in China.

The resultant time-dependent recovery rate ( ) is depicted in Figure 3 , where the rate fitted by the latest observation data point in the fitting period (i.e., 3/7/2020) remains stable in the following prediction period. The best-fit recovery rate is the highest for Shanghai (except for the impulse of ( ) for Wenzhou discussed below), which is expected since Shanghai has the best public health system of all of these cities. Contrarily, Hubei shows the lowest recovery rate, likely due to its delayed response and the relatively limited public health capability at the beginning of the outbreak compared with Shanghai. February 2020, resulting in the sudden increase in the total recovery rate. In addition, a relatively large number of people working in Wuhan returned to Wenzhou in late January, and it appears that the improved efficient screening process successfully identified the number of infected cases. The new cases were ~29,000 from 1/24/2020 to 1/31/2020 in Wenzhou (with an average of 3,600 new patients per day), who were then immediately centralized for treatment. It appears that this fast response helped to alleviate the spread of coronavirus in Wenzhou. 10 

The best-fit parameters of model (2) are listed in Table 1 , and the initial values for each group of people are listed in Table 2 . We reveal three behaviors in model parameters. First, the best-fit "S"-shaped ( ) (Figure 3) can be described by the sigmoid function ( ) = /(1 + )

(where a, b, and c are factors), showing that the recovery rate increases exponentially before reaching stable. This increase is likely due to the healthcare facility and experience improved with time (in an accelerating rate) before reaching their asymptote or maximum capacity.

Second, the rates and probabilities (r, , p, and ) affecting the COVID-19

transmission/infection slightly change in space and remain stable for a given site ( Table 1 ). The small spatial fluctuation of these rates may be due to the similar strategies of local government for We also introduce an index C to quantify the infection severity of COVID-19 at different places:

where is the maximum number of cumulative infected people at the given site. A smaller C represents a greater infection severity of coronavirus. There is a power-law relationship between 11 the regional population N and the maximum cumulated number of infected people ( Figure   4 ). This empirical formula may be used to approximate the largest cumulative number of infections, which will be applied below in predicting the COVID-19 evolution outside of China where the coronavirus infection has not yet reached its peak number of cases.

Different countries are applying different modes for slow the COVID-19 spread. In the next sub-sections, we discuss several representative countries and then fit/predict the virus spread there.

To 

To decrease the acceleration of COVID-19 spread, Italy's mode is now similar to China:

lockdown of the full population. The predicted COVID-19 spread in Italy is plotted in Figure 6a .

Although Italy has followed China's mode of national isolation, the number of infected people increased rapidly from 3/14/2020 to 3/19/2020 (~495 new cases per day). To account for the delayed national quarantine compared with China, we decrease the C index (when also increasing the upper limit of the cumulative infection). The COVID-19 evolution prediction results show that there may be a turning point in the next two weeks when the current infected cases begin to decline. We also separate the death toll from the number of recovered cases.

South Korea's mode of combating the spread of COVID-19 is fast detection and tracking of the disease. South Korea is using efficient mobile diagnoses tests and accurate tracing of infected cases, to maintain a low death rate even with a large infected population. The mobile method can test 20,000 people per day (the maximum capability on 3/12/2020), and apps for cell phones and/or credit cards can accurately track the routes of infected people with the help of local government (without invasion of privacy), so that warnings can be immediately delivered to the general population to obviate the places with high risk. The current infected population may have passed 13 its peak number of cases around 3/20/2020, and the prediction shows that the COVID-19 outbreak may be well controlled in ~35 days from 3/22/2020 (Figure 6b ).

Japan's mode of combating the spread of COVID-19 is compatible with that of the U.S., in addition to other changes such as enhanced education/outreach and fast treatment to the infected cases. Specific policies include social distancing (which might be a key barrier to the spread of the novel coronavirus), personal hygiene, and quarantine of the infected cases. The current data and modeling results (Figure 6c) show that Japan has an efficient way so far to limit the maximum population infected and slow the spread of COVID-19, while this outbreak may last for a while.

The model-predicted COVID-19 spread in U.S., using the fitted infection and recovery rates from China ( Figure 5) , reveals the impact of one possible mitigation scenario for COVID-19: state coronavirus lockdowns, which have now been implemented by some states in the U.S. such as California and New York. Other non-pharmaceutical options can and should also be evaluated using mathematical models, considering the recent surge of infected cases in the U.S.

When the number of infected persons is initially small compared to that of susceptible people, the infected and susceptible people are not well mixed and hence the system is not homogeneous.

Under such conditions, a stochastic model is needed, as the deterministic, continuum models (such as the SEIR model) assume well-mixing of components for a homogeneous system [21] [22] .

Hence, this section develops and applies a stochastic model to evaluate the non-pharmaceutical scenarios for mitigating COVID-19 with a small number of initial infections.

The random walk based stochastic model for COVID-19 spread is analogous to a mixing-limited bimolecular reaction-based mechanism/condition [23] . When a reactant A particle 14 (representing a susceptible individual A) meets a reactant B particle (representing an infectious person B), a chemical reaction may occur if the collision energy is large enough to break the chemical bond (meaning that the susceptible person A may be infected if satisfying additional criteria such as A and B are close enough, and B touches his/her face after receiving coronavirus from A). Therefore, the condition of A being infected is not deterministic but rather a random, probability-controlled process. This probability is related to various factors, such as the duration that A and B are in contact, the infectivity rate, and the distance between the two people, which may be characterized parsimoniously by the interaction radius R that controls the number of reactant pairs (susceptible + infectious) in a potential reaction (infection) [23] . Hence, the core of the random walk stochastic model for COVID-19 spread is to define the interaction radius R.

The analogous development and similarities between bimolecular reactions and the SIR model can also be seen from their governing equations. The time-dependent SIR model takes the form [24] :

where ( ) and ( ) denote the transmission rate and recovery rate at time t, respectively. Therefore, following the argument in Zhang et al. [23] and Lu et al. [25] , we derive analytically the interaction radius R for the SIR model (4):

where V denotes the volume of the domain, ∆ is the time step in random walk particle tracking,

= [ ] / is the mass (or weight) carried by each A particle, [ ] is the initial concentration of A (which can be assumed to be the normalized value 1 here), and denotes the initial number of susceptible people. After defining the interaction radius R, the particle tracking scheme proposed by Zhang et al. [23] and Lu et al. [25] can be applied to model the transmission of coronavirus between the susceptible and infectious people.

In addition to pharmaceutical strategies including vaccine and therapeutic drug development, and herd immunity that may either take a while or have a high risk, non-pharmaceutical scenarios can be tested. Several particle-tracking based stochastic models were proposed recently [26] to evaluate non-pharmaceutical scenarios to mitigate coronavirus spreading in a city. Here we evaluate three related scenarios (described below) using the stochastic model proposed above.

Scenario 1: No special constraints. Assuming a city with 10,000 people and 4 initial coronavirus cases. This scenario does not put any constraints for any population. We assume that 16 the population distributes in the city randomly at the beginning, and then moves randomly. The spread of COVID-19 is then simulated using the stochastic model proposed in Section 4.1.

Scenario 2: Social distancing. This scenario assumes that all people in the city maintain a social distance, significantly decreasing the probability of infection.

Scenario 3: Attempted quarantine. This scenario creates a forced quarantine for infected individuals before the outbreak (which is expected to be the most efficient way of quarantine).

In the stochastic model, we assume that 10 days after being infected, the person will be removed because of being cured or expired (dead). This is because the median disease incubation period was estimated to be 5.1 days [27] . For simplicity purposes, the interaction radius R (6) remains constant, since the constant interaction radius was found to be able to efficiently capture the temporal variation of effective reaction rates in mixing controlled reactions [23, 25] . The initial number of A and particles is 10,000 and 4, respectively. The Lagrangian solutions of the COVID-19 outbreak for the three scenarios are depicted in Figure 7 .

Scenarios 1, 2, and 3 have a peak in the curve of newly infected people at time t=28, 65, and 32, respectively, showing that the virus spreads the fastest for the scenario without mitigation constraints (i.e., scenario 1, where the number of the total infected increases by one order of magnitude every 10 days in the rising limb), as expected. However, the value for this peak (scenario 1, =198 people) is lower than that for scenario 3 (=267 people), although the total number of the infected people for scenario 1 (9,336) is slightly larger than scenario 3 (9,322). This may be due to a greater separation of infection cases for the higher number of initial coronavirus carriers in scenario 1, which causes a lower and relatively flatter COVID-19 evolution peak than scenario 3. Scenario 2 has the lowest peak value (=121 people) and the most-delayed peak in the curve of new cases, and the total infection time is almost doubled compared to the other two 17 scenarios, indicating that people living with strict social distancing may also suffer from a much longer period of COVID-19 threat. It is also noteworthy that the overall trend of the solution of scenario 1 (initial surge without special constraints, Figure 7a ) is similar to that for Italy which had delayed response to the COVID-19 outbreak initially (Figure 6a) , and scenario 3 solution (a lower peak value and a longer duration due to social distancing, Figure 7b ) is similar to that for Japan which has been taken social distancing actions (Figure 6c) .

The simulated particle plumes plotted in Figure 8 reveal the subtle discrepancy between the three mitigation scenarios. Scenario 1 assumes that four initial cases were initially located on the right side of the city, while the whole population (10,000 susceptible persons) was distributed randomly in the 11 domain ( Figure 8a ). The trajectory of each person is assumed to follow (two-dimensional) Brownian motion with retention, to capture the random vector for each displacement and the random waiting time between two consecutive motions. The virus moved quickly from east to west (Figures 8b and 8c) , spreading over the whole city before all the infected people were cured or expired (dead) at time t=69 (Figures 8d) . A total of 664 susceptible people (6.6% of the total population) distributed randomly around the city were never infected.

Scenario 2 assumes that social distancing can reduce the infection probability, which can be characterized by a smaller reaction rate or a smaller interaction radius in our Lagrangian approach. (Figure 8i) can still cause the spread of coronavirus (Figures 8j~8l) . Self-quarantine, therefore, may not be as effective as maintained social distancing.

The COVID-19 pandemic radically impacts our lives, altering our daily patterns and interactions in an unprecedented way and rate. In this study, the transmission/infection and 

Johns Hopkins Coronavirus Resource Center

A mathematical model for the novel coronavirus epidemic in Wuhan

COVID-19 -New insights on a rapidly changing epidemic

Feasibility of controlling COVID-19 outbreaks by isolation of cases and contacts

A contribution to the mathematical theory of epidemics

Epidemic analysis of COVID-19 in China by dynamical modeling

The lockdown of Hubei province causing different transmission dynamics of the novel coronavirus (2019-ncov) in Wuhan and Beijing. medRxiv

A time-dependent SIR model for COVID-19

An epidemiological forecast model and software assessing interventions on COVID-19 epidemic in China

Effectiveness of control strategies for Coronavirus Disease 2019: a SEIR dynamic modeling study

Transmission dynamics model of 2019 nCoV coronavirus with considering the weak infectious ability and changes in latency duration

Nowcasting and forecasting the potential domestic and international spread of the 2019-nCoV outbreak originating in Wuhan, China: a modelling study

Preliminary prediction of the basic reproduction number of the Wuhan novel coronavirus 2019-nCoV

The effect of travel restrictions on the spread of the 2019 novel coronavirus (COVID-19) outbreak

Risk assessment of novel coronavirus COVID-19 outbreaks outside

Estimation of country-level basic reproductive ratios for novel Coronavirus (COVID-19) using synthetic contact matrices

Linear model of dissipation whose q is almost frequency independent-II

Impact of absorbing and reflective boundaries on fractional derivative modes: Quantification, evaluation and application

Estimation of the transmission risk of the 2019-nCoV and its implication for public health interventions

An introduction to stochastic processes with applications to biology

Discrete and continuous SIS epidemic models: A unifying approach

Evaluation and linking of effective parameters in particle-based models and continuum models for mixing-limited bimolecular reactions

Networks: An Introduction

Lagrangian simulation of multi-step and rate-limited chemical reactions in multi-dimensional porous media

Why outbreaks like coronavirus spread exponentially and how to "flatten the curve

From Publicly Reported Confirmed Cases: Estimation and Application

Report 2: Estimating the potential total number of novel Coronavirus cases in Wuhan City