Stability analysis of a fractional ordered COVID-19 model


 The main purpose of this work is to study transmission dynamics of COVID-19 in Italy 2020, where the first case of Coronavirus disease 2019 (COVID-19) in Italy was reported on 31st January 2020. Taking into account the uncertainty due to the limited information about the Coronavirus (COVID-19), we have taken the modified Susceptible-Asymptomatic-Infectious-Recovered (SAIR) compartmental model under fractional order framework. We have formulated our model by subdividing infectious compartment into two sub compartments (reported and unreported) and introduced hospitalized class. In this work, we have studied the local and global stability of the system at different equilibrium points (disease free and endemic) and calculated sensitivity index for Italy scenario. The validity of the model is justified by comparing real data with the results obtained from simulations.


Introduction
The respiratory disease caused by novel Coronavirus was rst observed in Wuhan City, Hubei Province, China in December 2019. The outbreak of the disease is ongoing worldwide and the World Health Organization named it Coronavirus disease 2019 (COVID-19) on 11 February 2020. On February 21, 2020, the rst Italian patient with COVID-19 was diagnosed, a 38-year-old man hospitalized at Codogno Hospital, Lodi, in northern Italy. As of 3 September 2020, Italy has 28,915 active cases; during the elevation of the pandemic, Italy's number of active cases was one of the highest in the World. By 3 September, Italy had tested nearly 5,342,000 people. Due to the limited number of tests performed, the actual number of infected people in Italy, as in other nations, is estimated to be higher than the o cial count [37]. Due to this reason we have considered two sub classes, namely 'reported symptomatically infected class' and 'unreported symptomatically infected class'. Maria Van Kerkhove, the WHO's technical lead on the COVID-19 pandemic, made it very clear that the actual rates of asymptomatic transmission are not however known [36]. So, we have considered the fact that asymptomatic class does not transmit disease and the main spreading of the Corona virus is happened due to reported infected class and unreported infected class. Some proportion (critical conditions ) of reported infected class should be transferred to hospitalized class. Hospitalization capacity is limited for every state, especially ICU facilities and due to this hospitalization of COVID-19 patients have created a capital issue. Introducing hospitalized class, we want to predict the number of COVID-19 patients would be hospitalized which was a major issue in the beginning of pandemic in Italy 2020.
Fractional calculus can be regarded as the generalization of their order where fractional order is replaced with integer order [24,23]. In systematic study it has been observed that integer order model is a limited case of fractional order model where the solution of the fractional order system must converge to the solution of integer order system as the order approaches to one [31]. In that respect there are many elds where fractional order systems are more suitable than integer order systems. Phenomena, which are linked with memory property and a ected by hereditary property, cannot be expressed by integer order systems [11]. It is observed that the data garnered from real life phenomena t better with fractional order systems. We have already done some works in fractional order dynamics [5,6,7,8,9]. Many researchers have contributed remarkably on various models of COVID-19 [14,17,20,26,27,28,29]. Fractional order modeling is a bene cial approach that has been used to study the nature of diseases because the fractional derivative is a generalization of the integer-order derivative. Also, the integer derivative is local in nature, while the fractional derivative is global. This behavior is very useful for modeling epidemic problems. In addition, the fractional derivative serves to enhance the stability region of the system. The fractional order system provides extra parameter which is useful for better numerical simulations. The previous models are very much productive for studying the transmission of COVID-19 but these models do not contain hospitalized, reported and unreported classes and also the pandemic situation of Italy is not mentioned in the previous research works. These facts along with the advantages of fractional calculus motivate us to construct our model on COVID-19 in Caputo fractional di erential equation systems.
In this work, we have constructed a modi ed fractional order SAIR model for two sub compartment infected classes (reported and unreported) and hospitalized class (section 2). Next, we have calculated basic reproduction number of the model which is very important to know whether a virus is highly contagious or not (section 3). Followed by the uniqueness, non-negativity, boundedness of the solutions have been shown (section 4). We have analyzed local and global stability of our proposed system and also calculated sensitivity indices (section 4). Further we have numerically analyzed the dynamical system with respect to the parameter values connected with Italy scenario (section 5). Section 6 contains some important conclusions.

Model formulation
We have constructed the following six compartmental model under fractional order framework using Caputo fractional di erential equations: where < α < , and C t D α t is the notation due to Caputo fractional derivative, t ≥ is the initial time (it is assumed that t = ). System (I) is dimentionally correct as both sides have same time dimension time −α . We have ignored the upper script α of all parameters and the model becomes: ( ) It is assumed that δ ≤ δu , δv , δ h i.e., the natural death rate (δ) is lower than the disease induced death rates (δu , δv , δ h ). Here S(t), A(t), U(t), V(t), H(t) and R(t) represent the susceptible, asymptomatically infected, reported symptomatically infected, unreported symptomatically infected, hospitalized population and recovered or removed population at time t respectively. The description of all parameters is given in Table 1.

Equilibria and Basic reproduction number
System (1) has two equilibrium points.  Here For E to exist in feasible region R + , it is necessary and su cient that The basic reproduction number R provides the average number of secondary (infectious) cases produced by one infective individual. It is calculated as the maximum eigenvalue of the next generation matrix FV − at disease free equilibrium E [32], where Here the rst part is due to reported symptomatically infected individuals and the second part is for unreported symptomatically infected individuals. Clearly, each parameter depends on α and so R is a function of α. For analysis purposes, we have xed the value of α. If we change the value of α, then all other parametric values will be changed and this will change the value of R . According to Figure 8, it is clear that the basic reproduction number R will be decreased if the value of α is diminished.

Analysis of the model
In this section we have veri ed the existence, uniqueness, non-negative, boundedness criterion of the solution of system (1). Next we have performed stability analysis.

. Preliminaries
Let us recall some basic theories that are needed for dynamical analysis.

De nition 1 [23]
The Caputo fractional derivative with order α > for a function g ∈ C n ([t , ∞+),IR) is denoted and de ned as: where Γ(·) is the Gamma function, t ≥ t and n is a natural number. In particular, for α ∈ ( , ): Remark:
where F(s) = L g(t) denotes the Laplace transform of g(t).

Theorem 2 [16]
Let C be the complex plane. For any α , α ∈ R+ and M ∈ C, then , represents the real part of the complex number s, and Eα ,α is the Mittag-Le er function.

Theorem 3 [23]
Consider the following fractional-order system: The equilibrium points of this system are evaluated by solving the following system of equations: Ψ(X) = . These equilibrium points are locally asymptotically stable i each eigenvalue λ i of the Jacobian matrix J( ., x n ) calculated at the respective equilibrium points satisfy arg(λ i ) > απ .

. Existence and uniqueness
Lemma 2 [19] Consider the system: x) with respect to x, then there exists a solution of (4) on [t , ∞) × Ω which is unique.
To study the existence and uniqueness of system (1), let us consider the region For any X, X ∈ Ω: Hence L(X) satis es Lipschitz's condition with respect to X. Therefore, Lemma 2 con rms that there exists a unique solution X(t) of system (1) with initial condition The following theorem is the consequence of this result.

. Non-negativity and boundedness
In this section we have established the criterion for feasibility of the solution of system (1). Suppose IR+ stands for the set of all non-negative real numbers and Γ+ = (S, A, U, V , H, R) ∈ IR + . Proof: From Lemma 1, we can say S(t) is increasing in a neighborhood of time t where S(t) = and S(t) cannot cross the axis S(t) = . Hence, S(t) > for all t ≥ . Now, we claim that the solution A(t) starts from Γ+ and remains non-negative. If not then there exists τ such that A(t) crosses A(t) = axis at t = τ for the rst time and the following conditions hold . Now, two cases arise. Case1: If ω U(τ )S(τ ) + ω V(τ )S(τ ) > then by the remark of Lemma 1 we can say that A(t) is nondecreasing in a neighborhood of t = τ and which concludes A(τ + ) = . Hence, we have arrived at a contradiction. Case2: In this case we have ω U(τ )S(τ ) + ω V(τ )S(τ ) < , which implies any one or both of U(τ ), V(τ ) must be negative. Now we have two sub-cases. Sub-case1: If V(τ ) < , then there exists a τ such that < τ < τ with which contradicts our assumption that V(τ + ) < . Therefore, we have V(t) ≥ , ∀t ∈ [ , ∞).
Thus, all solutions of system (1) starting in Γ+ are con ned in the region Γ+.
Proof: From rst equation of (1), it is noted that Taking Laplace transforms on both sides, we have Taking inverse Laplace transforms (using Theorem 2): where M = max Λ δ , S( ) and since from the properties of Mittag Le er function [15]: represents the total population, then Applying Laplace transformation, we have (using Theorem 1): Taking inverse Laplace transforms (using Theorem 2): (6) and (7), we get

. Local stability
For stability analysis, let us reduce system (1) by discarding last two equations as H, R do not appear in rst four equations of system (1). If we study the dynamics of S, A, U, V, then the dynamics of H, R also be obtained from them. The reduced system is as follows: ( ) The equilibrium points of (8) are as follows: 1. Disease free equilibrium: E * = ( Λ δ , , , ) 2. Endemic equilibrium: E * = (S * , A * , U * , V * ). Here For E * to exist in feasible region R + , it is necessary and su cient that To analyze the local stability of endemic equilibrium E * , we need the followings: For n = , we have f (x) = x + α x + α x + α and f (x) = x + α x + α .
For n = , P(λ) = λ + a λ + a λ + a λ + a a  a  a  a  a  a  a  a  a  a  a  a  a  a  a  a  a  a  a a a
4. If ∇(P) < , a > , a > , a > , a > and a = a a a + a a , then the equilibrium point is asymptotically stable for α ∈ ( , ). 5. a > is the necessary condition for the equilibrium point E * is to be locally asymptotically stable.
At the disease-free equilibrium point E * , the Jacobian matrix is given by The following theorem is the consequence of these discussions.

Theorem 7
The disease free equilibrium E * of system (1) is asymptotically stable if any of the following conditions hold: where K i , i = , ..., are given in (10). Now, Jacobian matrix of system (1) at endemic equilibrium E * is given by Therefore, we have the following theorem: The endemic equilibrium E * of system (1) is asymptotically stable if any of the following conditions hold: where Q(λ) = λ + d λ + d λ + d λ + d is the Jacobian matrix at endemic equilibrium E * and d i 's ( i=1,2,3,4) can be found from (11) and (12).

. Global Stability
First, let us state the following important lemmas.

Lemma 6 [10](Uniform Asymptotic Stability Theorem) Consider the non-autonomous system
Let x * be an equilibrium point of the system (x * ∈ Ω ⊆ R n ) and Φ(t, x(t)) : [ , ∞) × Ω → R be a continuously di erentiable function such that Clearly L ≥ and L = only at ( , , ). Taking α order Caputo derivative C D α t of L along the solution of system (8), we have and ω M ≤ (η + ξ + δv). Therefore, using Lemma 6: lim Hence in the limit S(t) is given by the solutions of C D α t S(t) = Λ − δS. Since S( ) > , the theorem follows.
Proof: Consider a positive de nite function: It is observed that V ≥ and V = only at E * . Taking α order Caputo derivative C D α t of V and using Lemma 5, we have From steady-state of equilibrium point (9), we have (15) and (16), we have From relation (17) it is clear that C D α t (V) ≤ and thus C D α t (V) is negative de nite with respect to E * . Thus E * is globally asymptotically stable by Lemma 6.
. α α . Table 3: Sensitivity indices of di erent parameters of system (1) corresponding to Table 2 Parameters sensitivity index

. Sensitivity analysis
The basic reproduction number (R ) depends on several parameters and value of R = . according to Table 2. To examine the sensitivity of R to any parameter (say, θ), normalized forward sensitivity index with respect to each parameter has been computed as follows [3,32]: The sensitivity index may depend on some system parameters but also can be constant or independent of some parameters. These values are very much important to estimate the sensitivity of parameters which should be done cautiously, since a small perturbation in a parameter causes relevant quantitative changes. Merely in the estimation of a parameter with lower value of sensitivity index does not demand to deal cautiously, because a small perturbation in that parameter causes small changes. Since we cannot control δu , δ, δv, we have concentrated on the parameters ω , ω , σ, ϕ, ξ , ϵ, η, f . Now, Ω R θ = + indicates that increasing (decreasing) θ by a given proportion increases (decreases) R by same proportion and Ω R θ = − means that increasing (decreasing) θ by a given proportion decreases (increases) R by same proportion. The values of sensitivity indexes for the parameters ω , ω , σ, ϕ, ξ , ϵ, η, f corresponding to Table 2 is given in  Table 3.

Numerical simulations
For numerical simulations, we have used MatLab interface along with Predictor-corrector PECE method for fractional di erential equations introduced by Roberto Garrappa [12]. We have considered Table 2 for the scenario of Italy due to COVID-19 in the period 1st March 2020 to 20th April 2020. We have performed numerical simulations to compare the results of our model with the real data from various reports published by WHO [33] and worldometer [4]. The total population of Italy is around . × [35]. We have taken t = day as time unit and t = as the nal time. The recruitment rate has been calculated as . ×N(= . × ) × = per day. In Italy, birth rate per 1000 inhabitants is 7.2 [30]. Table 4 recommends the reported symptomatically infected cases and deaths from 1st March 2020 to 20th April 2020, we have taken 5 days gap between two successive reports [4]. We have considered system (I) rather than system (1) for numerical simulations as the parameters of system (I) contains α which gives better results. From Figure 2, it has been observed that in initial stage the time series of reported class is tted with actual data and for death rate the time series is relevant for rst 15 days ( Figure 3). Next, we have simulated our model with estimated parametric values ( Table  2) from 3rd September (present scenario in Italy) to 31st October (60 days) ( Figure 5-6). Figure 4 indicates that for α = . the time series of reported class is closer to the realistic scenario. The death cases have been computed at time t due to COVID-19 as δu U(t) + δv V(t) + δ h H(t). The number of death cases according to the time series on / / , / / , / / , / / are 29, 168, 327, 435 which are quite close to the real data. Similarly, on / / , / / , / / , / / number of reported cases according to time series are 36370, 59150, 78430, 93500. Figure 7 depicts that for slight changes in η the reported cases increases, but the unreported cases decreases rapidly and number of hospitalized class is also increased. The variation of R with order of derivative is depicted in Figure 8.

Conclusion
Fractional calculus plays an important role in real dynamical processes, including the cases of epidemic spreading. Here we have studied on the evolution of a modi ed SAIR epidemic model, incorporating memory e ects. We have taken the hospitalized class and subdivided symptomatically infected class into two subclasses, one is reported class and other is unreported class. Our model is giving a well approximation of the reality of Italy outbreak and predicting daily number of con rmed cases due to COVID-19. The number of hospitalized persons is relevant to provide an approximation of the Intensive Care Units (ICU) required. The model we have studied for Italy can also be employed to study the reality of the other countries. The value of order of derivative α may vary region to region. It is observed that the dynamics of the system (1) depends on  Table 2 and real data (from 1st march 2020 to 20th April 2020)  Table 2 di erent α = . , . , .   Table 2 the robustness of memory e ects, which is controlled by the order of fractional derivative α. The results will be di erent if we change the order of derivative for same set of parametric values (Figure 4). This shows that order of derivative plays an important role in simulations of system (1). From Table 2, it is clear that the values of parameter depends on order of derivative α. We will get di erent sets of parametric values if we change the value of α. In our context, we have xed the value of di erentiation at 0.85 which is more suitable for the real scenario. From Figures 5-6, we can conclude that the number of reported case will be decreasing at the end of October 2020. The morbidity rate will also be under control in the last week of October 2020. Rapid PCR (Polymerase Chain Reaction) will increase the reported cases which may convey the situation under control. No model is perfect for COVID-19 scenario. There are many factors required to be considered for modeling COVID-19. The situation in di erent countries in this pandemic is unlike. From analytical study of our model it seems to us that if we are able to decrease disease transmission rates (ω , ω ), then we will achieve disease-free state quickly. These are also consistent with the conditions of stability (see Theorem 9). The social distancing, mobility of asymptomatic individual, herd immunity, vaccination, improvement of health conditions are the key factors and essential in complex and accurate modeling. Though our model has these limitations, but even then this model can ray on the research on Mathematical study of COVID-19 spreading.