key: cord-0600332-rwuy73h4
authors: Lu, Zhenzhen; Yu, Yongguang; Chen, YangQuan; Ren, Guojian; Xu, Conghui; Wang, Shuhui; Yin, Zhe
title: A fractional-order SEIHDR model for COVID-19 with inter-city networked coupling effects
date: 2020-04-26
journal: nan
DOI: nan
sha: bd338ce20c06653892823f53330d4497af7ffb01
doc_id: 600332
cord_uid: rwuy73h4

In this paper, a mathematical model is proposed to analyze the dynamic behavior of COVID-19. Based on inter-city networked coupling effects, a fractional-order SEIHDR system with the real-data from 23 January to 18 March, 2020 of COVID-19 is discussed. Meanwhile, hospitalized individuals and the mortality rates of three types of individuals (exposed, infected and hospitalized) are firstly taken into account in the proposed model. And infectivity of individuals during incubation is also considered in this paper. By applying least squares method and predictor-correctors scheme, the numerical solutions of the proposed system in the absence of the inter-city network and with the inter-city network are stimulated by using the real-data from 23 January to $18-m$ March, 2020 where $m$ is equal to the number of prediction days. Compared with integer-order system ($alpha=0$), the fractional-order model without network is validated to have a better fitting of the data on Beijing, Shanghai, Wuhan, Huanggang and other cities. In contrast to the case without network, the results indicate that the inter-city network system may be not a significant case to virus spreading for China because of the lock down and quarantine measures, however, it may have an impact on cities that have not adopted city closure. Meanwhile, the proposed model better fits the data from 24 February to 31, March in Italy, and the peak number of confirmed people is also predicted by this fraction-order model. Furthermore, the existence and uniqueness of a bounded solution under the initial condition are considered in the proposed system. Afterwards, the basic reproduction number $R_0$ is analyzed and it is found to hold a threshold: the disease-free equilibrium point is locally asymptotically stable when $R_0le 1$, which provides a theoretical basis for whether COVID-19 will become a pandemic in the future.

In December 2019, the first case of respiratory disease caused by a novel coronavirus was identified in Wuhan City, Hubei Province, China. Contrary to the initial observations [8] , COVID-19 does spread from person to person as confirmed in [6] . It quickly spread to all parts of the country and parts of Southeast Asia, North America and Europe. As of 5 February, 2020, more than 24550 cases of coronavirus disease 2019 (COVID-19) had been confirmed, including over 190 cases outside of China, and more than 490 reported deaths globally. Several intervention strategies have been implemented in China in order to contain the epidemic. Whereas assessing the intervention measures of COVID-19 epidemic poses a major health concern.

Wuhan is the capital of Hubei Province, as the major air and train transportation hub of central China with huge population movement during the period of spring festival. The geographical factors plus self-sustaining human-to-human transmission in the community of the new coronavirus have added radical difficulties to the prevention and interference of this epidemic. The outflow of domestic passenger volumes from Wuhan was estimated to be 5 million during the spring festival, about 0.7 of which went to other cities in Hubei province. which made it difficult for the infected people to track back to the population. Government has progressively implemented metropolitan-wide quarantine of Wuhan since 23-24 January, 2020 which ultimately contributes a lot of the containment of virus. The epidemic was seeded by a force of frequent human contact generating from communication. Prasse et al. proposed a network-based discrete model to describe this condition and they found the networkbased modelling was beneficial for an accurate forecast of the epidemic outbreak [22] . Peng et al. constructed and analyzed a generalized SEIR model to analyze COVID- 19 and their work showed that the outbreak of COVID-19 in China could be dated back to the end of December 2019 [21] . In addition, the basic reproduction number R 0 is 'the expected number of secondary cases produced, in a completely susceptible population, by a typical infective individual'. If R 0 ≤ 1, on average an infected individual produces less than one new infected individual over the course of its infectious period, and the infection cannot grow. Conversely, if R 0 > 1, each infected individual produces, on average, more than one new infection, and the disease can invade the population [10]. There have been abundant articles estimating the basic reproduction number R 0 of COVID-19 to determine whether the epidemic disease will spread widely [15, 18, 24] . Meanwhile, in the process of infectious disease transmission, susceptible individuals contact infected individuals and become infected with a certain probability. There is a lot of evidence that the incidence rate is an important tool to describe this process [4, 5, 28] . It represents the infection ability of a infected individuals in per-unit time. So because of high infectiousness of COVID-19, the bilinear incidence rate is always considered. And the classical SEIR model assumes that the incubation of infected person is not infectious, but this assumption is quite different from the infection characteristics of the new coronavirus infection COVID-19 [27] . What's more, due to a variety of special conditions, patients with COVID-19 cannot be admitted to hospital for immediate treatment, and individuals infected with COVID-19 have a high mortality rate both in the incubation period and notreatment period.

It is worthing to point out that the state of epidemic models does not depend on its history in the classical integerorder epidemic model. However, in real life, the spread of infectious diseases depends not only on its current state, but also on its past state. Smethurst et al. find that the waiting times between doctor visits for a patient follow a power law model [26] . At the same time, a power law waiting time distribution P[J n > t] = Bx −α leads to Caputo fractional time derivative C t 0 D α t of the same order [19] . And the Caputo fractionalorder derivatives allows traditional initial and boundary conditions when dealing with real-world problems. Not only that, its derivative for a constant is zero. What's more, the precision of the Caputo fractional-order can supersede the integer order resulting from its changes at every instant of time and nonlocal behavior. Khan et al. considered a fractionalorder model to describe the brief details of interactional among the bats and unknown hosts, then among the people and the infections reservoir (seafood market) and found the fractional model can be helpful for the infection minimization [14] . Yu et al. establish a novel fractional time delay dynamic system (FTDD) to describe the local outbreak of COVID-19 and the reconstructed coefficients are used to predict the trend of the Corona-Virus [7] . Shaikh et al. estimated the effectiveness of preventive measures and various mitigations, predicting future outbreaks and potential control strategies using a Bats-Hosts-Reservoir-People transmission fractionalorder COVID-19 model [3] . Xu et al. proposed a generalized fractional-order SEIQRD model and this model has a basis guiding significance for the prediction of the possible outbreak of infectious diseases like COVID-19 and other insect diseases in the future [29] .

The purpose of this paper is to incorporate time-fractionalorder and inter-city networked coupling effects into a SEI-HDR epidemic model to investigate the dynamic behavior of COVID-19. There are a lot of research for COVID-19 that not only individuals are contagious, but those in the incubation period may have the same incidence rate as well as those with symptoms. So the infectiousness of the infected person's incubation period is considered in this paper. Also, considering an individual hospitalized with COVID-19, the different mortality rates are investigated in this model. Thus based on the above analysis and the memorability of Caputo fractional-order derivatives, a fractional-order SEIHDR model for COVID-19 with inter-city networked coupling effects are built. Should such a model can be fitted with data reasonably well, some epidemic parameters can be extracted. Then, based on the official data given by NHC (National Health Commission of the people's Republic of China [1]) every day, several numerical examples are exhibited to verify the rationality of the fractional-order model and effectiveness of inter-city networked coupling. Compared with integer-order system (α = 0), the fractional-order model without network is validated to have a better fitting of the data on Beijing, Shanghai, Wuhan, Huanggang and other cities. In the case of inter-city networked coupling effect, the results indicate that this system may be not a significant case to virus spreading for China, but may have an impact on preventive measures in other countries. Finally, the stability of disease-free equilibrium point is obtained by the basic reproduction number R 0 , which provide theoretical significance for the development of COVID-19 in futher.

The rest of the paper is organized as follows. In Section 2, a fractional SEIHDR epidemic model for COVID-19 is formulated. In Section 3, some dynamical behaviors of the proposed system are analyzed. In Section 4, some numerical simulations are presented to illustrate theoretical results according to the official dates. Finally, a brief discussion is given in Section 5.

Fractional-order calculus have been found wide applications to model dynamics processes in many well-known fields of science, engineering, biology, medicine, and many other [11, 12, 20, 23, 25] . Before presenting the epidemic model in fractional derivative, some necessary preliminaries are introduced.

In this subsection, some definitions and results are introduced firstly.

Definition 2.1. [30] A Gamma function of α > 0 is defined by

For any t > t 0 , the time Caputo fractional derivative of order α (n < α < n + 1) with the lower

A constant x * is an equilibrium point of the Caputo fractional dynamical system:

Lemma 2.1. [17] Consider the fractional-order system:

x) satisfies the local Lipschitz condition with respect to x, there exists a unique solution of the above system. Lemma 2.2. [13] Consider the following fractional-order system:

The equilibrium points of the above system are calculated by solving the following equation: f (x) = 0. These points are locally asymptotically stable if and only if all eigenvalues λ of the Jacobian matrix J = ∂ f ∂ x at x evaluated of the equilibrium points satisfy |arg(λ )| > απ 2 .

A novel coronavirus disease, named COVID-19, broke out in Wuhan, Hubei province, China, in December 2019. A great deal of epidemic models exists to describe spread of infectious diseases mathematically [2] . However, the outbreak of COVID-19 began during the period of spring festival in China, which large population movement made cities more interconnected. Meanwhile, Wuhan is the capital of Hubei Province, as the major air and train transportation hub of central China, which makes it necessary to built a inter-city networked to analyze the spread of COVID-19. And due to limited medical conditions, individuals in incubation and infection periods cannot be hospitalized immediately for treatment, which increases the mortality of COVID-19. Moreover, Tang et al proposed that the incubation of infected person of COVID-19 is infectious [27] . Based on the above analysis, a fractional-order SEIHDR model with inter-city networked coupling effects in this paper is developed as follows:

with the initial condition

The total population of city k is denoted by N k which is further classified into six subgroups where S k (t), E k (t), I k (t), H k (t), D k (t) and R k (t) (k = 1, 2, ..., n) present the number of the susceptible, exposed, infective (infected but not hospitalized), hospitalization, death and recovered individuals at time t and city k, respectively. The susceptible people S k infected after the interaction with E j and I j , given by

and hospitalization κ k (t)H k , where µ ik (i = 1, 2) and κ k (t) implies the disease-related death rate. The parameters λ k (t) be the recovery rate through hospitalization; r k denotes the transit rate of the exposed individuals E k ; δ k be hospitalization rate of the infective individuals; t = 0 represents 23 January, 2020. Furthermore, β k j , µ ik , δ k and r k (i = 1, 2) are positive constants; λ k (t) and κ k (t) are bounded function, i.e. |λ k (t)| ≤ M 1k and |κ k (t)| ≤ M 2k for all t ≥ 0, where M 1k and M 2k are positive constants.

In this section, some dynamical behaviors of the proposed system (1) are analyzed.

It is significant to demonstrate the existence, uniqueness and boundedness of a nonnegative solutions for system (1) before implementing its numerical process. Thus this subsection moves to the discussion of proprieties mentioned above.

Theorem 3.1. For any given initial condition (S k (0), E k (0), 0, 0, 0, 0, 0, 0) , there exists a unique nonnegative and boundedenss solution (S k (t), E k (t), I k (t), H k (t), R k (t), D k (t)) of system (1) for all k = 1, 2, ..., n and t > 0 where t = 0 represents 23 January, 2020.

Proof:

Obviously, one has

where L k = max(L 1k , L 2k , L 3k , L 4k , L 5k ), L 1k = 2β k j M N , L 2k = L 1k + µ 1k + r k , L 3k = r k + µ 2k + δ k , L 4k = M 1k + M 2k and L 5k = µ 1k + µ 2k + M 2k . So F k satisfies the local Lipschitz condition with respect to X k . Then there exists a unique boundedenss solution (S k (t), E k (t), I k (t), H k (t), R k (t), D k (t)) 1≤k≤n of system (1). Furthermore, consider the following auxiliary system:

Through the comparison theorem, it is not difficult to find that (S k , E k , I k , H k , R k , D k ) = (0, 0, 0, 0, 0, 0) is the lower solution of system (1). Thus, one has S k ≥ 0, E k ≥ 0, I k ≥ 0, H k ≥ 0, R k ≥ 0 and D k ≥ 0 for t ≥ 0 and k = 1, 2, ..., n.

Turning now to the exploration on local stability of system (1) by considering first the disease-free equilibrium and the basic reproduction number denoted by R 0 . According to system (1), the disease-free equilibrium for system (1) is E 0 = (S 1 (0), 0, 0, 0, 0, 0, ..., S n (0), 0, 0, 0, 0, 0) where S k (0) (k = 1, 2, ..., n) are the initial condition of system (1).

Theorem 3.2. The basic reproduction number of system (1) is

), V 10 = diag(µ 11 +r 1 , ..., µ 1n +r n ), V 20 = diag(−r 1 , ..., −r n ), V 21 = diag(δ 1 + µ 21 , ..., δ n + µ 2n ) and

V 02 = (δ k I k + µ 2k I k −r k E k ), V 03 = (λ k H k +κ k H k −δ k I k ). Then taking the derivative of F 0 , V 01 , V 02 and V 03 with respect to E k , I k and H k (k = 1, 2, ..., n) at the disease-free equilibrium point E 0 , one has

, j≤n , V 10 = diag(µ 11 +r 1 , ...., µ 1n +r n ), V 20 = diag(−r 1 , ...., −r n ), V 21 = diag(δ 1 + µ 21 , ...., δ n + µ 2n ), V 30 = diag(−δ 1 , ...., −δ n ) and V 31 = diag(λ 1 + κ 1 , ...., λ n + κ n ). Then according to the definition of the basic reproduction number [10], one has

Based on the basic reproduction number R 0 , the following theorem can be taken into consideration: Theorem 3.3. If R 0 ≤ 1, the disease-free equilibrium point E 0 = (S 1 (0), 0, 0, 0, 0, 0, ..., S n (0), 0, 0, 0, 0, 0) of system (1) is locally asymptotic stability.

Proof: The disease-free equilibrium point E 0 of system (1) is locally asymptotically stable if all eigenvalues of the Jacobian matrix of system (1) at E 0 namely,

satisfies |arg(s)| > απ 2 and unstable if for some eigenvalues s, |arg(s)| < απ 2 , where V 41 = diag(λ 1 (t), ...., λ n (t)), V 51 = diag(µ 11 , ...., µ 1n ), V 61 = diag(µ 21 , ...., µ 2n ) and V 71 = diag (κ 1 (t), ...., κ n (t)). One can calculate that the eigenvalues are s 1 = 0 and s 2 = −s(V 31 ) and

Thus the disease-free equilibrium point E 0 = (S 1 (0), 0, 0, 0, 0, 0, ..., S n (0), 0, 0, 0, 0, 0) of system (1) is locally asymptotic stability.

Remark 3.1. When there is no interaction between city k and city j (i.e. city k is isolated from other cities), the basic reproduction number in city k is given:

In this case, the disease-free equilibrium point E 0 k = (S k (0), 0, 0, 0, 0, 0) of city k is local asymptotic stability when the basic number R k 0 ≤ 1. Remark 3.2. It is obvious that the basic number R k 0 is dependent from the onset λ k and κ k . The sensitivity of R k 0 to the other parameters β kk , µ 1k , r k , µ 2k and δ k are calculated as follows:

,

where A β kk , A µ 1k , A µ 2k and A δ k denote the normalized sensitivity indexes with respect to β kk , µ 1k , µ 2k and δ k , respectively. It is remarkable that the m times' increase in β kk results in the m times' increase in R k 0 , but the m times' increase in µ 1k , µ 2k , δ k and r k results in the m times' decrease in R k 0 .

From the previous discussion, it can be seen that system (1) is locally stable at the disease-free equilibrium point E 0 , and it has one unique bounded solution, which may provide a theoretical basis for the coronavirus control. This section presents the numerical stimulation of the coronavirus system (1) by fitting the real-data dating from 21 January to 18 March, 2020, given by NHC every day. System (1) is solved by applying least squares method [9] and predictorcorrectors scheme. To evaluate the prediction accuracy, the data can be moved for a fixed number of days, say m, prior to March 18. The prediction model is determined upon the data from 23 January up to 18 − m March, 2020. Then we predict the course of the disease up to 18 March, and the number of omitted days m is equal to the number of prediction days. It is reasonable to assume like in [21] from COVID-19 that the cure rates λ (t) and the disease-related death rate κ(t) vary over time as

where λ 0 and κ 0 are initial cure and death rate. In the absence of the network (i.e. k = 1), one can see from Fig. 1,  Fig. 2, Fig. 3, Fig. 4 and Fig. 5 that the fitting effect of the fraction-order system (1) is better than that of the integerorder system (i.e. α = 1) both the peak number and time of confirmed individuals. Particularly, the fractional-orders are α = 1.1913, α = 1.2873, α = 1.4343, α = 1.2403 and α = 0.0332, respectively. Corresponding the basic reproduction number R 0 k are 0.8855, 0.8833, 0.9848, 1.1491, 0.9568, respectively. Then the sensitivity of the basic reproduction number is presented from Fig. 6 , which is consistent with Remark 3.2. In addition Fig. 7 give a observation iconically on that as the disease transmission coefficient β kk and the hospitalization rate of infective individuals δ k change, the peak time and final infection size are given (if control measures before January 22 are to continue). The numerical results reveal that strict control of human contact, increased accounting reagents and clinical testing plays a critical role in reducing the number of confirmed case and postponing the occurrence of peak in Hubei province. However, it also can be found from Fig. 7 that when the disease transmission coefficient β kk increases, the number of the infected reaches its peak time earlier, meanwhile, its peak size increases. It can be induced that the measures taken by some countries are feasible and effective.

In the case of the inter-city network effect involving the cities in Hubei and their interactions (eg. traffic flow), the precise interactions between cities are known and must be inferred from system (1) by applying least squares method. Then the relationship between Beijing, Shanghai, Wuhan and Huanggang are described in Fig. 8 and Fig. 9 . It can be seen that the network-based model is no longer suitable for China, since China imposed city closures on 23 January. Nevertheless, with regard to population mobility in other countries, it can be hypothesized that the networkbased model shows a better performance in describing the disease situation which will be the ensuring research focus.

Finally, based on the data of symptoms, hospitalizations, deaths and recoveries in Italy from 24 February to 31 March, the short-term prediction for future outbreaks is realized in Fig. 10 . Then by calculation, the basic reproduction number R 0 k = 1.3368 × 10 3 is significantly greater than 1, which suggests that COVID-19 of Italy will not be eliminated in a short term. In addition, it can be seen from Fig.10 that the number of symptomatic individuals peaked on 30 June and the peak day of the number of hospitalizations occurred on 20 April, without surprise, the number of death and recovery individuals have increased.

In this paper, incorporating inter-city networked coupling effects for Beijing, Shanghai, Wuhan and Huanggang, a fractional order SEIHDR epidemic model is proposed. By applying least squares method and predictor-correctors scheme, the numerical solution of system (1) is compared with the real-date from 23 January to 18 March, 2020 about Beijing, Shanghai, Hubei, Wuhan and Huanggang. One of the most significant findings to emerge from this investigation is that the fractional-order system has a better fitting effect than the integer-order system. In addition, in the case of the inter-city network effect, the results indicate that the network system may be not a significant case to virus spreading for China. Meanwhile, system (1) better fits the data from 24 February to 31, March in Italy and the peak number of confirmed individuals in Italy is predicted according to system (1). Moreover, the existence, uniqueness, and positivity solution by the initial-value problem of system (1) are establish. Then the local stability of disease-free equilibrium point are studied by the basic reproduction number R 0 . And in the absence of network, the sensitivity of R k 0 to the other parameters is analyzed, which provide theoretical basis for disease control.

While this study did not confirm the epidemic model based on inter-city networked effects in China cases due to their strict lock down measures, it did partially substantiate that the fractional-order model based on network may yield more accurate predictions than modeling the epidemic for other countries independently. In spite of its efficiency in predicting the case of COVID-19 in China, this study has inspired several questions in need of further investigation: if the network model can be adapted for predicting the epidemic in other regions; how the factors like medical treatment affect the tendency of virus spreading and so on. It is recommended that further research be undertaken in predicting outbreaks in other countries with networked effects included in the dynamic model. 

A deterministic epidemic model for the emergence of COVID-19 in China

A Mathematical Model of COVID-19 Using Fractional Derivative: Outbreak in India with Dynamics of Transmission and Control

Infectious diseases of humans: dynamics and control

The mathematical theory of infectious diseases and its applications

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

The reconstruction and prediction algorithm of the fractional TDD for the local outbreak of COVID-19

-novel Coronavirus: Where We are and What We Know

Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission

Stability analysis of a fractional order model for the HIV/AIDS epidemic in a patchy environment

The effect of vaccines on backward bifurcation in a fractional order HIV model

Fractional-order Nonliear Systems: Modeling Analysis and Simulation

Modeling the dynamics of novel coronavirus (2019-nCov) with fractional derivative

Severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) and coronavirus disease-2019 (COVID-19): The epidemic and the challenges

Mittag-Leffler stability of fractional order nonlinear dynamic systems

Stability of fractional-order nonlinear dynamic systems: Lyapunov direct method and generalized Mittag-Leffler stability

The reproductive number of COVID-19 is higher compared to SARS coronavirus

Stochastic models for fractional calculus

Mathematical modelling and analysis of love dynamics: A fractional approach

Epidemic analysis of COVID-19 in China by dynamical modeling

Network-based prediction of the 2019-ncov epidemic outbreak in the Chinese province Hubei

Analysis of a fractional SEIR model with treatment

Estimating the reproduction number of COVID-19 in Iran using epidemic modeling

Diffusion process modeling by using fractional-order models

Are hospital waiting lists selfregulating?

Prediction of New Coronavirus Infection Based on a Modified SEIR Model. Cold Spring Harbor Laboratory

Dynamics of an SEIR epidemic model with nonlinear incidence and treatment rates

Forecast analysis of the epidemics trend of COVID-19 in the United States by a generalized fractional-order seir model. Nonlinear Dynam (Special Issue

Fractional differential equations

We declare that we have no conflict of interest.