key: cord-0859964-9xhtlo1h
authors: Liu, Jijun; Wang, Liyan; Zhang, Qiang; Yau, Shing Tung
title: The dynamical model for COVID-19 with asymptotic analysis and numerical implementations
date: 2020-08-08
journal: Appl Math Model
DOI: 10.1016/j.apm.2020.07.057
sha: 07bb2be25f332680da2afe909df998ed232c20bc
doc_id: 859964
cord_uid: 9xhtlo1h

The 2019 novel coronavirus (COVID-19) emerged at the end of 2019 has a great impact on China and all over the world. The transmission mechanism of COVID-19 is still unclear. Except for the initial status and the imported cases, the isolation measures and the medical treatments of the infected patients have essential influences on the spread of COVID-19. In this paper, we establish a mathematical model for COVID-19 transmission involving the interactive effect of various factors for the infected people, including imported cases, isolating rate, diagnostic rate, recovery rate and also the mortality rate. Under the assumption that the random incubation period, the cure period and the diagnosis period are subject to the Weibull distribution, the quantity of daily existing infected people is finally governed by a linear integral-differential equation with convolution kernel. Based on the asymptotic behavior and the quantitative analysis on the model, we rigorously prove that, for limited external input patients, both the quantity of infected patients and its variation ratio will finally tend to zero, if the infected patients are sufficiently isolated or the infection rate is small enough. Finally, numerical performances for the proposed model as well as the comparisons between our simulations and the clinical data of the city Wuhan and Italy are demonstrated, showing the validity of our model with suitably specified model parameters.

The spread of COVID-19 all over the world from the end of 2019 was a critical incident. By strong interventions and powerful treatments, the spread of COVID-19 in China has been effectively controlled, and most patients have been successfully cured. However, after an outbreak in South Korea, Japan and other Asian countries, the epidemic has now spread rapidly and severely to European and American countries. It was reported that on March 20, 2020, there were 5986 newly confirmed cases and 47,021 cumulative confirmed cases of COVID-19 across Italy, to say noting of 4032 died patients who are positive-tested for COVID-19. The severe situations of COVID-19 have raised universal concerns.

As a new type of epidemic disease, the transmission mechanism of COVID-19 is very complex, without any confirmed conclusion. Due to its huge impact on the whole society, in addition to the active researches on pathogenesis and effective treatment drugs by medical communities, the establishment of corresponding transmission model of COVID-19 based on the fundamental epidemiological law is another important way to the understanding of this new respiratory disease. Moreover, the effective transmission mechanism considering the interactive effects of different factors can validly predict the spread of COVID-19, and consequently help people to prevent the epidemic spread and to reduce its harm to human society. In general, as a new but typical epidemic disease, the spread of COVID-19 depends on the following factors:

• The harmful nature of the viruses themselves, such as the incubation period, infection rate and mortality rate;

• The cognition and therapeutic effects of contemporary medicine on viruses, including the cure rate, diagnosis rate and cure time of patients;

• Measures to control the viruses spread, such as travel restrictions, isolating patients and controlling the imported cases.

For COVID-19, there is little confirmed information about the nature of the virus and the therapeutic effect. It is generally believed that the COVID-19 may have some potential relation with SARS in 2003 [5] , so there have been some research using SARS-related data (e.g. latent period of disease, infection rate) to study the transmission law of COVID-19 [1, 2, 4] . On the other hand, the strict and efficient control measures by the Chinese government since January 2020 (e.g., banning the flow of people in infected areas, treating all the suspected cases, and controlling the external input cases) have had a clear effect on blocking the spread of the disease. Since March 2020, the number of COVID-19 patients in mainland China has decreased significantly. The officially issued information about the infected people can partially describe the propagation property of the virus and help us to evaluate the impact of various control measures from a data perspective [3] . However, the issued data are very limited, which are not enough for giving an effective portrayal of virus transmission in principle.

There have been some classic models on the diffusion of general epidemic diseases, such as Logistic model, SIS model, SIR model, SEIR model, etc [7, 9] , and also the time delay models for diseases transmissions [6, 8] . However, these models take into account neither the human mobility nor the lag effect of the virus incubation, and consequently are not completely applicable to reveal the COVID-19 transmission. Mathematically, the delayed effect of disease transmissions can be represented by some integral terms with convolution kernels. On the other hand, the obervable existing epidemic data can compensate the lack of virus transmission mechanism in some sense. So by modeling the transmission process mathematically embedding the issued outbreak data, we can not only depict the transmission mechanism of the virus itself, but also make use of the mathematical tools to give a forecast on the epidemic trends, and to propose some suggestions on the implementable control measures that should be taken in the future.

This paper considers the combined effects of the crucial factors of COVID-19 such as external input patients, isolation of infected patients, incubation of disease, cure time, mortality rate, and establishes a dynamic model for the transmission process of the epidemic in terms of the quantity of daily existing infected people. Under the assumption that the incubation period, cure time, death time obey the Weibull distribution, the dynamic influence of untreated patients and isolated patients on the epidemic was considered, and the binding relationship between the daily existing cases, the daily new confirmed cases and the daily cured cases was analyzed, and the influence of various factors was finally expressed in terms of the number of daily existing infected cases. Then, an integro-differential equation with time-delay integral kernels is established to govern the distribution of daily existing infection people. Based on the rigorous mathematical analysis, the distribution of the initial outbreak, the number of external imported patients, social control measures and the other quantitative effects on the epidemic infection are involved in our model comprehensively. For the limited imported patients, we obtain the following conclusions. 1. Under sufficient isolation of infected patients or lower transmission rate of infected patients, the number of infections will eventually disappear. This conclusion shows that the effective isolation of the infected people plays a decisive role in blocking the spread of the epidemic, eliminating the possibility that the disease can coexist with human beings for a long time; 2. By strictly controlling the input cases such that the imported cases are at a lower level of growth, and maintaining a small infection rate or sufficient isolation, the variation of daily existing infected people will eventually come to zero, indicating that the spread of the epidemic will be relatively stable after a long time; 3. The numerical simulations of the model are given in different configurations, and the simulations are also compared with the actual situations in Wuhan and Italy by checking the number of daily confirmed cases. The validity of the model is verified for the appropriate specified model parameters.

Suppose that COVID-19 in some city was caused by the infected patients at some initial time (e.g., December 1, 2019), and then all the infected patients are generated by the external input patients and the spread of the local patients, which intensified the epidemic. Denote by I(t) the cumulative infected people at time interval [0, t], which comprises of three ingredients: infected but cured patients I c (t), infected but dead patients I d (t), patients I s (t) who are still infected at time t. Introduce the following quantities:

• i(t): newly infected people due to internal infection at time t (suppose that they just get infected);

• s(t): imported cases at time t (also suppose that they just get infected);

• c(t): cured people at time t, including people cured in hospital c 1 (t), self-cured patients outside the hospital c 2 (t);

• d(t): dead people at time t, including people dead after treatments in hospital d 1 (t), untreated patients died of deterioration d 2 (t).

Then the fundamental relations governing the transmission of COVID-19 are as follows: (2.4) and the initial status I s (0) = I 0 s . (2.5)

From this system, the cumulative cured people and the dead people in hospital in time interval [0, t] are

respectively. In practical situations, we consider the discrete variable t = k(k = 1, 2, · · · ) and then dIs(t) dt = I s (t + 1) − I s (t). In the transmission process, some suspected cases are isolated in the medical institutes, part of them become confirmed cases after the latent period, then these persons get further medical treatments. Notice, the latent period can be considered as 0, if the patients went to hospital directly. In a given region, the officially issued outbreak data include:

• Daily confirmed cases in hospitals j(t);

• Daily cured cases in hospitals c 1 (t) and cumulative cured cases in hospitals I 1 c (t) (these two quantities are determined each other);

• Daily dead cases in hospitals d 1 (t) and cumulative dead cases in hospitals I 1 d (t) (these two quantities are determined each other).

The modelling process by the above general rules should reflect the evolution law for (i(t), c(t), d(t)) and I s (t), which are subject to the constraints each others. Due to the complexity of disease transmissions and medical treatments, we have to establish the evolution process under some assumptions, with the basic consideration that we express the individual differences of patients in latent period, cure period and other indicators by the distributions of some random variables. To this end, we make the following assumptions:

1. The asymptomatic and the symptomatic in the incubation have the same infectivity. Denote by β > 0 the infection rate of the disease, i.e., the number of newly infected people from one patient at unit time.

2. The patients are no longer of infectivity once they are isolated or treated in the hospitals. Suppose that the isolation ratio at time t for all the infected people is κ(t) ∈ (0, 1), which reflects the governmental control capability for the infected people.

Based on these two assumptions, the number of isolated patients at time t is κ(t) × I s (t), and the number of exposed patients is e(t) := (1 − κ(t)) × I s (t). Only the infected people e(t) exposed in the society are infectious, except the daily imported patients s(t). Then the number of newly infected people caused by the exposed patients is

(2.7)

Note that newly infected people at time t are at the beginning of incubation period. To describe the daily cured patients c(t), we have to understand the treatment process after isolation and diagnosis. Denote by j(t) the number of newly confirmed patients at time t, i.e., the daily confirmed cases in the officially released news.

The incubation for each infected patient is different. Introduce the average incubation period τ 1 for all the infected people, and the probability density function for the incubation period τ is h s (τ ; τ 1 ). Further, we assume that the incubation period of all individual cases in the newly infected population is independent and identically distributed. Then h s (τ ; τ 1 )dτ is the proportion of patients who are of τ incubation period. Obviously, the density function satisfies

To describe the infection propagation by the infected people at initial time t = 0, namely, I s (0) = I 0 s , we consider it equivalently as a continuous and uniform input s 0 (t) in the time interval [0, 0 ] for some small 0 > 0. Then the infection source in the transmission process is

which satisfies 0 0 s 0 (t)dt = I 0 s . In the actual discrete model, we can take 0 = 1 or 0 = 2, 3. Since the isolation ratio is κ(t), then for all the newly infected people i(t) + S(t) at time t, the number of isolated patients is κ(t)(i(t) + S(t)). These isolated patients will get confirmed once they have passed the incubation period. For any random τ ∈ [0, t] with probability density function h s (τ ; τ 1 ), the isolated part of newly infected people at time t − τ , i.e., κ(t − τ )(i(t − τ ) + S(t − τ )), will get confirmed at time t. Owing to the evolution law, κ(t )(i(t ) + S(t )) for t > t makes no contribution to the status at time t. So, for the transmission process initialized from t = 0, we have

In the case that t = K is discrete positive integer variable, it becomes

indicating the probability of patients with incubation k in newly infected people i(K − k) + S(K − k) on the day K − k. Thus, among all these isolated patients, p(k)κ(K − k)(βe(K − k) + S(K − k)) will get confirmed on the day K.

Note that the newly confirmed people j(t) at time t is different from newly infected people I s (t + 1) − I s (t). On the other hand, there are two parts for the newly cured patients: one comes from the isolated people, i.e., c 1 (t), and the other comes from the exposed people, i.e., c 2 (t). The former reflects the treatment ability from the medical communities, while the latter depends on the resistance from the patients themselves. For the former, there are two indices to describe the treatment effect, namely, the average cure time τ 2 after diagnosis and the recovery rate κ 2 . The cure time τ with average τ 2 is a random variable, and we suppose the density function of a patient's cure time is h iso c (τ, τ 2 ). Similarly to the analysis on j(t), the newly confirmed patients j(t − τ ) at time t − τ will be cured at time t. Then the number of cured patients from isolation is

For the latter, assume that the recovery time of exposed patients among newly infected patients is a random variable with density function h exp c (τ ; τ 4 ), and the recovery rate is κ 4 . Then we also have

Thus the number of newly cured patients at time t is

Just like c(t), the newly dead patients can also be separated into two parts: d 1 (t) from the isolated patients and d 2 (t) from the exposed patients. The former also reflects the treatment level, and the latter reflects the patients' own resistance. We denote the average time between diagnosis and final death of isolated patients by τ 3 , with density function h iso d (τ ; τ 3 ); while the average time to death of the exposed people is τ 5 , with density function h exp d (τ, τ 5 ). Therefore it follows that

Obviously, there should be the relations κ 2 > κ 4 , τ 2 < τ 4 , τ 3 > τ 5 , which means that the mortality rate decreases and the average survival time increases due to the medical treatments.

Based on the above analysis and assumptions, now we have established the following coupled system for COVID-19 transmission in t > 0:

This model is essentially an integro-differential equation with convolution kernels. Once we solve I s (t) from this system, (I c (t), I d (t), I(t)) can be obtained easily. For the above model, we can also describe the number of patients recovered and died based on another assumption. The number of patients exposed at time t is e(t), which are infectious, causing newly infected patients i(t). The other part of the newly infected are the imported patients S(t). For all these newly infected patients, the number of isolated patients is κ(t)(i(t) + S(t)). Denote by the recovery rate κ 8 and the time τ from infection to cure with the density function h iso 8 (τ ; τ 8 ) of average time τ 8 . Meanwhile, for the exposed part (1 − κ(t))(i(t) + S(t)) of newly infected patients, we introduce the recovery time τ with density function h exp c (τ ; τ 4 ) for average recovery time τ 4 , and the recovery rate κ 4 . Then the number of newly cured patients is

Similarly, for d(t),we introduce the density function h iso 10 (τ ; τ 10 ), then we have

Notice, the second terms on the right hand of (2.17),(2.18) and the second terms of c(t) and d(t) in (2.16) are the same. In other words, these two different considerations, with different density functions (one is after the incubation period, the other is just after the infection), differ only in the terms for the isolated part. On the other hand, if there is no incubation period for the newly infected patients i(t) + S(t), the isolated patients κ(t)(i(t) + S(t)) will be confirmed immediately, yielding h s (t − τ ; τ 1 ) = δ(t − τ ). Then we have

That is, (2.17) and (2.18) can be considered as the simplification of c(t) and d(t) in (2.16). For all the probability density functions of random variable τ , i.e., h s (τ ;

with mean value µ > 0, we can also consider the Weibull distribution, which has been considered as a standard distribution in the epidemiological investigation:

where λ > 0 is the shape parameter. It is the exponential distribution for k = 1, and the Rayleigh distribution for k = 2. By straightforward calculations, we have

for different (k, λ), leading to various density functions in our model. The mathematical models and the numerical simulations with other density functions are the same, but the corresponding theoretical analysis needs to be carried out separately.

The coupled system (2.16) is an integro-differential equation system with convolution kernels, which reveals the nonlocal property of the epidemic spread. First, we need to analyze the qualitative properties of the solution. Owing to the complicated coupling relationship among the components of (I s (t), i(t), j(t), c(t), d(t)), we eliminate the last four components. By direct computations and exchanging the order of integrations, the governed equation for

together with the initial value

Obviously, for given density functions as well as the related parameters and the imported cases s(t), (2.21)-(2.22) is an integro-differential equation, which can be further transformed into a linear Volterra integral equation of the second kind with continuous kernel, thus it is well-posed, i.e., I s (t) can be uniquely solved. Then for known I s (t), (i(t), j(t), c(t), d(t)) can be determined from the other equations in (2.16), which implies that all the relevant quantities describing the spread of the COVID-19 can be obtained.

The last item in the right hand side of (2.21) presents the impact of spontaneous healing and natural death of the exposed people on the newly infected patients, which is an important indicator to describe the development of the disease without intervention. This part of the cured and dead patients are no longer infectious, which reduces the growth rate of the infections. However, such a positive impact is not due to the medical treatments. On the other hand, the second term on the right side of (2.21) indicates the effect from the isolated patients, those who died of invalid treatments are no longer infectious, which also reduce the growth rate of the infections. This term reveals the influence of medical intervention. Therefore, the ratios κ 2 and κ 4 ∈ (0, 1) together with (τ 2 , τ 3 , τ 4 , τ 5 ) are the fundamental values characterizing the spread of the epidemic.

For given infectious diseases, the epidemic situation should tend to some stable status after a long time. Any reasonable mathematical model describing the spread of COVID-19 should meet this requirement. The following result ensures that COVID-19 transmission governed by our proposed model will finally disappear, provided that some a-priori reasonable assumptions on the isolation ratio, imported cases are specified.

Theorem 1. For given T > 0, there exists a unique solution I s (t) ∈ C[0, T ] to (2.21)-(2.22). Assume that the corresponding random variables describing the disease property (incubation period, cure period, death period, etc) obey the Weibull distribution f k λ (t), and the quantity of the imported cases satisfies the growth condition i.e., for a closed system with limited imported cases, the dynamical system for COVID-19 tends to be stable after a long time, and all the infected patients will eventually disappear.

The biological explanation of the two conditions in this result is clear: the limited external input condition requires that the number of the daily imported patients decay exponentially (it is obviously reasonable for practical transmissions, e.g., s(t) is of finite support); while the condition 0 < β(1 − κ(z)) 1 requires lower transmission rate or large isolation ratio.

Proof. Define 

Since h j (t) is of the form of Weibull function (2.19) with τ j = Γ(1 + k −1 )λ j (j = 4, 5, 8, 10), we have from direct computations that

which lead to

. Since the support of s 0 (t) is [0, 0 ] and 0 0 s 0 (z)dz = I 0 s , then (3.3), (3.5) and (3.6) lead to for t > 0 that

Now we estimate I s (t) for large t for two Weibull distributions with k = 1, 2. Notice, we have I s (t) ≥ 0 by ( 

By variable translation t − 0 = t , this estimate says

Consider the general form of (3.11): 

So the induction arguments yield

) n−1 λ 0 ) 2 dτ, n = 1, 2, · · · (3.14)

Inserting (3.14) into (3.13), we have for βκ 0 λ 0 ∈ (0, 1) that

So, for any given t ∈ [0, ∞), we have by taking N → ∞ that

Hence, (3.11) for k = 2 leads to the estimate

According to the definition of F (t) in (3.10), we have for t > 0 that

where E 0 = max{ 0 , T * } > 0, C * = C * ( 0 , T * , I 0 s , E 0 ). Then we obtain

We assert that the series in the right hand side tends to zero as t → ∞. Without loss of generality, we consider the case λ 0 = 1. By the estimation

Therefore the integral term in (3.16) tends to 0 as t → ∞ from the definition of limitation. Using (3.17) again, we have lim t→∞ I s (t) = 0. Since we have clarified the asymptotic behavior of I s (t), the asymptotic behavior of (i(t), c(t), d(t)) can be shown immediately. The proof is complete.

An indicator to describe the epidemic is the turning point t * of I s (t), thus we need to consider the zero point of I s (t). Since I s (t) > 0 at the initial stage of the outbreak, once the zero point of I s (t) appears, the number of daily infected patients will start to decrease (although may not continuously decrease after this turning point t * ). Define

Then the straightforward computations of (3.3) say

Proof. Under the condition lim t→∞ G(t) = G 0 , it can be verified directly that

Applying this limit in the identity Now let us show, by our proposed model, I s (t) always has a zero point for bounded imported patients, i.e., the COVID-19 must have a turning point for the number of infected people (not the number of daily added infected ones). Moreover, I s (t) is very small for large t > 0, which means that the variation of the number of the infected people will be very smooth after large time.

Theorem 3. Assume that the isolation ratio satisfies lim t→∞ κ(t) = κ ∞ ∈ (0, 1), and the imported cases s(t), which is continuous on (0, ∞), meets the growth condition of Theorem 1. Then we have 1. lim t→∞ I s (t) = 0; 2. There exists some zero pointt * > 0 of I s (t) on (0, ∞).

Next we prove that there exists some zero point of I s (t) on [0, ∞). If I 0 s = 0, there must be 0 ≤ s(t) ≡ 0, which ensures I s (t 0 ) > 0 for some t 0 > 0. Hence, the epidemic spread can be considered with the initial instant t 0 . So we can assume I 0 s > 0 without loss of generality. Since s(t) and I s (t) are uniformly bounded continuous functions on [0, ∞), we have I s (t) ∈ C[0, ∞) and

We assert that there must exist some t * ∈ (0, ∞) such that I s (t * ) < 0. Otherwise, I s (t) ≥ 0 on [0, ∞). Then it holds for all t > 0 that

However, by Theorem 1, for β max [0,∞) (1 − κ(z)) > 0 small enough, we have lim t→∞ I s (t) = 0, which is a contradiction. Hence, by (3.22) , there exists somet * ∈ (0, t * ) such that I s (t * ) = 0. The proof is complete.

Remark 4. This result asserts two things. Firstly, there will be a turning point in the spread of the epidemic. Secondly, after a long time, the change of I s (t) (i.e., I s (t)) will be sufficiently small. However, we cannot assert that I s (t) has only one zero point on (0, ∞), that is, the number of infected patients may not monotonously decreasing. The reason is that we do not specify too many restrictions on the external input function s(t).

Remark 5. According to Theorem 1 and Theorem 3, the daily existing infected patients will eventually disappear completely, in the sense that This conclusion is based on the following two premises: 1. There are strict restrictions on the daily imported cases s(t); 2. The infection rate is sufficiently small or the isolation rate is sufficiently close to 1. Because the patient's infection rate β is not clear at present stage, The efficient controlling of the input of s(t) or improving the isolation rate κ(t) is the crucial measure to eliminate the disease finally.

By the analysis in above sections, we have described the transmission of the epidemic spread by a self-enclosed system with respect to the intermediate function I s (t). If all the parameters in this model and inputs are known, the system can completely determine I s (t). The model parameters and the input sources can be divided into the following five categories:

• sources leading to the spread of the disease: the number of initial patients I 0 s and the imported patients s(t);

• the qualitative hypothesis of disease: the specific form of the density functions;

• the quantitative description of the epidemic disease: τ j for j = 1, 2, 3, 4, 5 in the density functions;

• the quantitative parameters of diseases: κ 2 , κ 4 , β;

• social ability to control the spread of the epidemic: the isolation rate κ(t).

This section analyzes the trend of COVID-19 by numerically studying the performances of I s (t) with different types of inputs for our model. Obviously, the behavior of I s (t) depends on the density functions for the epidemic situations and also the parameters in these functions. The determination of these unknown parameters in governed system by the observed data is a nonlinear inverse problem for differential equation system, and the further prediction of future epidemic trend based on these obtained parameters is a kind of initial value problem for a non-local integro-differential equation.

Based on the given model parameters, now we verify the reasonability of our proposed model. The task making use of the observed data related to the transmission of disease to determine these model parameters should be the further work in the future. The Weibull distribution and Gaussian distribution of the density function for our model are used to simulate I s (t) numerically for the external input s(t) in different forms. Moreover, the newly confirmed cases j(t) and the corresponding cumulative confirmed cases J(t) on [0, t] are also computed. By comparing the numerical results for (j(t), J(t)) from our proposed model with the officially issued data of Wuhan from January 23 to March 18, 2020 and also of Italy, the simulant data are in good agreement with the actual data, when the model parameters and the initial data are properly specified.

Firstly, we consider the numerical behavior of I s (t) for different forms of s(t), with the initial data I s (0) = I 0 s = 100. For the Weibull distribution functions f k λ (t) with k = 1, 2, we take the parameters in the model as β = 0.2, κ(t) = 0.85, κ 2 = 0.9, κ 4 = 0.6, τ 1 = 5.2, τ 2 = 9, τ 3 = 15, τ 4 = 18, τ 5 = 10.

Consider s(t) in the following three forms: For specified density function f 1 λ (t), we show the numerical behaviors corresponding to s(t) in three different forms specified by (4.1) in Fig.1, Fig.2 and Fig.3 , respectively. In Fig.1 , it is easy to find that the resulting I s (t) is continuous but with a sudden sharp point due to the differential equation in (2.16) for discontinuous s(t). However, due to the smoothing process of the convolution, even for discontinuous input s(t), (j(t), c(t), d(t)) is still smooth. On the other hand, with the convolution kernel, our model is a good illustration of the non-local effect of COVID-19 transmission: the spread of the epidemic is a trailing process describing slow diffusion phenomenon, similarly to the slow diffusion depicted by time fractional derivative.

The numerical performances in Fig.1 remind us that the concrete form of s(t) may have some impact on I s (t), and will change the spread of epidemic dramatically. To show this phenomenon clearly, consider the following two distributions of s(t): For these two configurations with the same total input ∞ 0 s(t)dt = 3300, s(t) always has some periodic structure in finite time interval, simulating the periodic variation of the imported cases. The numerical behaviors of (I s (t), j(t), c(t), d(t)) corresponding to s(t) in (4.2) and the density function f 1 λ (t) are shown in Fig.4 and Fig.5 . It is found that the monotone behaviors of I s (t) on two sides of the peak point as shown in Fig.1-Fig.3 do not exist anymore, i.e., the spread of epidemic has some rebounding phenomenon. On the other hand, we also find that the corresponding (I s (t), j(t)) are similar either in the distributions or in their peak values for these two configurations. This numerical performance illustrates the necessity of controlling the imported patients.

In the case of Rayleigh distribution with k = 2, we obtain similar results as shown in Fig.6 for s(t) in the second form of (4.2). The difference is that (I s (t),j(t),c(t),d(t)) are smoother than the case for k = 1. Comparing Fig.5 and Fig.6 , it can be seen that the size and trend reflected by these two density functions in different forms are almost similar. In addition to the Weibull distributions (k = 1, 2), the distribution function

can also be considered, which is an modification of the Gaussian distribution. For given τ * , σ > 0, we can first calculate c * > 0 by the normalization of the density function, and then calculate its mean value τ > 0 for the specified density function f * (t). For density function in this form, we show the corresponding numerical performances for the second form of s(t) in (4.2) in Fig.7 , which are almost the same as those for the Weibull density functions. That is, our proposed model for COVID-19 is of some robustness for the concrete forms of density functions, and consequently may be applicable for different density functions. Next, we will check our proposed model using some officially issued COVID-19 data, by comparing some quantitative values. As mentioned earlier, the daily existing infected patients I s (t) is a crucial index of epidemic disease, so we should compare it produced by our model with the practical data of COVID-19 in principle. Unfortunately, there is no information about I s (t) officially issued. On the other hand, I s (t) is affected by various factors, especially the imported cases as shown in Fig.1-Fig.7 . Due to these two reasons, instead of I s (t), we consider the newly confirmed case j(t) from our model with the practical data. Since Wuhan, the earliest COVID-19 outbreak city in the world, goes into lockdown from January 23, 2020 and curbs population flow, we can take its external input s(t) = 0 to remove the uncertainty of the input cases. The other advantage for checking j(t) in Wuhan is that the newly confirmed data j(t) (also J(t) := t 0 j(τ )dτ ) is officially issued. Our model transforms the dynamical system of epidemic spread described by multi index (I s (t), j(t), c(t), d(t)) into a linear integro-differential equation with respect to I s (t). Once we solve I s (t), we can get j(t) as well as (c(t), d(t)). However, when solving I s (t), in addition to the model parameters in the density functions, we also need the initial value I s (0) = I 0 s , which is not available in the officially issued data and difficult to be fixed. To check our model avoiding this difficulty, we choose I 0 s as well as the model parameters by trial and error. The numerical simulations show that, for properly chosen I 0 s and the model parameters, the newly confirmed cases j(t) (especially the cumulative value J(t) on [0, t]) from our model is consistent with the number of confirmed cases officially issued. This fact ensures qualitatively that our model can describe the spread of the epidemic. Moreover, if we can identify the model parameters well from the partially existing epidemic data, which is a nonlinear inverse problem in terms of the governed system, our proposed model can give an effective prediction of the epidemic or provide some epidemic data hard to be measured directly, which is helpful for decision-making.

We take January 23, 2020 as the initial time with the assumed existing infections I 0 s = 43300. Since the continuous system is discretized by time step 1, we take 0 = 1.01 for considering the propagation of the initial infected persons. We check our model for the Weibull distribution with k = 2 and the generalized Gaussian distribution respectively, with the infection rate, isolation rate, cure rate of isolated patients and self-healing rate of the exposed patients taken as β = 0.467, κ(t) ≡ 0.989, κ 2 = 0.9, κ 4 = 0.6, while the mean value of each random index is taken as τ 1 = 15, τ 2 = 11, τ 3 = 15, τ 4 = 18, τ 5 = 10.

For the Weibull distribution with k = 2, we calculate j(t) from January 23 to March 18 by our model using these parameters, and compared them with the corresponding officially issued data. The results are shown in Fig.8(right) , and the comparisons of the cumulative confirmed infected cases J(t) are also shown in Fig.8(left) , while the corresponding comparisons using the generalized Gaussian distribution are shown in Fig.9 . By comparing Fig.8 and Fig.9 , we observe that our simulation data are not sensitive to the form of density functions. be found that our simulation data (red) are higher than the officially data (blue) between January 25 and February 13. However, we notice that the officially issued data on February 13 is of some sharp jump, indicating that the government has an overall compensations to the lower data released early. If we add these positive error into the issued data between January 25 and February 13, our simulations will match the actual data better.

We can further give the quantitative comparisons on (j(t), J(t)) with officially data in Table 1 . As already analyzed, there was a large error on the daily confirmed data at the initial stage. However, after the adjustment by the officially issued data on February 13, the numerical data match with the actual data well after 13 February. We also apply our proposed model to simulate the spread of COVID-19 in Italy, noticing that the main cities of Italy also went into lockdown. For the Rayleigh distribution, we take the parameters β = 0.475, κ(t) = 0.951, κ 2 = 0.9, κ 4 = 0.6 with initial data I s (0) = 110000 at 10 March. The average periods are assumed to be τ 1 = 18, τ 2 = 11, τ 3 = 15, τ 4 = 18, τ 5 = 10.

With these specified parameters, our model for (j(t), J(t)) can match the officially issued data very well from 10 March to 20 April as shown in Fig.10 . Just like we already explain, if we smoothen the random oscillation of daily confirmed case j(t) and consider its total quantity J(t) = t 0 j(s)ds, our predictions match the issued data well as shown in Fig.10  (left) . We expect that our predication data after 20 April can also match with the officially issued data of Italy up to some extent, which should be verified by the actual situations.

Since the existing infected patients consists of many ingredients such as the patients at different stages, patients either in isolation or in exposure, the quantity I s (t) is almost impossible to get from the statistic process. However, since the simulation results for (j(t), J(t)) by our model with specified model parameters and initial data match with the officially issued data in Italy, as shown in Fig.10 , it is reasonable to make predications on I s (t) from our model with the same model parameters. Using the same density function, our simulations on I s (t) in Italy from 10 March to 17 June are shown in Fig.11 , although these data are not officially available at present stage. It is not so optimistic from our simulations, since that COVID-19 in Italy will continuously spread after 17 June according to Fig.11 .

our model in generally can reflect the number of confirmed cases j(t) (right), especially the cumulative data J(t) = t 0 j(s)ds (left)

A time delay dynamic system with external source for the local outbreak of 2019-nCoV

Modeling and prediction for the trend of outbreak of NCP based on a time-delay dynamic system

When will be the resumption of work in Wuhan and its surrounding areas during COVID-19 epidemic? A data-driven network modeling analysis (in Chinese)

Modified SEIR and AI prediction of the epidemics trend of COVID-19 in China under public health interventions

Updated understanding of the outbreak of 2019 novel coronavirus (2019-nCoV) in Wuhan

Permanence of a discrete SIRS epidemic model with time delays

Global behavior and permanence of SIRS epidemic model with time delay

Global asymptotic stability of an SIR epidemic model with distribured time delay

Acknowledgements. This work is supported by NSFC (Nos. 11971104, 11531005).