1.  Introduction {#s0005}
===============

After diagnosed in Wuhan of Hubei province, COVID-19 spread quickly from Wuhan to surrounding cities, regions outside Hubei in mainland China.  By January 22nd, 2020, 571 cases were reported by 25 provinces in mainland China ([@b0055]).  From January 23rd, 2020 Chinese government triggered a series of prevention and control strategies to prevent disease from spreading, including the lockdown of Wuhan, travel restrictions and calling for staying at home during Spring Festival holidays.  However, lots of infections or latent individuals have already been imported into regions outside Hubei because of the large volume of travels during Spring festival migrations, resulting in spatially spreading of COVID-19 epidemic in mainland China. 

Recently the COVID-19 epidemic in mainland China is weakening with enhanced prevention and control strategies.  When to return to work/study and how to manage the orderly work resumption become the challenging problems.  Some local governments have already divided the provinces into different areas according to the epidemic situation of the region itself, and the areas with low epidemic level are allowed to restore production or normal social order.  However, determining the date of work resumption is not only dependent on the epidemic situation of the region itself but also the epidemic of other regions as well as migration rate.  Then it is essential to accurately estimate the epidemic situations of a region and other regions and investigate the impacts of various dates of work resumption with weakened control measures on the potential outbreak. 

Predicting the epidemic of Hubei, estimating the number of cases in regions outside Hubei imported from Hubei province, and investigating the time of work resumption and lifting the lockdown in Hubei by linking reported data of Hubei and regions outside Hubei as well as the mobility of population between Hubei and regions outside Hubei are very important and fall within the scope of this study.  There are some early studies which have estimated the epidemic in Hubei or Wuhan by using the reported data and population movement data ([@b0105], [@b0010], [@b0030], [@b0070]).  Some mathematical modelling studies have estimated the epidemic in Hubei or mainland China without considering the mobility of population between Hubei and regions outside Hubei or without considering the time varying interventions ([@b0080], [@b0085], [@b0090], [@b0075], [@b0100]).  Note that the Wuhan government began to screen all the suspected cases on February 8 and increased the diagnosis rate from suspected cases greatly ([@b0025]).  How to quantify the changing or time-dependent detection rate or contact rate and evaluate the epidemic more precisely is challenging. 

In this study, we shall formulate a two-patchy mathematical model with one patch denoting Hubei province and the other representing regions outside Hubei in mainland China.  The two patches connected to each other because of the travelling of population between Hubei and regions outside Hubei in mainland China.  Based on the daily reported data for the epidemic in Hubei and regions outside Hubei, as well as the mobility data between Hubei and regions outside Hubei we shall focus on the estimation of the epidemic in Hubei and the effect of work resumption and lifting the lockdown at Hubei province on the epidemic in Hubei and regions outside Hubei in mainland China. 

2.  Methods {#s0010}
==========

2.1.  The model {#s0015}
--------------

A two-patchy model is formulated with one patch denoting Hubei province and the other being the regions outside Hubei in mainland China.  Before Jan 23th, 2020, these two patches are connected with each other due to the mobility of population between Hubei and regions outside Hubei.  Then, the two patches become isolated as the lockdown measure is implemented.  After lifting the lockdown in Hubei province, the two patches will be connected again.  For convenience, the subscript *h* and *c* represent the corresponding quantity in Hubei and regions outside Hubei in mainland China.  Susceptibles in Hubei $S_{h}$ (or regions outside Hubei $S_{c}$), may be infected when contacting with the infectious population in Hubei $I_{h}$ (or regions outside Hubei $I_{c}$).  After infected, the individual will firstly stay in the latent class, $E_{h}$ and $E_{c}$, and then become infectious $I_{h}$ and $I_{c}$. The infectious population may be diagnosed and confirmed with COVID-19 ($Q_{h}$ and $Q_{c}$) directly, or move to the suspected class $Y_{h}^{I}$ and $Y_{c}^{I}$, and then be confirmed with COVID-19.  Such kind of suspected individuals are in fact already infected, thus this class is also called as suspected and infected class.  The susceptibles may move to suspected but uninfected class $Y_{h}^{S}$ and $Y_{c}^{S}$ because of common fever or pneumonia, and the suspected but uninfected population may also move back to the susceptibles class after screening.  The confirmed population will either recover ($R_{h}$ and $R_{c}$) or dead.  The susceptibles, latent population and infectious population may travel between Hubei and regions outside Hubei at a certain probability.  Since the confirmed and suspected population will be quarantined, we do not consider the mobility rate for people in these classes. [ Fig.  1](#f0005){ref-type="fig"} shows the flow chart of the model, while Eqs. [( 1)](#e0005){ref-type="disp-formula"}, [(2)](#e0010){ref-type="disp-formula"} describe the model equations for Hubei and regions outside Hubei, respectively.$$\begin{array}{ll}
\frac{dS_{h}}{dt} & {= - \frac{\beta C_{h}(t)S_{h}I_{h}}{N_{h}} + \delta_{h}^{S}(t)Y_{h}^{S} - \mathit{mS}_{h} + r_{i}(t)S_{c} - r_{e}(t)S_{h},} \\
\frac{dE_{h}}{dt} & {= \frac{\beta C_{h}(t)S_{h}I_{h}}{N_{h}} - \sigma E_{h} + r_{i}(t)E_{c} - r_{e}(t)E_{h},} \\
\frac{dI_{h}}{dt} & {= \sigma E_{h} - \delta_{h}^{I}I_{h} + r_{i}(t)I_{c} - r_{e}(t)I_{h},} \\
\frac{dY_{h}^{s}}{dt} & {= \mathit{mS}_{h} - \delta_{h}^{S}(t)Y_{h}^{S},} \\
\frac{dY_{h}^{I}}{dt} & {= \delta_{h}^{I}(1 - q)I_{h} - \delta_{h}^{Y}(t)Y_{h}^{I},} \\
\frac{dQ_{h}}{dt} & {= \delta_{h}^{I}\mathit{qI}_{h} + \delta_{h}^{Y}(t)Y_{h}^{I} - \gamma_{h}Q_{h} - \alpha_{h}Q_{h},} \\
\frac{dR_{h}}{dt} & {= \gamma_{h}Q_{h}\text{.}} \\ 
\end{array}$$ $$\begin{array}{ll}
\frac{dS_{c}}{dt} & {= - \frac{\beta C_{c}(t)S_{c}I_{c}}{N_{c}} + \delta_{c}^{S}Y_{c}^{S} - \mathit{mS}_{c} + r_{e}(t)S_{h} - r_{i}(t)S_{c},} \\
\frac{dE_{c}}{dt} & {= \frac{\beta C_{c}(t)S_{c}I_{c}}{N_{c}} - \sigma E_{c} + r_{e}(t)E_{h} - r_{i}(t)E_{c},} \\
\frac{dI_{c}}{dt} & {= \sigma E_{c} - \delta_{c}^{I}I_{c} + r_{e}(t)I_{h} - r_{i}(t)I_{c},} \\
\frac{dY_{c}^{s}}{dt} & {= \mathit{mS}_{c} - \delta_{c}^{S}Y_{c}^{S},} \\
\frac{dY_{c}^{I}}{dt} & {= \delta_{c}^{I}(1 - q)I_{c} - \delta_{c}^{Y}Y_{c}^{I},} \\
\frac{dQ_{c}}{dt} & {= \delta_{c}^{I}\mathit{qI}_{c} + \delta_{c}^{Y}Y_{c}^{I} - \gamma_{c}Q_{c} - \alpha_{c}Q_{c},} \\
\frac{dR_{c}}{dt} & {= \gamma_{c}Q_{c}\text{.}} \\ 
\end{array}$$ Fig.  1The flow chart of the model.  The red dotted line from suspected and infected class to confirmed class in Hubei patch denotes the extra confirmation rate after Feb 8th, 2020 in Hubei province. 

The outbound travel rate from Hubei province to regions outside Hubei is $r_{e}(t)$, and the inbound travel rate from regions outside Hubei to Hubei province is $r_{i}(t)$ at time t. $\beta$ is the transmission probability per contact.  The contact rate in Hubei and regions outside Hubei are both piecewise functions which are proposed by Tang et al. ([ @b0085], [@b0090]), and denoted by $c_{h}(t)$ and $c_{c}(t)$, as shown in Eqs. [( 3)](#e0015){ref-type="disp-formula"}, [(4)](#e0020){ref-type="disp-formula"}, respectively.$$c_{h}(t) = \left\{ \begin{array}{ll}
{c_{h}^{0},} & {{before}\mspace{6mu}{Jan}\mspace{6mu} 23{rd},} \\
{(c_{h}^{0} - c_{h}^{b})e^{- r_{h}(t - t_{23})} + c_{h}^{b}.} & {{ after}\mspace{6mu}{Jan}\mspace{6mu} 23{rd}\text{.}} \\ 
\end{array} \right)$$ $$c_{c}(t) = \left\{ \begin{array}{ll}
{c_{c}^{0},} & {{before}\mspace{6mu}{Jan}\mspace{6mu} 23{rd},} \\
{(c_{c}^{0} - c_{c}^{b})e^{- r_{c}(t - t_{23})} + c_{c}^{b}.} & {{ after}\mspace{6mu}{Jan}\mspace{6mu} 23{rd}\text{.}} \\ 
\end{array} \right)$$

It is assumed that before Jan 23rd, the contact rate for Hubei and regions outside Hubei are both constants, denoted by $c_{h}^{0}$ and $c_{c}^{0}$, respectively.  After Jan 23rd, both contact rates decrease exponentially with the exponential decreasing rates being $r_{h}$ and $r_{c}$. And the smallest contact rates are $c_{h}^{b}$ and $c_{c}^{b}$ for Hubei and regions outside Hubei.  For convenience, initial and smallest transmission rates are further defined as $\beta_{h}^{0}\ \triangleq\ \beta c_{h}^{0},\beta_{h}^{b}\ \triangleq\ \beta c_{h}^{b}$ and $\beta_{c}^{0}\ \triangleq\ \beta c_{c}^{0},\beta_{c}^{b}\ \triangleq\ \beta c_{c}^{b}$ for Hubei and regions outside Hubei.  The latent period is assumed to be $1/\sigma$ after infected.  The diagnosis rates for infectious population are $\delta_{h}^{I}$ and $\delta_{c}^{I}$, among those diagnosed population a proportion *q* is confirmed with COVID-19, and a proportion $1 - q$ moves to the suspected and infected class.  Due to the fact that Hubei government amplifies the screening of suspected population greatly after Feb 8th, there will be an extra confirmation rate from the suspected and infected population to confirmed population.  Thus, the confirmation rate for the suspected and infected population in Hubei province is also a piecewise function, which is described as follows.$$\delta_{h}^{Y}(t) = \left\{ \begin{array}{ll}
{\delta_{h}^{Y_{0}},} & {{before}\mspace{6mu}{Feb}\mspace{6mu} 8{th},} \\
{\delta_{h}^{Y_{0}} + \delta_{h}^{Y_{b}}.} & {{ after}\mspace{6mu}{Feb}\mspace{6mu} 8{th}\text{.}} \\ 
\end{array} \right)$$

As the screening policy for regions outside Hubei is stable, the confirmation rate is assumed as a constant $\delta_{c}^{Y}$. Since it is difficult to differentiate whether the suspected population is real infected, it is assumed that the exclusion rate for the suspected but uninfected population (denoted by $\delta_{h}^{S}(t)$ and $\delta_{c}^{S}$) are consistent with those suspected and infected population, that is $\delta_{h}^{S}(t) \equiv \delta_{h}^{Y}(t)$ and $\delta_{c}^{S} = \delta_{c}^{Y}$. The recover rate and disease-induced death rate for the confirmed cases are $\gamma_{h},\gamma_{c}$ and $\alpha_{h},\alpha_{c}$, respectively.  Definitions of all the parameters involved in the model are listed in [Table 1](#t0005){ref-type="table"} .Table 1Parameter and initial values. Parameter orDefinitionsMean valueSourceInitial valuesHubeiOther regions$\beta c_{h}(t)$$(\beta c_{c}(t))$Transmission coefficient----see text$\beta$Transmission probability per contact----see text$c_{h}^{0}$($c_{c}^{0}$)Initial contact rate----see text$c_{h}^{b}$($c_{c}^{b}$)Smallest contact rate----see text$\beta_{h}^{0}$ ($\beta_{c}^{0}$)Initial transmission rate0.71820.9202Estimated$\beta_{h}^{b}$ ($\beta_{c}^{b}$)Smallest transmission rate0.00210.0021Estimated$r_{h}$($r_{c}$)Exponential decreasing rate of contact rate0.16090.2377Estimated$1/\sigma$Latent period4.24.2([@b0070])*m*Suspected rate for the susceptibles5.9857e−075.9857e−07Estimated*q*$\begin{matrix}
 & {{Confirmation}\mspace{6mu}{proportion}\mspace{6mu}{immediately}\mspace{6mu}} \\
 & {{after}\mspace{6mu}{diagnosed}} \\
 & \\
\end{matrix}$0.09280.0928Estimated$\delta_{h}^{I}$($\delta_{c}^{I}$)Diagnosis rate0.20010.2820Estimated$\delta_{h}^{Y}(t)$($\delta_{c}^{Y}$)$\begin{matrix}
 & {{Confirmation}\mspace{6mu}{rate}\mspace{6mu}{from}\mspace{6mu}{suspected}\mspace{6mu}{and}} \\
 & {{infected}\mspace{6mu}{class}} \\
 & \\
\end{matrix}$--0.3005Estimated$\delta_{h}^{Y_{0}}$$\begin{matrix}
 & {{Initial}\mspace{6mu}{confirmation}\mspace{6mu}{rate}\mspace{6mu}{from}\mspace{6mu}{suspected}} \\
 & {{and}\mspace{6mu}{infected}\mspace{6mu}{class}} \\
 & \\
\end{matrix}$0.1391--Estimated$\delta_{h}^{Y_{b}}$$\begin{matrix}
 & {{Extra}\mspace{6mu}{confirmation}\mspace{6mu}{rate}\mspace{6mu}{from}\mspace{6mu}{suspected}} \\
 & {{and}\mspace{6mu}{infected}\mspace{6mu}{class}} \\
 & \\
\end{matrix}$0.3104--Estimated$\delta_{h}^{S}(t)$($\delta_{c}^{S}$)$\begin{matrix}
 & {{Exclusion}\mspace{6mu}{rate}\mspace{6mu}{from}\mspace{6mu}{suspected}\mspace{6mu}{but}} \\
 & {{uninfected}\mspace{6mu}{class}} \\
 & \\
\end{matrix}$$\delta_{h}^{Y}(t)$$\delta_{c}^{Y}$Estimated$\gamma_{h}$($\gamma_{c}$)Recover rate1/11.51/11.5([@b0070])$\alpha_{h}$($\alpha_{c}$)Disease induced death rate0.00488.7240e−04Estimated$r_{i}(t)$move in mobility ratevariablevariablesee text$r_{e}(t)$move out mobility ratevariablevariablesee text$S_{h}(0)$($S_{c}(0)$)Initial Susceptibles$5.9170e + 07$1.3362e + 09([@b0050])$E_{h}(0)$($E_{c}(0)$)Initial latent population84.99920Estimated$I_{h}(0)$($I_{c}(0)$)Initial infectious population1040([@b0020])$Y_{h}^{I}(0)$($Y_{c}^{I}(0)$)Initial suspected and infected population00Data$Y_{h}^{S}(0)$($Y_{c}^{S}(0)$)Initial suspected but uninfected population00Data$Q_{h}(0)$($Q_{c}(0)$)Initial confirmed cases270Data$R_{h}(0)$($R_{c}(0)$)Initial recovered cases00Data

The effective reproduction numbers for Hubei and regions outside Hubei can be calculated without considering the mobility of population between these two patches.  As the model structures for both patches are identical, in the following, only the detailed derivation procedure of the effective reproduction number for Hubei patch is given.  We firstly calculate the basic reproduction number when all the involved parameters are constants by using the next generation matrix method ([@b0095], [@b0015]).  The linearized infection subsystem for Hubei patch is$$\begin{array}{ll}
\frac{dE_{h}}{dt} & {= \beta c_{h}I_{h} - \sigma E_{h},} \\
\frac{dI_{h}}{dt} & {= \sigma E_{h} - \delta_{h}^{I}I_{h},} \\
\frac{dY_{h}^{I}}{dt} & {= \delta_{h}^{I}(1 - q)I_{h} - \delta_{h}^{Y}Y_{h}^{I},} \\
\frac{dQ_{h}}{dt} & {= \delta_{h}^{I}\mathit{qI}_{h} + \delta_{h}^{Y}Y_{h}^{I} - \gamma_{h}Q_{h} - \alpha_{h}Q_{h}\text{.}} \\ 
\end{array}$$

Thus, the transmission and transition matrix can be given as$$F = \begin{pmatrix}
0 & {\beta c_{h}} & 0 & 0 \\
0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 \\
 & & & \\
\end{pmatrix},V = \begin{pmatrix}
\sigma & 0 & 0 & 0 \\
{- \sigma} & \delta_{h}^{I} & 0 & 0 \\
0 & {- \delta_{h}^{I}(1 - q)} & \delta_{h}^{Y} & 0 \\
0 & {- \delta_{h}^{I}q} & {- \delta_{h}^{Y}} & {\gamma_{h} + \alpha_{h}} \\
 & & & \\
\end{pmatrix}\text{.}$$

The basic reproduction number is the dominant eigenvalue of matrix $\mathit{FV}^{- 1}$,$$\rho(\mathit{FV}^{- 1}) = \frac{\beta c_{h}}{\delta_{h}^{I}}\text{.}$$

Thus, the effective reproduction number for Hubei patch and regions outside Hubei patch are$$R_{t} = \left\{ \begin{array}{ll}
{\frac{\beta c_{h}(t)}{\delta_{h}^{I}},} & {{Hubei}\mspace{6mu}{province},} \\
{\frac{\beta c_{c}(t)}{\delta_{c}^{I}},} & {{Regions}\mspace{6mu}{outside}\mspace{6mu}{Hubei}\text{.}} \\ 
\end{array} \right)$$

2.2.  The data {#s0020}
-------------

**COVID-19 reported data.**  The data is collected from the website of National Health Commission of the People's Republic of China([@b0055]) and Health Commission of Hubei Province ([@b0025]), on which the number of daily new/cumulative cases and new/cumulative deaths for both Hubei province and Mainland China are reported every day from Jan 11th, 2020.  By which we can calculate the number of such cases for regions outside Hubei easily.  Since there are no new cases reported from Jan 11th to 15th, the data from Jan 16th are used in this article.  The number of daily reported new confirmed cases and daily new deaths for Hubei and regions outside Hubei from Jan 16th to Feb 21st are described in [Fig.  2](#f0010){ref-type="fig"} (A) and (B).  Meanwhile, the number of existing suspected cases in Hubei province from Feb 8th to Feb 21st are collected, which is illustrated in [Fig.  2](#f0010){ref-type="fig"} (C).  The number of daily new suspected cases can also be obtained for mainland China from Jan 22nd, as shown in [Fig.  2](#f0010){ref-type="fig"} (D).Fig.  2Daily data for Hubei and regions outside Hubei. ( A) new reported cases, (B) new deaths, (C) existing suspected cases in Hubei province, (D) new suspected cases. 

**Mobility data** We collect the mobility index between Wuhan and other regions except for Wuhan from Baidu Migration ([@b0005]), [Fig.  3](#f0015){ref-type="fig"} (A) shows the mobility index between Wuhan and other regions from Jan 1st to Jan 23rd. [ Fig.  3](#f0015){ref-type="fig"}(B) illustrates the proportions of within-province travel among the total outbound and inbound travel of Wuhan. Fig.  3Mobility indices from and to Wuhan city. ( A) the mobility index. ( B) the proportion of mobility within Hubei province. 

2.3.  Parameters and initial values {#s0025}
----------------------------------

We firstly determine some parameter and initial values by references and reported data, then estimate the remaining parameters and initial values by using Least Square Method.  So as to avoid the temporal autocorrelations in the cumulative data ([@b0035]), the daily new reported cases, deaths in Hubei and regions outside Hubei, the existing suspected cases in Hubei province from Feb 8th to Feb 21st, and the daily new suspected cases in mainland China are used to estimate the unknown parameters.  Since most of the infectious cases before Dec 31st, 2019 are linked to Huanan market ([@b0040]), we consider the disease spreading after Dec 31st, 2019.  Thus, the initial value for each class should be the corresponding value on Dec 31st, 2019.  According to the reported data, we can determine the initial values for the suspected class, confirmed class and recovered class for Hubei province.  The Epidemiology Working Group for NCIP Epidemic Response of Chinese CDC ([@b0020]) published that 104 cases are symptomatic before Dec 31st, 2019, that is $I_{h}(0) = 104$. Meanwhile, the initial values for all the classes other than the susceptible class in regions outside Hubei should be very small, thus are assumed to be 0.  The initial values for susceptible classes for Hubei and regions outside Hubei can be regard as the total population in Hubei and regions outside Hubei, respectively, which are obtained from China Statistical Yearbook ([@b0050]).  The mean duration of the latent period is determined according to the estimation by Sanche et al. ([ @b0070]), i.e. $\sigma = 1/4.2$. Sanche et al. ([ @b0070]) also estimated that the duration from initial hospital admittance to discharge is 11.5 days, thus the recovery rates are $\gamma_{h} = 1/11.5,\gamma_{c} = 1/11.5$. As strong correlations exist, it is difficult to estimate the transmission probability per contact and the contact rate simultaneously, thus we estimate the transmission rates $\beta_{h}^{0},\beta_{h}^{b},\beta_{c}^{0},\beta_{c}^{b}$ instead. 

By the mobility data we can calculate the number of travelers between Wuhan and regions outside Hubei ([@b0070]).  The mobility rate for people in latent and infectious class between Hubei and regions outside Hubei is calculated according to the mobility data of Wuhan.  The reason is described as follows: 1, Wuhan is the capital of Hubei province and is also the transportation hub of China, most of the travel between Hubei and regions outside Hubei are in fact the travel between Wuhan and regions outside Hubei; 2, most infected cases are located in Wuhan before Jan 23rd, thus the travel between other cities in Hubei and regions outside Hubei has little effect on the disease transmission from Hubei to regions outside Hubei, especially for the importation of infected cases of regions outside Hubei.  We calculate the outbound and inbound mobility rate at time *t* as $r_{e}(t) = \frac{m_{e}(t)}{N_{1}}$, and $r_{i}(t) = \frac{m_{i}(t)}{N_{2}}$, where $m_{e}(t)$ and $m_{i}(t)$ are the number of persons move out from Wuhan and move into Wuhan from regions outside Hubei, and $N_{1}$ and $N_{2}$ are the total population in Wuhan city and regions outside Hubei, respectively.  As Wuhan and nearby cities in Hubei province are quarantined successively from Jan 23rd, we do not consider the population mobility between Hubei and regions outside Hubei after Jan 23rd.  Considering the fact that almost the same number of susceptibles move into or move out from Hubei every day, thus has little effect on the number of susceptibles in Hubei and regions outside Hubei, we ignore the mobility between susceptibles in Hubei and regions outside Hubei. 

2.4.  Simulations {#s0030}
----------------

We implement simulations to study the effect or work resumption on the epidemic of Hubei and regions outside Hubei, and the effect of time for lifting the lockdown in Hubei on the epidemic of regions outside Hubei.  Work resumption will inevitably lead to the increase of contact rate, but the increase of effective contact rate (i.e. the contact that can lead to infection) could be controlled at a certain low level by preventive measures.  Therefore, we suppose that the contact rate will increase to different levels after work resumption, then investigate the work resumption time on the epidemic of Hubei province.  Similarly, we study the work resumption time on the epidemic of regions outside Hubei. 

To study the effect of lifting the lockdown in Hubei on the epidemic of regions outside Hubei, we determine the mobility rate of the population between Hubei and regions outside Hubei according to the mean mobility rate before the Chinese new year.  As *Chunyun* begins at Jan 10th, 2020, the daily number of people move out from Hubei in the first 14 days after lifting the lockdown is same to the mean daily number of people move into Hubei 14 days before the Jan 23rd, i.e. Jan 10th to Jan 23rd during *Chunyun*.  From the 15th day after lifting the lockdown, the daily number of people move out from Hubei is the same with those move in Hubei before Jan 10th, 2020.  Also, we study the effect of lifting the lockdown in Hubei on the epidemic of regions outside Hubei with different contact rate levels. 

We utilized Matlab R2019 b for windows to implement the parameter estimation as well as the numerical simulations. 

3.  Results {#s0035}
==========

The data fitting results are illustrated in [Fig.  4](#f0020){ref-type="fig"} , where (A), (B) and (C) describe the results for daily new reported cases, deaths and existing suspected cases in Hubei province, (D), (E) show the results for daily new reported cases and deaths in regions outside Hubei, and (F) describe the daily new suspected cases in mainland China.  The estimated parameter values and the initial values of the model are listed in [Table 1](#t0005){ref-type="table"}.  As shown in [Fig.  4](#f0020){ref-type="fig"}, the model fits the real reported data well.  Based on the estimation results, the effective reproduction numbers for Hubei and regions outside Hubei can be calculated, which are described in [Fig.  5](#f0025){ref-type="fig"} .  The basic reproduction number (i.e. the effective reproduction number before Jan 23rd, 2020) for Hubei and regions outside Hubei are 3.59 and 3.26.  The high value of the reproduction number indicates that the COVID-19 have relatively strong infectivity.  Note that the effective reproduction numbers decrease quickly after Jan 23rd with strengthening prevention and control strategies, and they drop below 1 on Jan 31st and Jan 28th for Hubei and regions outside Hubei, respectively.  Our estimation and model fitting results are verified by comparing the predicated cumulative number of cases with the real data, as shown in [Fig.  6](#f0030){ref-type="fig"} .  It is predicated that the epidemic may be almost controlled in mid- March, and the final accumulative number of cases in Hubei and regions outside Hubei will be about 68,500 cases and 13,900 cases.  After verifying our model estimation we can also estimate the epidemic before Jan 23rd 2020 when the novel coronovirus is hardly known and many infections could not be confirmed.  It was estimated that about 818 cases in regions outside Hubei, either in the latent or infected classes, were imported from Hubei province before Jan 23rd, 2020.Fig.  4Data fitting results. ( A) daily new reported cases in Hubei province, (B) daily new deaths in Hubei province, (C) existing suspected cases in Hubei province, (D) daily new reported cases in regions outside Hubei, (E) daily new deaths in regions outside Hubei, (F) daily new suspected cases in mainland China.  The orange stars illustrate the real data, the blue curves show the model predicting results. Fig.  5The effective reproduction number in Hubei and regions outside Hubei. Fig.  6Comparison between model prediction with real reported data. ( A) cumulative number of cases in Hubei province, (B) cumulative number of cases in regions outside Hubei province.  The orange stars illustrate the real data, the blue curves show the model predicting results. 

To investigate the effect of work resumption on epidemic, we initially consider variation in contact rate within each patch (Hubei province or regions outside Hubei as decoupled patch) on the number of daily new infectious cases. [ Fig.  7](#f0035){ref-type="fig"} shows the effect of work resumption in Hubei province, thus the variation in contact rate within Hubei patch. [ Fig.  7](#f0035){ref-type="fig"} (A) illustrates the result of no work resumption, (B), (C) and (D) describe the results if the contact rate becomes 5%, 10%, 20% of the level before Jan 23rd, 2020.  While in each sub-figure, dark, cyan, red, green and blue curve denotes that the work resumption begins on Mar 2nd, Mar 9th, Mar 16th, Mar 23rd and no work resumption, respectively.  From [Fig.  7](#f0035){ref-type="fig"} (B) we can get that if the contact rate becomes 5% of the level before Jan 23rd, all the curves do not show significant increase.  That is, if the contact rate can keep at a relative low level, work resumption has little effect on the epidemic in Hubei province as long as the time of work resumption is later than March 2nd.  If the contact rate becomes 10% of the level before Jan 23rd, as shown in [Fig.  7](#f0035){ref-type="fig"} (C), the dark cure and cyan curves increase quickly and then decrease quickly.  While if the contact rate becomes 20% of the level before Jan 23rd, as shown in [Fig.  7](#f0035){ref-type="fig"} (D), the dark and cyan curves increase significantly to much higher levels compared with that of [Fig.  7](#f0035){ref-type="fig"} (C), and the red and green curves increase slightly.  But all the curves decrease very slowly in [Fig.  7](#f0035){ref-type="fig"} (D). [ Fig.  7](#f0035){ref-type="fig"} (C) and (D) indicate that if the contact rate increase to a higher level after work resumption, there may be another outbreak with both the scale and duration of the outbreak increasing with the contact rate level.  In this situation, work resumption should be postpone to Mar 16th or even more later. Fig.  7Effect of work resumption time on the epidemic in Hubei province. ( A) shows the result of no work resumption. ( B)--(D) show the effect of work resumption.  After work resumption, it is assumed that the contact rate becomes (B) 5%, (C) 10%, (D) 20% of the level before Jan 23, 2020.  The work resumption time are Mar 2, Mar 9, Mar 16 and Mar 23, which are illustrated in black, cyan, red and green curve, respectively. 

Similarly, [Fig.  8](#f0040){ref-type="fig"} describes the effect of work resumption on epidemic in regions outside Hubei, i.e, the variation in contact rate within regions outside Hubei patch. [ Fig.  8](#f0040){ref-type="fig"} (A) illustrates the result of no work resumption, (B), (C) and (D) describe the results if the contact rate becomes 5%, 10%, 20% of the level before Jan 23rd.  All the curves show similar trends to that of [Fig.  7](#f0035){ref-type="fig"}.  But the difference is that the amount of daily new infectious cases is sufficiently small, i.e. less than 2 cases, compared with that of [Fig.  7](#f0035){ref-type="fig"}, whose daily new infectious cases may increase to more than 40 cases.  This indicates that the epidemic in regions outside Hubei has already been controlled at a relatively low level after Mar 2nd, thus work resumption has very little effect on the epidemic in regions outside Hubei if the contact rate could be kept at less than 20% of the level before Jan 23rd. Fig.  8Effect of work resumption time on the epidemic in regions outside Hubei. ( A) shows the result of no work resumption. ( B)--(D) show the effect of work resumption.  After work resumption, it is assumed that the contact rate becomes (B) 5%, (C) 10%, (D) 20% of the level before Jan 23, 2020.  The work resumption time are Mar 2, Mar 9, Mar 16 and Mar 23, which are illustrated in black, cyan, red and green curve, respectively. 

To examine the impact of lifting the lockdown in Hubei province, we consider the mobility of population between Hubei and regions outside Hubei.  At the meantime, the contact rate also increases to different levels after lifting the lockdown in Hubei.  We then check the variation of the new infectious cases with the variation of the beginning time of lifting the lockdown in Hubei.  Lifting the lockdown in Hubei has little effect on the epidemic of Hubei province, as fewer infected individuals may move into Hubei compared with that move out from Hubei.  Here, we mainly focus on the effect of lifting the lockdown in Hubei on the epidemic in regions outside Hubei.  The results are shown in [Fig.  9](#f0045){ref-type="fig"} with (A) describing the results of continuing the lockdown in Hubei and no variation in contact rate, (B), (C) and (D) describing the results if the contact rate after lifting the lockdown in Hubei in both patches increasing to 5%, 10% and 20% of the level before Jan 23rd.  Dark, cyan, red, green and blue curves in each sub-figure denote that lockdown in Hubei is lifted on Mar 2nd, Mar 9th, Mar 16th, Mar 23rd and lockdown continues at Hubei, respectively.  From [Fig.  9](#f0045){ref-type="fig"} (B)-(D) we can get that all the curves will increase slightly, with the peak of curves decreasing with the postpone of lifting the lockdown in Hubei.  Furthermore, the number of new infectious cases increases with the increasing of contact rate after lifting the lockdown in Hubei, which indicates that keeping relatively low contact rate is crucial to avoid a second outbreak of the epidemic in regions outside Hubei after lifting the lockdown in Hubei province.  Comparing [Fig.  9](#f0045){ref-type="fig"} with [Fig.  8](#f0040){ref-type="fig"}, we can get that though lifting the lockdown in Hubei in March do affect the epidemic in regions outside Hubei, showing as more new infectious cases after lifting the lockdown, the daily number of new infections can be controlled at low level if the contact rate could be kept below 20% of the level before Jan 23rd.  If we intuitively choose the daily number of new infectious cases being no more than 5 cases as a criterion, lockdown in Hubei should lift after Mar 9th if the contact rate could be kept at 20% of the level before Jan 23rd.  In other words, if the contact rate can be kept at no more than 20% of the level before Jan 23rd, lift the lockdown in Hubei after Mar 9th will not induce large scale second outbreak in regions outside Hubei. Fig.  9Effect of lifting the lockdown in Hubei on the epidemic in regions outside Hubei. ( A) shows the results of continuing lockdown measure. ( B)--(D) show the results for lifting the lockdown in Hubei.  After lifting the lockdown in Hubei, it is assumed that the contact rate becomes (B) 5%, (C) 10%, (D) 20% of the level before Jan 23, 2020.  The time of lifting the lockdown in Hubei are Mar 2, Mar 9, Mar 16 and Mar 23, which are illustrated in black, cyan, red and green curve, respectively. 

As there are lots of unknown parameters involved in the model, variations of the parameters may affect the trend of the epidemic greatly.  To study the sensitivity of the unknown parameters on the model prediction results, we conduct Latin Hypercube Sampling (LHS) and Partial Rank Correlation Coefficients (PRCCS) analysis.  All the parameter values are sampled from uniform distributions with the ranges given by experience, which are listed in [Table 2](#t0010){ref-type="table"} .  The results of the sensitivity analysis are illustrated in [Fig.  10](#f0050){ref-type="fig"} , where (A) and (B) denote the variations of the base-10 logarithm of the cumulative number of confirmed cases in Hubei and regions outside Hubei, respectively.  From Figure (A) and (B) we can get that the cumulative number of cases both in Hubei and regions outside Hubei will reach the steady stage before early April, and the base-10 logarithm of the final sizes vary greatly from 0 to 6 and 0 to 4 for Hubei and regions outside Hubei, respectively. [ Fig.  10](#f0050){ref-type="fig"} (C) and (D) illustrate the PRCCs for the dependence of cumulative number of confirmed cases in April 1st in Hubei and regions outside Hubei on each parameter.  As shown in [Fig.  10](#f0050){ref-type="fig"} (C), the initial and smallest transmission rates $\beta_{h}^{0},\beta_{h}^{b}$, exponential decreasing rate of the contact rate $r_{h}$, diagnosis rate $\delta_{h}^{I}$ and disease induced death rate $\alpha_{h}$ affect the epidemic in Hubei province greatly, while all the parameters involved in regions outside Hubei patch has relatively little effect on the epidemic in Hubei province.  As shown in [Fig.  10](#f0050){ref-type="fig"} (C), the transmission rates ($\beta_{h}^{0},\beta_{h}^{b}$) and disease related death rate ($\alpha_{h}$) have positive effect, while the decreasing rate of the contact rate ($r_{h}$) and diagnosis rate ($\delta_{h}^{I}$) have negative effect on the epidemic, which indicates the importance of reducing the transmission rates and disease related death rate, and increasing the diagnosis rate for disease control.  Similar results are obtained for the epidemic in regions outside Hubei, as shown in [Fig.  10](#f0050){ref-type="fig"} (D), the initial and smallest transmission rates $\beta_{c}^{0},\beta_{c}^{b}$, exponential decreasing rate of the contact rate $r_{c}$, diagnosis rate $\delta_{c}^{I}$ and disease induced death rate $\alpha_{c}$ within regions outside Hubei patch have great effects on the epidemic.  It needs to be noticed that the initial transmission rate $\beta_{h}^{0}$ and diagnosis rate $\delta_{h}$ in Hubei patch also affect the epidemic in regions outside Hubei greatly.  This happens due to the fact that lots of infected individuals move out from Hubei to regions outside Hubei before Jan 23rd, the epidemic in Hubei before Jan 23rd plays an important role in the disease transmission in regions outside Hubei. Table 2Ranges of the parameters for sensitivity analysis. ParameterRangesParameterRanges*m*\[1e−8,1e−6\]*q*\[0.05, 0.5\]$\beta_{h}^{0}$\[0.6, 1\]$\beta_{c}^{0}$\[0.6, 1\]$\beta_{c}^{0}$\[1e−4, 0.1\]$\beta_{c}^{b}$\[1e−4, 0.1\]$r_{h}$\[0.05, 0.5\]$r_{c}$\[0.05, 0.5\]$\delta_{h}^{Y_{0}}$\[0.1, 0.4\]$\delta_{c}^{Y}$\[0.1,1\]$\delta_{h}^{Y_{b}}$\[0.1, 0.6\]----$\delta_{h}^{I}$\[0.125, 0.6\]$\delta_{c}^{I}$\[0.125, 0.6\]$\alpha_{h}$\[1e−4,1e−2\]$\alpha_{c}$\[1e−4,1e−2\]Fig.  10Sensitivity analysis. ( A) variations of the logarithm cumulative number of confirmed cases in Hubei province. ( B) variations of the logarithm cumulative number of confirmed cases in regions outside Hubei. ( C) Partial Rank Correlation Coefficients (PRCC) for the dependence of cumulative number of cases in April 1st in Hubei province on each parameter. ( D) Partial Rank Correlation Coefficients (PRCC) for the dependence of cumulative number of cases in April 1st in regions outside Hubei on each parameter. 

4.  Discussion {#s0040}
=============

It is essential to examine the early epidemic situation of COVID-19 by using late and relatively accurate reported data.  This can not only provide the evaluation of the effectiveness and timeliness of public health intervention strategies but also predict the transmission risk of population flow due to returning back to work/study.  A two-patchy model, which describes the mobility of population between Hubei and regions outside Hubei in mainland China, is formulated to estimate the epidemic in Hubei province and try to investigate the time of work resumption and lifting the lockdown measure of Hubei province.  Similar model approach for considering population movement between Hubei and regions outside Hubei can be found in former study ([@b0070]), in which the authors investigated the rapid spread of COVID-19 in early stages.  Our model reflects the effect of time varying interventions on the transmission rate, thus may characterize the epidemic more precisely.  The model is parameterized by fitting 6 columns of reported data from Jan 16th to Feb 21st 2020 including daily new reported cases, new deaths and etc for Hubei province and regions outside Hubei.  As the number of daily reported cases begin to include the clinically diagnosed cases on Feb 12nd, 2020, both the number of new and cumulative reported cases have big jumps on that day.  We do not consider this changing in the model in this article, thus it is difficult to fit the number of new/cumulative reported cases both before and after Feb 12nd well.  To predict the epidemic accurately, more efforts are made to fit the epidemic after Feb 12nd well, thus the estimations before Feb 12nd are slightly larger than the reported data.  Especially, the estimated epidemic situation in Hubei before Jan 23rd, 2020 is more serious than the reported data indicated, as shown in [Fig.  4](#f0020){ref-type="fig"} (A) and [Fig.  6](#f0030){ref-type="fig"} (A).  The effective reproduction numbers for both Hubei and regions outside Hubei are very high before Jan 23rd, and decrease quickly after Jan 23rd, and become less than 1 after Jan 31st and Jan 2th.  It is estimated that the reproduction number before Jan 23rd is as high as 3.59 and 3.26 for Hubei and regions outside Hubei, respectively, which follow the ranges obtained by former researchers ([@b0105], [@b0010], [@b0080], [@b0075], [@b0040], [@b0065], [@b0060], [@b0110], [@b0045]).  It is estimated that the epidemic in Hubei and regions outside Hubei may be almost controlled at early- to mid-March, the final size of the epidemic in Hubei and regions outside Hubei are estimated to be about 68,500 cases and 13,900 cases, and about 818 cases in regions outside Hubei are exported from Hubei before Jan 23rd. 

We investigate the impact of the work resumption and lifting the lockdown in Hubei province on the potential second outbreak.  We consider two scenarios: 1) changing contact rates within each patch at different time points, and 2) different time of lifting the lockdown in Hubei with various contact rates.  Main results can be concluded as follows.  In Hubei province, if the strick interventions are implemented after work resumption, i.e. keep the contact rate at relatively low level, work resumption after March 2nd has little effect on the epidemic of Hubei province.  However, if the contact rate increases to higher levels after work resumption, another outbreak may appear, and the scale and duration of the outbreak increases with the contact rate level.  In regions outside Hubei province, work resumption may increase the daily number of new infectious cases but at a relative low level if the intervention measures implemented suitably after work resumption.  Lifting the lockdown in Hubei province in March may cause an increase of the new infectious cases in regions outside Hubei.  But the contact rate after lifting the lockdown plays a vital role on the scale of second outbreak.  If relatively strict interventions are implemented after lifting the lockdown, the daily number of new infectious cases can be controlled to less than 5 cases if lift the lockdown after Mar 9th.  In conclusion, strict interventions for maintaining the contact rate at a relatively low level is critical to avoid a second outbreak of the epidemic after work resumption or lifting the lockdown in Hubei province. 

As many unknown parameters involved in the model, sensitivity analysis is implemented to see the variation of the epidemic curve when parameters vary, which further indicates the importance for reducing the transmission rate and increasing the diagnosis rate on the disease control.  Nevertheless, there are also limitations for this study.  We do not fit the data for daily new cured cases but use the recovery rate estimated by former study ([@b0070]), as the cure rates in different regions and stages, and for different age groups change a lot, it is difficult to fit the daily data for cured cases and deaths simultaneously.  Furthermore, it is assumed that all the infected cases will not transmit the disease after diagnosis, no matter of immediate confirmation with COVID-19 or moving to the suspected class, thus the cure rate in fact has little effect on the disease transmission. 

CRediT authorship contribution statement {#s0045}
========================================

**Xiaodan Sun:** Methodology, Software, Investigation, Formal analysis, Writing - original draft, Writing - review & editing. ** Yanni Xiao:** Conceptualization, Methodology, Supervision, Writing - review & editing. ** Xiangting Ji:** Methodology, Investigation. 

The authors are funded by the 10.13039/501100001809National Natural Science Foundation of China (Grant Nos.:  11701442 (XS), 11631012 (YX)), the 10.13039/501100009996Shaanxi Province Postdoctoral Science Foundation (XS), the 10.13039/501100012226Fundamental Research Funds for the Central Universities (Grant Nos.:  xzy032020028 (XS), xjj2018259 (XS)).  We would like to thank Dr. Hao Wang and Zhe Liu in Baidu Inc. for very helpful suggestions and discussions which greatly improved this manuscript.