key: cord-0961483-63hcyb9e
authors: Liu, Chenzhengyi; Zhao, Jingwei; Liu, Guohang; Gao, Yuanning; Gao, Xiaofeng
title: D(2)EA: Depict the Epidemic Picture of COVID-19
date: 2020-04-07
journal: J Shanghai Jiaotong Univ Sci
DOI: 10.1007/s12204-020-2170-7
sha: 6ea7e4881093cc0587485ef75150a83055942ae1
doc_id: 961483
cord_uid: 63hcyb9e

The outbreak of coronavirus disease 2019 (COVID-19) has aroused a global alert. To release social panic and guide future schedules, this article proposes a novel mathematical model, the Delay Differential Epidemic Analyzer (D(2)EA), to analyze the dynamics of epidemic and forecast its future trends. Based on the traditional Susceptible-Exposed-Infectious-Recovered (SEIR) model, the D(2)EA model innovatively introduces a set of quarantine states and applies both ordinary differential equations and delay differential equations to describe the transition between two states. Potential variations of practical factors are further considered to reveal the true epidemic picture. In the experiment part, we use the D(2)EA model to simulate the epidemic in Hubei Province. Fitting to the collected real data as non-linear optimization, the D(2)EA model forecasts that the accumulated confirmed infected cases in Hubei Province will reach the peak at the end of February and then steady down. We also evaluate the effectiveness of the quarantine measures and schedule the date to reopen Hubei Province.

The city of Wuhan in China has gained global attention due to an ongoing outbreak of viral pneumonia associated with a novel coronavirus -coronavirus disease 2019 . Reported illnesses have ranged from people with little to no symptoms to people being severely ill (fever, cough, and shortness of breath) and dying [1] . While the infection closely links to contact with bush meat from wild or captive sources at the Wuhan seafood market, human-to-human transmission has been discovered [2] , indicating the risk of a much wider spread. The situation could be worse because the outbreak happens to occur during Chunyun, a largescale conventional passenger transport within China for the Chinese Lunar New Year. Confirmed infected cases have also been found out of China through flights connecting Wuhan and international cities. As of February 22, 2020, according to World Health Organization (WHO) [3] , China (including Hong Kong, Macau, and Taiwan) has announced 76 392 accumulated cases of laboratory-confirmed COVID-19 pneumonia (including 20 673 cured cases) and 2 348 deaths. Another 1 402 accumulated confirmed cases are reported in 26 countries, territories, and areas out of China, with 11 death reported from Korea, Japan, Philippines, France, and Iran.

In order to suppress further infection, the Chinese government and health authorities have been dedicated to a range of response measures, such as locating and quarantining close contacts, blocking Wuhan and a number of other severely infected cities, extending the Lunar New Year vacation, and establishing makeshift hospitals. However, it is found that transmission can occur during the incubation period [4] , which can be as long as 12.5 d [5] with no symptoms. As a result, carriers could have spread the coronavirus nationwide and globally without awareness, and the epidemic picture is still worsening on a daily basis. Under this circumstance, an authoritative analysis of the pneumonia infection is in need to guide life and work schedules, as well as to prevent the public from getting panicked at untrustworthy information sources.

This article develops the Delay Differential Epidemic Analyzer (D 2 EA), a mathematical model that depicts the epidemic picture of COVID-19 comprehensively. The D 2 EA begins with a traditional Susceptible-Exposed-Infectious-Recovered (SEIR) epidemic model, where everybody freely contacts with each other. Then we consider the quarantine policies adopted by Chinese health authorities in the later stage and introduce a set of new states (Quarantined) into the model, which makes it capable to represent the true epidemic dynamics. Finally, we consider the potential future variations of three practical factors, namely, the contact rate (during Chunyun return journey), the recovery rate (due to emergence of vaccine or mutation of the virus), and the quarantine rate (associated with dynamics between hospital resource and spread of infection), and analyze how these factors influence the future trends of the epidemic. In the experiment part, we apply the D 2 EA to Hubei Province, the source of the epidemic. The D 2 EA model is fitted to daily epidemic data (infected, recovered, and dead cases) reported by the National Health Commission of the People's Republic of China, the Health Commission of Hubei Province, and the Centers for Dis-ease Control and Prevention (CDC). We also compare the D 2 EA with the traditional SEIR model to illustrate the virtues of quarantine states and reveal the D 2 EA's superiority. Through further sensitivity analysis, we prove that the D 2 EA is not sensitive to some important parameters so that can comprehensively depict the epidemic picture of the COVID-19 pneumonia.

Since initial cases of clustered pneumonia infection in Wuhan were first reported to the WHO on December 31, 2019 [6] , worldwide medical researchers have been closely watching the situation and engaged in scientific analysis from a lot of dimensions. While some of them study the genome and protein characteristics of the COVID-19 for a better understanding [7] [8] [9] [10] [11] , others look into its potential primary reservoirs to prevent future zoonotic disease [12] [13] . There are also research works upon medicine treatment [14] , vaccine development [15] , effectiveness of airport screen detecting [16] , etc.

Among the various works, those focusing on epidemiologic analysis and transmission dynamics are more related to our work. Among them, Zhao et al. [17] accounted for the impact of the variations in disease reporting rate and modeled the epidemic curve of the COVID-19 cases time series. With the collected data, Wu et al. [18] nowcasted and forecasted the potential domestic and international spread of the COVID-19 outbreak. Riou and Althaus [19] looked into the transmission patterns and indicated clues of human-to-human transmission. Read et al. [20] took the flight connection between Wuhan and other cities into consideration and established a transition model to analyze the transmission dynamics of COVID-19. In principle, the above works exploit typical epidemic models and probabilistic methods to analyze the epidemiologic characteristics of pneumonia statistically. However, most of the above works mainly focus on the early epidemic stage and aim to analyze the intrinsic characteristics of COVID-19 itself. Many essential factors such as the effectiveness of quarantine and the impact of the Chunyun return journey are not quantitatively considered or incorrectly estimated.

The SEIR model is frequently referred to in the aforementioned works and other literature. It is a classical mathematical model to analyze the epidemiologic characteristics [21] . By estimating basic reproductive number, incubation period, and other parameters in interest, the SEIR can describe the early stage of epidemics with different seasons, ages, and heterogeneities on virus transmission [22] [23] [24] . However, in most of the real situations, the traditional SEIR model is too simple to forecast future pictures due to its assumption that the individuals freely contact with each other. Some previous works model physical epidemics by adding parameters or compartments into the traditional SEIR and achieve good results. Funk et al. [25] considered infectious of Ebola virus on the dead and split infectious people into those that seek healthcare and those that do not. Wang and Ruan [26] used limited data to simulate the SARS outbreak in Beijing by simplifying the model to a two-compartment suspect-probable model and a single-compartment probable model. Gilberto et al. [27] suggested spatiotemporal models in order to reproduce the time series data and the spread of influenza A. These works go even further and build their models based on reality, but are therefore limited by their specific situations which are not able to be applied to analyze COVID-19.

In this section, we introduce the scheme of the D 2 EA model, including the basic assumptions, the main components of model, and the practical factors the model considers. Briefly speaking, the D 2 EA starts with the traditional SEIR model and integrates itself with newly designed quarantine states and corresponding transition rules. The D 2 EA also accounts for potential variations of three factors, namely, the contact rate, the recovery rate, and the quarantine rate. An overall scheme of the D 2 EA model is shown in Fig. 1 , where, the center part is the state transition graph of the model; the left part is the practical factors considered by the model; the right and bottom parts are the optimization and evaluation process the model takes up. 

The D 2 EA model is established based on the following basic assumptions.

(1) The susceptible population is homogeneous in terms of age group. The assumption of homogeneity can largely facilitate our modeling. Although death cases are more likely to be found among elderly people, the spread of infection has not shown an affinity to the population at a certain age yet. In fact, the confirmed cases cover all age groups, so it is valid to make such an assumption.

(2) Recovered cases are not going to be infected again. Although it is reported that recovered cases can still be infected, we believe such cases can be neglected because the recovered population will surely pay more attention to their health during the epidemic period, and thus are less likely to get exposed.

(3) Birth rate and natural death rate are neglected. According to our survey, the annual birth rate and annual natural death rate of Hubei Province are 1.15% and 0.70%, much smaller than the death rate (3.55%) and cure rate (21.4%) of COVID-19 in Hubei as of February 22, 2020. Furthermore, it can be inferred that the time span of the epidemic will be far less than a year. As a result, the influence of birth rate and natural death rate can be neglected in our model.

(4) A quarantined human will not further infect other people. In our modeling, as long as one transfers to a quarantine state, we assume he/she will not infect other people anymore. This assumption makes the D 2 EA model more tractable, and it is also practical in reality. (5) No case dies in the incubation period. There is no report about death in the incubation period yet, and we do not consider this probability.

(6) The quarantine duration is long enough. In the D 2 EA model, after a quarantine duration of τ days, close contacts will be divided into positive (infected) cases and negative cases. We assume τ is long enough for the judgment of division to be accurate.

As one of the most classic epidemic models, the SEIR models the flows of the population between four states: S (Susceptible), E (Exposed), I (Infectious), and R (Recovered). Since death cases have been reported, we add a new state D (Death) to better depict the epidemic picture and try to predict the increase of death cases. This traditional model is illustrated in Fig. 2 , and the dashed box over which indicates the free spatial flow of the populations.

Unlike previous cases of SARS and MERS, the exposed population of the COVID-19 pneumonia is also infectious. This fact urges us to modify the conventional model to adapt it to the real case, and we introduce two parameters, β I and β E to represent the probabilities of disease transmission per contact with the infectious population I(t) and with the exposed population E(t), respectively. We also account for the fact that a few of the exposed cases could recover on their own, and introduce two recovery rates γ I and γ E for I(t) and E(t), respectively. We also introduce the death rate α I for I(t). The ordinary differential equations of this traditional model are therefore shown by Eqs. (1)-(6). They describe the transition between these states and can be easily derived according to Fig. 2 .

with

However, this traditional-stage model can only describe the early stage of the epidemic, because the SEIR model assumes the individuals freely contact with each other. To make our model more powerful, we need to consider more factors and revise the model.

After the outbreak of the COVID-19 pneumonia, the Chinese government and health authorities have been dedicated to a range of prompt response measures, such as locating and quarantining close contacts, blocking Wuhan and a number of other severely infected cities, extending the Lunar New Year vacation, and establishing makeshift hospitals. Therefore, we believe that the most important factor to consider is the influence of the quarantine measures.

To account for this factor, we add a range of quarantine states into the model and establish the corresponding differential relations. The revised model is illustrated in Fig. 3 , and the new equation set is explained as follows.

To make the explanation clearer, we introduce general symbols X in (t) and X out (t) for state X(X can be S, E, I, R, D, Q S , Q N , Q E , Q P , or Q I ) to represent the total of the population having reached and left state X up to time t respectively. Correspondingly, the derivatives of X in (t) and X out (t) represent the number of humans coming into and leaving state X per unit time. Therefore we have Eq. (7) to express the dynamics of state X. All the equations we derive in the following are based on this basic equation.

There are totally 5 quarantine states in the revised model, i.e., Q S , Q N , Q E , Q P , and Q I . When infectious cases of state I are constantly sent to the hospital and get quarantined at rate δ, they transition from state I to state Q P , which represents quarantined suspected and positive cases. Here "positive" refers to being exactly infected with the COVID-19. Conversely, "negative" refers to suffering from suspectable symptoms yet due to other diseases (e.g., the common cold), with which we define state Q N , quarantined suspected but negative cases. We define η as the proportion of humans exhibiting suspectable symptoms yet due to other diseases among all susceptible population and assume these humans can also get quarantined at rate δ. Both Q N and Q P are suspected cases and need time T to be confirmed. Then Q P transitions to Q I , the state of confirmed infected cases, while Q N stops being quarantined and returns to S. Hence humans leaving state Q N (or Q P ) at time t are those who enter Q N (or Q P ) at time t − T , and we have

Please notice that these equations are delay differential equations with confirming time T as the delay factor. Meanwhile, humans in state Q I are the confirmed cases being treated in hospital and they will be cured at rate γ or die at rate α. Therefore we have Eqs. (10)-(12) for quarantine state Q N , Q P and Q I .

During time dt, δ(t)I(t) cases from infectious population I(t) and δ(t)ηS(t) cases from susceptible population S(t) are quarantined, transitioning to Q P (t) and Q N (t), respectively. At the same time, all the close contacts of them during the last τ days will also be quarantined. These contacts will be quarantined for τ days, where τ is the longest incubation period for the COVID-19. Among the quarantined contacts, those who are not infected with the COVID-19 (Q S ) will return to S and, according to our assumptions, the exposed ones who do have been infected (Q E ) will either become suspected cases (enter Q P ) at rate ε or recover (enter R) at rate γ E . The number of all susceptible humans contacted with one potentially infected person during the last τ days c(t) is estimated by

Therefore, we have the number of quarantined humans entering state Q S , Q P , and Q E during time dt:

Then the differential equations representing dynamics of quarantine state Q S and Q E are

Till now, we have finished designing and quantitatively representing the dynamics of the five new states which take the quarantine measures into consideration. Since we constitute the quarantine states with cases transitioned from the original S, E, I, and R states, the previous traditional equations also require revision. After proper modification, we obtain Eqs. (19)-(29), which are the complete differential equation set for the D 2 EA model. Please notice that the dQ in S (t) dt term is already shown in Eq. (14) .

with

where X can be S, E, I, R, D, Q S , Q N , Q E , Q P , Q I . Since part of the equations are delay differential equations, we let X(t) = X(0) for t < 0 in Eq. (29), which guarantees the initial condition is continuous on the interval [− max{τ, T }, 0]. E(0) = 1 is because we assume there is only one infected human at first. For convenience, we let C(t) be the primitive function of c(t) and we have

To enhance the reliability of the D 2 EA, we consider potential variations of three factors which could bring out great influence, namely, the contact rate, the recovery rate, and the quarantine rate. The variation patterns of these factors and their influences are defined in the followings. Because new data come out every day, in this section we only use the data before February 8, 2020 in Hubei Province as an example to show how to estimate these rates. In our experiment part, we try to use the newest data and reestimate these rates in the same way introduced in this section.

Contact rate refers to the number of contacts with one person per unit time, i.e., r(t). While the nationwide passenger transport during Chunyun and its return journey can temporarily enhance r(t), the city blocking instructions from authority can keep r(t) low.

We use the idea in the kinetic theory of gases to estimate r(t). We regard a person p as a circle with diameter d and others as points. 0.5d is the longest distance the virus can transmit in air. For COVID-19, d = 3 m. If one person (as a point) enters the circle, we regard this person contacts p. Therefore we can estimate the contact rate r(t) of p:

where, u(t) is the ratio of active time (the time duration people stay outside) and total time; n is the average number of humans per area; v is the mean relative velocity of two persons.

As for Hubei Province, we have n = 300 person/km 2 . And reasonably, v = 1 km/h and u(0) = 0.25, which means a person goes outside for 6 h every day and moves 1 km every hour on average. Then by Eq. (32), we have r(0) ≈ 5 person/d.

Equation (32) gives the general idea of estimation for the contact rate. However, influenced by the city blocking policy and Chunyun, the change of r(t) is complex in recent days. We have to take all these into account. In the following, we let r 1 (t) be the contact rate influenced only by the city blocking policy and r 2 (t) be the rate influenced only by Chunyun. First, we quantitatively estimate the variation of r 1 (t):

where, t 1 = 47 d is the index of the day the city of Wuhan is blocked; t 2 is the day when Wuhan is reopened. Initially, r 1 (0) = r(0). Before t 1 , u(t) slowly decreases because the public is more and more wary about the potential epidemic. As u(t) is an average value, it should not decrease too fast, so we let the average active time reduce around 20 s every day. That is the first part in Eq. (33). After t 1 , the contact rate decreases sharply due to powerful authoritative management and control, especially the city blocking policy.

Considering that other cities in Hubei Province are also quickly blocked after t 1 , we use exponential function to estimate the changing of r 1 (t) of Hubei Province, and it exponentially decreases to 0 after t 1 . Although t 2 has not come yet, it is reasonable to deduce r 1 (t) will return to normal also exponentially after t 2 .

By analyzing the data of passenger flows during Chunyun of the last few years, we find passenger flows will reach the peak for Chunyun and its return journey 4 d before and 10 d after Lunar New Year on average, respectively. We use these results to estimate r 2 (t):

where, t a = 43 d and t 4 = 57 d are the index of the day when passenger flows reach the peak for Chunyun and its return journey, respectively. Initially, r 2 (0) = r(0) as well, and we assume that the active time doubles at the peaks, i.e., u(t) = 0.5, and use two Gaussian functions to fit the peaks for Chunyun and its return journey. Finally, we apply Eq. (35) to get our final representation for the variation of contact rate r(t), which is the geometric mean of r 1 (t) and r 2 (t).

In terms of the recovery rate, here we mainly discuss the recovery rate of Q I , i.e., γ(t). As time goes by, medical workers will get more familiar with the COVID-19 in terms of their treatment and care. Therefore γ(t) tends to be growing along the time and approaches a certain saturation value. We know that γ(t)Q I (t) is the number of cured cases per unit time at time t, so γ(t) can be estimated by the number of cured cases during t -th day, where · is the down rounding operator. To let γ(t) be a continuous function of t and make the form as simple as possible, we try to fit a Logistic function with the number of cured cases every day and approximate the result with piecewise linear function to get the estimation of γ(t), as illustrated by

where, t 5 = 14 d is the length of period during which we have yet no therapy for COVID-19; t 6 = 80d is the day when γ(t) reaches saturation according to our estimation. However, two events could greatly influence γ(t). One is the emergence of the vaccine, which can instantly increase γ(t); the other is the mutation of the virus, which can instantly decrease γ(t) and increase β E , β I , α, and α I . We will briefly analyze the concrete influence these events will cause and offer corresponding suggestions to the public and government in Section 3.4.

The quarantine rate δ(t) refers to the speed of people getting quarantined among the infectious population I(t) per unit time. It is associated with the bulk of medical resource, the number of infectious population to be quarantined, and the number of recovered population to leave the hospital. Please notice that δ −1 (t) means the average time from developing symptoms to getting quarantined, which we refer to as the quarantine buffer, we can use everyday quarantine buffer of the cases to estimate δ(t). From official reports, the quarantine buffer in January is around 15 d and around 6 d on February 8. We get the quarantine buffer from the reports and use Segmented Least Squares algorithm to get the piecewise linear estimation for δ(t). The concrete variation pattern of δ(t) is illustrated by

where, t 7 = 54 d, t 8 = 80 d, and t 9 = 120 d are the piecewise points derived by the algorithm; δ(t) = 0 when t < 0 means nobody is quarantined before t = 0.

We can see that δ(t) increases faster since t 7 . This is due to the sharp increase of infected people in the early stage of the epidemic when the medical resource is still enough. At t 8 , however, the medical resource becomes critical and δ(t) starts decreasing. Finally, at t 9 , a balance is established among limited medical resources, newly infected cases, and recovered cases. t 8 and t 9 have not come yet and they are predicted by the data generated by us. In fact, due to the special measures since February 11, δ(t) increases sharply (from official data, quarantine buffer δ −1 (t) ≈ 2 d), which does not match the result we show here. This does not matter because we just give an example here to show the way we do the estimation. We try to use the latest real data collected from official reports in our experiment to ensure accuracy; moreover, we conduct a sensitive analysis for δ(t) in Section 4 to show that our model is not sensitive to this factor.

In this section, we use data from Hubei Province to experiment with the D 2 EA model and analyze the results. Firstly, we set part of the parameters in the D 2 EA by estimating with statistical data, and then we apply non-linear optimization to fit four key parameters that cannot be easily estimated. Our experiment data come from the National Health Commission of the People's Republic of China, the Health Commission of Hubei Province, and the CDC. We collect the number of confirmed cases, recovery, and death before February 20, 2020 as our input. After that, we use the D 2 EA with fitted parameters to make reasonable predictions, evaluate the quarantine measures, and discuss the influence of possible specific medicine (e.g., the vaccine) and mutation of the virus as well.

The D 2 EA is associated with a number of parameters whose values are to be assigned. While some of these values can be estimated with statistical data, others have to be obtained by fitting the model to real data. The whole parameters associated with the model are listed in Table 1 . 

In Table 1 , the parameters with vacant values are to be assigned through model fitting; those with values in parenthesis are variables whose values can change with time, and the assigned values in the parenthesis are the initial or final values; other parameters are constants with values estimated statistically. Concretely, N (0) is the population of Hubei Province; γ I is estimated by the data in preliminary stage of the epidemic; γ E is a small constant since few cases will recover during incubation period; ε −1 is the average incubation period; η is estimated by the number of cases leaving Q N (t) every day; τ is the duration of quarantine; T is the time needed to confirm a suspected case. r(t), γ(t), and δ(t) are estimated by the methods in Section 2.4 using the newest data, especially considering the special measures since February 11, 2020. More concretely, the variations of these factors are shown in Fig. 4 . We see the quarantine rate sharply increases at February 11, 2020 (t = 65 d), the day when clinically diagnosed cases start to be regarded as confirmed cases. 

In the D 2 EA model, as introduced in Table 1 , there are four parameters to fit. Namely, these parameters are:

α, death rate of quarantined confirmed cases Q I ; α I death rate of free infectious population I; β E , probability of transmission per contact with free exposed population E; β I , probability of transmission per contact with free infectious population I.

We explore non-linear optimization libraries in MAT-LAB to fit the D 2 EA, represented by the delay differential equation set, with the real data collected from official reports about the ongoing epidemic. Concretely, we constrain the model to Hubei Province to relieve the influence of the spatial heterogeneity, and we collect the daily updated number of accumulated confirmed infected cases, recovered cases, and death cases in Hubei up to February 20, 2020 as our first-handed data. Our experiment data come from the National Health Commission of the People's Republic of China and the CDC. The detailed fitting results are shown in Table 2 . Since the D 2 EA model considers five different quarantine states, it can better describe real situations with quarantine policies and make full use of the data to get better results. We can see that the confidence interval (CI) listed in Table 2 is small, which means the fitting is successful.

To this end, we have obtained all the necessary parameter values for the model. In the next stage, we will use the model to evaluate the current situation of the epidemic and make forecasts for the future trends of the ongoing epidemic.

The D 2 EA's forecast curves for quarantine infectious population Q I , recovered population R and deaths D are shown in Fig. 5 . The dark red dots along the curves are the up-to-date ground-truth data we collect from the official source mentioned above.

From Fig. 5 , we see that our Q I , R, and D curves fit the reported data well. The origin of the x-axis (time-axis) of our prediction curves in Fig. 5 represents December 8, 2019, i.e., the day when the first infected cases were sent to the hospital. When this article is finished (on February 22, 2020), we are on the 76th day after that and the model's prediction value fits well the current number of reported accumulated confirmed cases of Hubei Province, i.e., 63 454 cases. The same index for China and the globe (excluding China) is 76 392 and 1 402. According to the further prediction, the accumulated infected cases in Hubei will reach its peak at the end of February (i.e., around February 29, 2020), increasing to around 65 000, and then steady down. 

It is clear that the quarantine policies adopted by the Chinese government and health authorities have made a great difference in preventing the epidemic from wildly spreading. A natural question is that, how well have they been doing? To our knowledge, no literature has discussed this topic quantitatively. With the D 2 EA, however, it is easy to do so.

In the first place, the D 2 EA offers us with the epidemic prediction curves under the real current circumstance because we explicitly take the factor of quarantine measures into consideration. At this stage, we can also neglect the role of all quarantine states. In fact, without the quarantine states, the D 2 EA degenerates to the classic SEIR model. To evaluate how effective the quarantine measures are, we compare the S, E, I, and R curves between the D 2 EA and the degenerated SEIR. The results are shown in Fig. 6 .

For each of the curves in Fig. 6 , the left vertical axis is set for susceptible population S; the right vertical axis is set for exposed population E, infectious population I, and recovered population R. Please notice that the right vertical axis of Fig. 6(a) is 10 4 and that of Fig. 6(b) is 10 5 . What it indicates is that the number of exposed and infectious population is only 1/10 of what it could have been, given that the quarantine measures are conducted. In other words, the quarantine measures taken by the Chinese government and health authorities suppress the expansion of the epidemic within 1/10 of what it would be without any control measures. The quarantine measures are indeed very effective.

In Section 2.4, we introduce variation patterns of three practical factors, namely, the contact rate r(t), the recovery rate γ(t), and the quarantine rate δ(t). Section 3.1 depicts the concrete trend estimations with regard to the latest data in Fig. 4 , and we have observations and discussion as follow:

(1) While the Chunyun return journey will inevitably bring larger chances of contact, as long as Hubei Province, especially the city of Wuhan is not reopened too early, the influence of Chunyun can be suppressed. In Section 4.1, we will quantitatively discuss what the most appropriate time is for Hubei Province to reopen. Anyway, life in Hubei Province should come back to normal before May.

(2) Currently, we can expect to see a higher recovery rate of the COVID-19, as shown in Fig. 4(b) , and it will eventually approach 0.07/d if no external events occur. This trend could be greatly influenced by two events, however, namely, the emergence of specific medicine (e.g., the vaccine) and the mutation of virus. We simulate these two situations and the results are shown in Fig. 7 . For the first event, the emergence of specific medicine can enhance the recovery rate γ(t). We assume the medicine is invented at t = 78 d and let γ(t) increase quickly after that. As simulation shown in Fig. 7(a) , the emergence of specific medicine can accelerate the termination of the epidemic and transfer more cases from state Q I to R. For the second event, if the mutation takes place, the probability of transmission (β E and β I ) and death rate (α I and α) will increase. We assume the mutation takes place at t = 78 d and let the probability of transmission and death rate become higher after that. An increase in the number of deaths can be easily foreseen ( Fig. 7(b) ); however, the scale of infection can still be retained as long as the quarantine measures are properly executed (Fig. 7(c) ). As a result, the wise option for the public to take is to stick to the quarantine policies even when the situation just turns about to be better. Meanwhile, the government should weigh more funds to vaccine production in order to subside this epidemic as early as possible and leave lower expectations for the mutation to take place.

(3) Due to limited medical resources, it is impossible for each confirmed case to get quarantined instantly. The smaller the quarantine rate, the longer the epidemic will last. In fact, according to our estimation in Fig. 4(c) , δ(t) approaching 0.5/d reflects a two-day time delay before an infection can actually be quarantined. Considering the severe shortage of medical resources in Hubei Province, especially in Wuhan, although we have forecasted the number of confirmed cases will steady down soon, the overall situation is not totally optimistic.

Time to reopen Hubei Province t 2 is very important in the whole model. We assume Wuhan is the last city to be reopened and t 2 is just the time to reopen Wuhan. If the authority stops blocking Wuhan too early (i.e., t 2 is too small), the epidemic will be likely to break out again. If it is too late, the economic loss will be huge. Hence the government should carefully choose t 2 to minimize the loss, and t 2 will greatly influence the results of our model. To find out the best value of t 2 , we change the value of it and conduct several simulations. The result is shown in Fig. 8(a) . We can see the best time to reopen the province is at least 100 d after the outbreak. That is to say, the authority should reopen Hubei Province after March 16, 2020. For prudential reasons, we recommend that Hubei, especially the city of Wuhan, should be reopened no earlier than in late March, or the epidemic may break out again with a high chance.

Quarantine rate δ(t) is associated with the population size of Q S , Q E , Q N , Q P , and Q I . If the quarantine rate is high, it means the government can quickly quarantine infected cases and thus more infection can be avoided. However, the quarantine rate cannot grow very high because the medical resource is limited. We believe δ(t) will finally reach an equilibrium point δ(+∞). This point is decided by the speed of resource being sent to the epidemic area, resource being released by humans having finished quarantine, and the demands of newly confirmed infected cases. However, due to the special measures on February 11, 2020, the changing of δ(+∞) is difficult to estimate accurately. In our previous experiment, we assume δ(t) reaches δ(+∞) after February 12, 2020 and set δ(+∞) = 0.5/d, which is derived by the quarantine buffer of 2 d but is likely to be inaccurate. To analyze the sensitivity of the model towards different values of δ(+∞) we change the expression of δ(+∞) and conduct corresponding simulations. The results are shown in Fig. 8(b) , where δ(+∞) = 5/d means the quarantine buffer is 0.2 d. We can see that the model is not sensitive to δ(+∞) as the accumulated number of Q I does not change significantly.

In this article, we introduce our model, the D 2 EA (Delay Differential Epidemic Analyser) to depict the epidemic picture of the ongoing COVID-19 pneumonia. The D 2 EA is revised from the SEIR model, with ample consideration into the quarantine measures and variation of practical factors. Through comparison and analysis, we see that the D 2 EA model depicts a com-prehensive picture of the ongoing epidemic, and thus its forecast is reliable.

In our experiment part, we fit the D 2 EA to the collected data from official reports of the epidemic in Hubei Province, the source of the epidemic, using nonlinear optimization methods and conduct forecast of the epidemic. According to the D 2 EA's forecast, accumulated confirmed infected cases in Hubei will reach the peak at the end of February and then steady down. We also quantify that the currently adopted quarantine measures keep the epidemic 1/10 the scale it could have been, and recommend that Hubei Province, especially the city of Wuhan relieves the lockdown state no earlier than in late March.

novel coronavirus symptoms

A familial cluster of pneumonia associated with the 2019 novel coronavirus indicating person-to-person transmission: a study of a family cluster

WHO. Coronavirus disease 2019 (COVID-19): Situation report-24

Transmission of 2019-ncov infection from an asymptomatic contact in Germany

Early transmission dynamics in Wuhan, China, of novel coronavirus-infected pneumonia

Genomic and protein structure modelling analysis depicts the origin and infectivity of 2019-nCoV

Functional assessment of cell entry and receptor usage for SARS-CoV-2 and other lineage B betacoronaviruses

Detection of 2019 novel coronavirus (2019-nCoV) by real-time RT-PCR

Highly distinguished amino acid sequences of 2019-nCoV

Nucleotide analogues as inhibitors of viral polymerases

Host and infectivity prediction of Wuhan 2019 novel coronavirus using deep learning algorithm

A pneumonia outbreak associated with a new coronavirus of probable bat origin

Drug treatment options for the 2019-new coronavirus (2019-nCoV

Insights into cross-species evolution of novel human coronavirus 2019-nCoV and defining immune determinants for vaccine development

Effectiveness of airport screening at detecting travellers infected with 2019-nCoV

Preliminary estimation of the basic reproduction number of novel coronavirus (2019-nCoV) in China

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

Pattern of early human-tohuman transmission of Wuhan

Novel coronavirus 2019-nCoV: Early estimation of epidemiological parameters and epidemic predictions

Mathematical modeling and research of infectious disease dynamics

Mathematical tools for understanding infectious disease dynamics

The mathematics of infectious diseases

Modeling infectious diseases in humans and animals

The impact of control strategies and behavioural changes on the elimination of Ebola from Lofa County

Simulating the SARS outbreak in Beijing with limited data

Modeling the epidemic waves of AH1N1/09 influenza around the world