Study of fractional order dynamics of nonlinear mathematical model

This manuscript is devoted to establishing some theoretical and numerical results for a nonlinear dynamical system under Caputo fractional order derivative. Further, the said system addresses an infectious disease like COVID-19. The proposed system involves natural death rates of susceptible, infected and recovered classes respectively. By using nonlinear analysis feasible region and boundedness have been established first in this study. Global and Local stability analysis along with basic reproduction number have also addressed by using the next generation matrix method. Upon using the fixed point approach, existence and uniqueness of the approximate solution for the mentioned problem has also investigated. Some stability results of Hyers-Ulam (H-U) type have also discussed. Further for numerical treatment, we have exercised two numerical schemes including modified Euler method (MEM) and nonstandard finite difference (NSFD) method. Further the two numerical schemes have also compared with respect to CPU time. Graphical presentations have been displayed corresponding to different fractional order by using some real data.


Introduction
Fractional calculus has got much attention in recent times. The foundation of the said area had been provided by Newton and Lebnitz during seventeenth century. Later on Reimann, Liouville, Fourier, Abel and Euler developed this branch very well [1]. The first notable definition had been given by Reimann-Liouville in 1832. Also the mentioned definition has been mod-ified by Caputo in 1967 to elaborate the concept more clearly in geometry [2,3]. Since fractional derivative infact is a definite integrals which can be defined in numbers of ways. Geometrically the said operator gives accumulation of a function which includes its integer counter part as a special case. Therefore, various researchers have given different definitions for fractional differential operators [4]. The said operators have numerous application in modeling of various real world problems. As in most of the situations where memory terms or index laws are involved cannot be well explained via classical operators. For comprehensive descriptions and well explanations fractional calculus is powerful tool to use (see [5,6]). Researchers have used the concept of fractional calculus for detailed explanations of many real world problems including physical, biological, dynamical and chemical phenomenons (we refer few article as [7][8][9][10][11]).
Recently various authors have used other kids of fractional differential operators like Hilfer, nonsingular kernel type, Caputo-Fabrizio operator to investigate various models. Here we state that each operator has own merits and de-merits. Here we remark that using various fractional differential operators numerous results relating controllability, existence theory, optimization and numerical analysis have been published in last few years. The concept of Hilfer derivative of fractional order has been applied to investigate some controllability results for delay type problems in [12,13]. In same line authors [14] have established existence and controllability results for nonlocal mixed Volterra-Fredholm type fractional delay integro-differential equations. Some advanced results in controllability using fractional order derivatives have been established in [15,16]. In same line authors have used neutral fractional order derivative to investigate a class of partial differential equations in [17]. Also authors have developed very applicable results for controllability and existence theory for neutral fractional integro-differential equations. For detail work, we refer [18][19][20]. Various results regarding analysis of various problems by using generalized proportional fractional integral operators, Atangana-Baleanue-Caputo operators have been established in last few years. Some recent contribution can be read as [21][22][23][24][25]. The analysis in aforesaid work has also applicable and can be further extended to COVID-19 or other infectious disease model. Here we remark that nonsingular kernel type and Caputo-Fabrizio operators suffer from initialization effect and extra assumptions need to be imposed on right hand side of a differential equations involving such like operators. This is not good effect and therefore now reducing the use of the said operators in applications. Although significant work has been done by using the aforesaid operators in epidemiology. But to the best of our knowledge in epidemiological models, the most powerful operators are those introduced by Reimann-Liouville and Caputo. Because they have clear geometrical interpretation like integer order differential equations and upon using integer value, the concerned equations reduce to their corresponding classical form.
The research area devoted to mathematical models for analysis of various physical phenomenons in real life is very interesting and has received wide attention among researchers in last several years. The idea of mathematical model had been originated by Bernoulli in 1776 (see [26]). Later on the concept adopted by Mekendric and Kermark in 1927 to describe an infectious disease model which is known as SIR. The said model then further modified and extended to form new models for various infectious disease like TB, Cholerea, Typhoid, Dengue, HIV and AIDs, etc (see [27,28]). Recently the Coronavirous disease (COVID- 19) which has been originated from China and transmitted all over the world. The aforesaid disease which reported in an outbreak in 2019 in Wuhan, Hubei province, China, is caused by the SARS-CoV-2 virus (for further etiological information see [29]). According to researchers of virology, the said virus belongs to the class of beta-Coronavirus and like Middle East Respiratory Syndrome Coronavirus and Severe Acute Respiratory Syndrome Coronavirus. The novel virus began to cause pneumonia, and was named as Coronavirus disease (2019) on 11 February 2020 by WHO (see detail in [30]). Also many researchers have formulated various mathematical models for COVID-19. This disease has greatly affected the whole globe. Due to COVID-19, more than fifty millions people have been died. Nearly five hundred million people have got infectious throughout the world (see some detail in [31,32]). The said disease has greatly destroyed the life style and economical situation of many countries around the globe. Recently some countries have worked on to prepare vaccine for permanent cure of the disease in which they have got some success like UK, USA, China and Germany, etc (see [33,34]). Researchers from bioengineering side, virology, epidemiology and those working on mathematical biology have worked significantly in last two years to investigate various procedure how to eradicate or control the disease from further spreading in community. Those who are working on computational biology have developed several mathematical models to understand the transmission dynamics of the said disease. Therefore large numbers of mathematical models have been formulated for COVID-19. Authors [42] have formulated the following SIRD type model to demonstrate the diffusion of the infection in various communities where S stand for susceptible, I for infected, R for recovered and D for death class. Further, d denote death rate, c is used for infection rates and a represents rate of conversion to susceptible. Recently various analysis, investigations and procedure have been conducted to deal the disease from various aspects. The contribution done in this regards, we refer few article as [35][36][37][38][39][40][41]. Since the concept of fractional calculus has got significant attention during last few decades. Therefore mathematician as well as researchers in biomathematics have given much attention to study said models for the mentioned infection under various type derivatives of arbitrary orders (see [43][44][45][46][47][48][49]). The investigation of biological model of infection disease under fractional derivatives is an interesting study as compared to classical order derivative. Because fractional calculus provides dynamical interpretations of real world phenomenon with more degree of freedom (see few as [50][51][52][53]). Therefore inspired from the mentioned applications and importance of fractional calculus and biological models, we update model (1) by incorporating natural death rates of various compartments as well as recruitment rate to form the following model where k is recruitment rate, d i ; i ¼ 1; 2; 3 denote natural death of susceptible, infected and recovered class respectively. Also the notation C 0 D b t stands for standard Caputo derivative of fractional order b 2 ð0; 1. Also l is recovery rate from disease. Here we present the flow chart of the model (2) in Fig. 1 as Here we will establish feasible region and verify boundedness of the solution for the corresponding model (2). Further using the next generation matrix method, we will develop global and local stability results also. Keeping the importance and real approach analysis of fractional calculus, we will investigate model (2) under formal arbitrary order derivative of Caputo type numerically. We establish some sufficient results for existence and uniqueness of solution to the model(2) by using some fixed point theorem due to Banach and Schauder. Some stability results about H-U type are also established for our problem. Further we will establish numerical algorithm based on NSFD scheme to simulate the results for various fractional order against the available real initial data. The concerned NSFD scheme has already used to handle many problems (see [54][55][56][57][58]). Also we will use modified Euler method (MEM) to describe the numerical results. The said method has also used in many papers for instance [50][51][52]. A comparison between both classical and fractional order model to conclude which one is more better for numerical interpretation of the proposed problem is given. We will use Matlab for our numerical analysis in this research work.
Here we organized our work as: In first part we give introduction. Then we give some elementary results in Section second. Third section is devoted to some fundamental results about the model under consideration. Fourth section is devoted to qualitative analysis. Fifth section is related to numerical analysis which is further divided into four subsections. Sixth section is devoted to comparison of both schemes. Further last section seven is related to a brief conclusion.
provided that integral on right exists pointwise on ð0; 1Þ.
( provided that integral on right exists pointwise on ð0; 1Þ.

Feasible region and stability theory
Here we derive some conditions for feasible region as well as for stability theory using the next generation matrix method.
Lemma 3.1. The solution ðS; I; R; DÞ 2 R 4 þ is bounded and attracted towards the feasible region defined by Proof. Summing all the equations of the Model (2) yields By solving (3), we get when t ! 1, then one has NðtÞ 6 k ðd 1 þd 2 þd 3 Þ . Hence the required result. h Theorem 3.2. The disease free equilibrium and pandemic equilibrium points of the model (2) are given by E 0 ¼ ð k d 1 ; 0; 0; 0Þ and such that Fig. 1 Flow chart of the proposed model (2).
Proof. For local and global stability analysis, the equilibrium points are important and can be computed from model (2) as Hence using (7) in model (2), one can easily computed the equilibrium points which are described in (5) and (6) (2) is computed as Proof. Let X ¼ I and considering second equation of (2) as Since R 0 is equivalent to the maximum eigenvalue of FV À1 at disease free equilibrium point E 0 ¼ ð k d 1 ; 0; 0; 0Þ. Therefore Hence the required quantity is given by Thus we arrived at our intended result.
Theorem 3.4. The proposed model (2) is locally asymptotically stable at E 0 if R 0 < 1, while the the model (2) is locally asymptotically stable at E Ã if R 0 > 1.
Proof. Since the jacobian matrix is computed from all four equations of the model (2), but here we take the first three equations of the model as these are independent of D. Therefore the Jacobian matrix for the model (2) can be computed as On setting the corresponding values, one has After putting the value of E 0 in (10), we have Now the characteristics equation can be computed as detðJ À gIÞ ¼ Thus the eigen values are given by Further g 2 can be written as Since we see that g 2 < 0, if R 0 < 1. Hence the needful results recived. h

Theoretical analysis
Here in this part of our work, fixed point theory is used to establish sufficient conditions for existence and puniness of solution to the considered model. Therefore, the existence of approximate solution and its uniqueness are investigated by using Schauder and Banach fixed point results [59,60]. We write our proposed model (2) with 0 < b 1 as Thank to Lemma 2.3, the system (11) is equivalent to the given system of nonlinear integral equations as For establishing theoretical results, we need the given hypothesis to be hold: ðE1Þ There exists constant L G > 0, for each UðtÞ; UðtÞ 2 R Â R, such that jGðt; UðtÞÞ À Gðt; UðtÞÞj 6 L G jUðtÞ À UðtÞj; ðE2Þ There exist constants C G > 0 and M G > 0, with jGðt; UðtÞÞj 6 C G jUj þ M G : If B : A ! A be the operator, then inview of (13), one has For any U 2 A, one has jBðUÞðtÞj 6 jU 0 j þ 1 CðbÞ R t 0 ðt À qÞ bÀ1 jGðq; UðqÞÞjdq; 6 jU 0 j þ 1 which implies that kBðUÞk 6 r: Therefore kBðUÞðt 2 Þ À BðUÞðt 1 Þk ! 0; as t 1 ! t 2 : Hence B is equi-continuous operator. Hence model (2) has at least one solution.
For existence of unique solution., we have the following result.
Cðbþ1Þ L G kU À Uk: ð18Þ Hence B is a contraction mapping, so by the use of Banach theorem, the considered system has a unique solution. h Let f 2 C½0; T with fð0Þ ¼ 0 independent of U as.

Numerical solution
This part is devoted to establish two numerical schemes. We apply these schemes to our model one by one and compared the results.

Numerical simulation by MEM
For the model (2), the numerical approximations are performed in this section by using MEM. Let ½0; T be the set of points, on which we must have to evaluate the series solution of the model (2). Upon further subdivision of ½0; T into m sub-intervals ½t p ; t pþ1 of equal difference j ¼ p m between consequent points using t P ¼ pj with p ¼ 0; 1; Á Á Á ; m, then obviously SðtÞ; IðtÞ; RðtÞ; DðtÞ t ½DðtÞ are continuous on ½0; T. Applying the modified Euler's or Taylor's method about t ¼ t 0 to the considered model expressed in (??) and for each value of t taking a 2 ð0; TÞ, the expressions for t 1 is given as Cð2bþ1Þ Á Á Á ; Cð2bþ1Þ Á Á Á ; Cð2bþ1Þ Á Á Á ; By taking the difference between consecutive points as j very very small, such that we may ignore terms containing higherorder derivatives to get On repeating the procedure, we get sequence of point to approximate ðSðtÞ; IðtÞ; ðRðtÞ; ðDðtÞÞ is formed. A generalized formula in this regard at t pþ1 ¼ t p þ h is given by where r ¼ 0; 1; 2; Á Á Á ; n À 1.

Numerical results and Discussion
Here we take real data of Pakistan as the total pollution of the country N ¼ 220 millions [61] and the other values are given in Table 1: Here we now present numerical interpretation of the proposed model (2) using MEM in Figs. 2-5. From Fig. 2, we see that susceptibility is decreasing at various rate due to different fractional orders. As fractional order tends to integer order the solution converges to the integer order solution. In same line the infected class is increasing in initial 100 days with various rate of transmission due to various fractional order then it turns to decrease with same behaviors as in Fig. 3. Consequently the increase in death class cause increase in recovered class. The recovery class dynamics raises at various fractional order as shown in Fig. 4. In Fig. 5, the dynamics of death class is also increasing until become stable due to faster infection rate.

Numerical interpretation of the Model (2) by NSFD scheme
Using NSFD scheme under the concept of fractional order derivative, we approximate the proposed model (2). We use Gru¨nwald-Letnikov approximation for Caputo derivative. About some detail for this scheme, we refer [54][55][56][57][58]. Since for nonlinear systems the investigation of exact solution in most cases is impossible, so we focuss on best approximation of the solution. In this regards, some NSFD schemes have been introduced. The said scheme has the ability to avoid the full implicit scheme and preserve positivity, monotonicity and con-vergency which make this method popular. For detail applications of the method see [62].
To construct the scheme, consider Gru¨nwald-Letnikov approximation for the fractional order derivative given by with t ¼ @h, such that h is the step size. Consider the following FDE as

Study of fractional order dynamics
Using (27), from (12) one has where t n ¼ nh and K b i are the Gru¨nwald-Letnikov coefficients calculated as and Based on the above definition, our considered model (2) can be discritized as

Numerical interpretation and explanation
Here, we now present the transmission dynamics of the proposed model (2) using NSFD scheme in Figs. 6-9. The dynamical behaviors presented in Figs. 6-9 is also same as given in Fig. 2-5 respectively.

Comparison of both methods
Here we compare graphs of both methods for fixed time t ¼ 300 in Figs Table 2.

Conclusion
In this research work, we have established a modified type COVID-19 model by incorporating natural death rates of susceptible, infected and recovered classes respectively. The model under investigation has been studied under fractional order derivative of Caputo type. Further by using fixed point theory, we have established existence theory for numerical solutions. Also we have developed various results for global and local stability by using the next generation matrix method. The basic reproduction number has been computed. Moreover, the feasi-ble region has also established for the proposed model along with its boundedness. The numerical interpretations have been performed by using two different numerical approaches based on MEM and NSFD methods. Both methods provide nearly same results for our model. Therefore we have compared both procedures in CPU time to see which one is most expensive with respect to time. In this regards NSFD method which is slightly general than MEM. Further MEM is slightly expensive in time than MEM for same number of iteration corresponding to different range of time. Graphical presentations have been provided for taking some real data about COVID À19 transmission in Pakistan. Both scheme have been compared     with real data. The concerned graphs have been plotted against different fractional order. Hence we concluded that fractional calculus provides more best explanations to real world problems for understanding their dynamics. In future, we will investigate some mathematical models under fractional order stochastic differential equations. Also piecewise concept will be applied for investigating epidemiological models.

Funding
There does not exist any funding source.

Availability of data
All data used in this paper is included within the article.

Authors contributions
All authors played their role equally. First and second authors designed the models and did theoretical analysis. Third author did numerical. Last two authors edited and drafted the article.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.