A New Numerical Scheme for Time Fractional Diffusive SEAIR Model with Non-Linear Incidence Rate: An Application to Computational Biology

: In this paper, we propose a modified fractional diffusive SEAIR epidemic model with a nonlinear incidence rate. A constructed model of fractional partial differential equations (PDEs) is more general than the corresponding model of fractional ordinary differential equations (ODEs). The Caputo fractional derivative is considered. Linear stability analysis of the disease-free equilibrium state of the epidemic model (ODEs) is presented by employing Routh – Hurwitz stability criteria. In order to solve this model, a fractional numerical scheme is proposed. The proposed scheme can be used to find conditions for obtaining positive solutions for diffusive epidemic models. The stability of the scheme is given, and convergence conditions are found for the system of the linearized diffusive fractional epidemic model. In addition to this, the deficiencies of accuracy and consistency in the nonstandard finite difference method are also underlined by comparing the results with the standard fractional scheme and the MATLAB built-in solver pdepe. The proposed scheme shows an advantage over the fractional nonstandard finite difference method in terms of accuracy. In addition, numerical results are supplied to evaluate the proposed scheme ’ s performance.


Introduction
Nowadays, the world faces the contagious outcomes of a fatal pathogen called coronavirus [1], which belongs to the family of SARS-CoV-2. Viruses of this family can transfer their pathogenicity to humans and animals. In humans, they affect the respiratory tract, injure the digestive lining, and cause severe gastroenteritis. Lung disorders may vary from mild cold to acute degenerative respiratory syndrome that ultimately becomes a cause of death. Precautionary and preventive measures have been taken at the government level to contain the virus to save the economy and health system from collapse. This situation of uncertainty between life and death has opened a new gospel in the field of quantitative mathematical modeling. Civil society and philanthropists are turning towards numerical modeling to gain insight to control the spread of the virus and predict possible measures to reduce cases. Many researchers have imparted fruitful already had the virus in the body [13]. Trend analysis of the coronavirus revealed that this disease has become endemic in some parts of the world, like Lassa fever in Nigeria [14] and Ebola in the Democratic Republic of the Congo [15]. This shows that if the infection is long forgotten, it may become endemic in a particular place. In this regard, epidemiologists examined whether certain infections are prevalent in society or completely lost [16][17][18]. They proposed that if 0 < 1, a disease can be eradicated from the public, whereas an epidemic occurs when 0 > 1 [19].
Alexander et al. [20] and Arino et al. [21] proposed that the concept of 0 < 1 is insufficient to describe the control of the spread of infection and introduced a backward bifurcation process. The backward bifurcation comprises three equilibriums for a certain range of parameters with 0 < 1. First is the stable disease-free equilibrium, second is a large stable endemic equilibrium, and third is a small unstable endemic equilibrium that acts as a boundary between two stable equilibriums. Sometimes, bifurcation leads to instability, making an infection endemic in the population, producing the first exponential outbreak. References [22][23][24][25] describe the bifurcation of COVID-19 models. This bifurcation phenomenon raised a question of epidemiology and risk for disease management since its existence implies that the basic reproduction number of the disease should be reduced to a value much lower than 1 to ensure the eradication of the epidemic. The classical epidemic model shows the direct relationship between the rate of treatment and the number of infected individuals. In the case of COVID-19, medical equipment such as oxygen machines, masks, patient beds, and healthcare staff are too dwindled compared to the number of patients. This imparts great pressure on the healthcare system. While developing a mathematical model for infectious diseases, such as the coronavirus, it is more appropriate to assume a saturation incidence treatment rate, which increases more slowly as the infected population grows larger.The present study comprises the SEAIR COVID-19 epidemic model considering ABC fractional derivatives due to the surged trend of their use in biological models.
The epidemic model presented in this research is based on five categories of susceptible, exposed, asymptomatic, infectious, and recovered individuals. Susceptible individuals do not have the pathogen. Exposed individuals have the pathogen but cannot infect others. The asymptomatic individuals have the pathogen, and these people spread the disease without being aware they are infected. The infectious individuals have the disease; these people spread the disease and have symptoms of the disease. The recovered individuals are those who are healed after having the disease.
A fractional model of susceptible-infected-recovered has been presented in [26]. The model considers the spread of infectious diseases. In addition, another model of the system of the incommensurate fractional differential equation has been used. The equilibrium points and their stability analysis have been computed and investigated. Another fractional model [27] for the treatment of cancer using stem cells and chemotherapy has been presented and local stability with respect to equilibrium points has been discussed. In the study [27], memory and hereditary traits were considered, which emphasized the advantages of the non-integer order model. Instead of a fractional model, a stochastic model comprising four human classes has been reformulated in [28]. It was shown that the solution of the proposed model existed and was unique. The stochastic Lyapunov function has been used for stationary distribution. To solve the model, a first-order stochastic Runge-Kutta method has been employed. In literature, various methods exist for solving fractional differential equations. Among these, the singular boundary method [29] has been utilized to solve variable-order time diffusion equations. The inclusion of the diffusion effect in fractional time differential equations has been described in [29]. Three types of methods have been employed to solve the problem, namely, the singular boundary method, dual reciprocity method, and finite difference method. For the discretization fractional time derivative, which was given in Caputo sense, the finite difference method was employed, and the remaining two methods were implemented for spatial discretization. The fundamental solution of the problem was transformed into inhomogeneous Helmholtz-type and then singular boundary method approximation was implemented and dual reciprocity method was considered to find the particular solution. Another numerical method was implemented in [30] to solve noninteger-order time reaction-sub diffusion equations, and the Riemann-Liouville derivative was considered. The method was based on the central difference approach for the discretization of the time derivative. For spatial discretization, hybridization of the cubic and Gaussian kernels was deliberated. The method consisted of developing a kernel that gave benefits that were achieved from the advantages of two different kernels. It avoided their limitations, and additionally, it maintained the global collocation method. The unconditional stability and convergence of the method have also been discussed. Another fractional advection-diffusion equation on the finite domain was considered in [31] and both time and space derivatives were non-integer order derivatives in Caputo sense. For the time-fractional direction, a quadratic interpolation was applied, and the spatial derivative was discretized by the Chebyshev collocation method of the fourth kind. The energy method was also employed to prove the unconditional stability and gain convergence order.
The standard susceptible-infected-recovered (SIR) model gives the interaction between susceptible and infectious individuals. This interaction term in the context of the Covid-19 pandemic can be studied with the safety measure of susceptible people and the quarantine of infected people. The effective interaction of susceptible and infectious individuals decreases by the increment of infectious people under some circumstances. A category of asymptomatic people describes the existence of a pathogen in a group of people which can spread the disease but these people are not aware of their disease and this category of people has also been considered in the literature.
The main aim of this work is to propose a numerical scheme for solving timedependent epidemic models. This scheme has the main advantage of providing a more accurate solution than the existing nonstandard finite difference method. The nonstandard method is named for its use of the nonstandard difference formula for time derivatives with a particular way of discretizing the differential equation that gives the positive and unconditionally stable solution. However, it has the drawback of providing an inaccurate solution, which is also shown in the presented simulations. In literature, there exist some ways to choose time step size in nonstandard finite difference method(s), but these ways do not describe the loss of the order of accuracy. it should be chosen so that the resulting difference equation must give at least first-order accuracy. Thus, there will not be more than one way to choose the time step size in the category of nonstandard finite difference methods when the scheme takes care of the order of accuracy and also provides a guarantee to get the positive solution. This is one of the issues in most of the existing nonstandard finite difference methods for addressing sufficient order of accuracy, due to having the main focus on developing nonstandard schemes for obtaining the positive solution without checking the order of accuracy. Obtaining the positive solution is a useful property in the category of models whose solution must be positive to make physical sense, but it must also achieve minimum first-order accuracy. In this work, the comparison of the proposed fractional numerical scheme is also made with the existing corresponding fractional nonstandard finite difference method.
Let be the infecting rate of infectious individuals, 1 the bilinear incidence rate, the transfer rate of exposed people to become asymptomatic, the transferring rate of asymptomatic people to become infectious, and the healing rate of both asymptomatic and infectious people. It is assumed that the exposed people have the pathogen but cannot spread it and that the asymptomatic people have the pathogen and can spread it, but are not aware that they have the disease. The infectious people have pathogens, they are aware that they have the disease, and they can spread the disease. The final category of people is recovered people, which is the class of people who are healed from the disease. The proposed epidemic model can be expressed as: One of the trivial solutions of the system (1)-(5) is the disease-free state (1, 0, 0, 0, 0). It corresponds to a situation when the pathogen is absent.

Theorem 2. The system (1)-(5) is in linear stability of the disease-free equilibrium states if it
Proof. To prove this theorem, the Jacobian matrix of system (1)-(5) is expressed as: At disease-free equilibrium point (1,0,0,0,0), the Jacobian can be expressed as: The eigenvalue of Jacobian | can be found from the following equation: where . is an identity matrix of 5 × 5.
∎ The model given in (1)-(5) is a set of ordinary differential equations. To propose a more general model, a diffusion effect can be considered. Since the model of ordinary differential equations is constructed only on time variables, the model of partial differential equations is given constructed on both time and space coordinates. The diffusion effects describe the spread of susceptible, exposed, asymptomatic, infected, and recovered people over the spatial coordinate. If a diffusive space term is introduced, the fractional diffusive epidemic model can be expressed as: subject to the boundary conditions, and the initial conditions are expressed as: where 1 , 2 , 3 , 4 , and 5 are constants.

Proposed Fractional Scheme
To solve Equations (12)-(18), a fractional numerical scheme is proposed based on the fractional Taylor series approach. The first Equation (12) in the diffusive model is discretized as: where 1 is unknown, to be determined later. We expand +1 and −1 using fractional Taylor series as: Substituting the fractional Taylor series expansion for +1 and −1 from (20) and (21) into Equation (19) results in: By substituting Equation (12) into Equation (22), we obtain: Comparing coefficients of and on both sides of Equation (23) results in: . Solving Equation (24) gives 1 : The difference equation for numerically solving Equation (12) can be expressed as: Similarly, the difference equations for solving the remaining Equations (13)-(16) are expressed as: . Proof. We rewrite Equation (26) as: This is positive because every term is positive.
Similarly, we rewrite Equations (27)-(30) as: Since all terms on the left-hand side of Equations (32)- (35) are positive, the components of the solution at the th grid point and at the ( + 1)th time level will be positive. Since the same difference Equations (31)-(35) will be used to find each component of each solution at every grid point and at every time level, the positive solution will be obtained subject to the satisfaction of conditions given in the statement of the theorem. Therefore, the proposed scheme is a conditionally positivity-preserving scheme.

Stability of Proposed Fractional Scheme
In this section, the stability of the proposed scheme for the first fractional equation in the considered model (12)-(16) is presented. Since the first Equation (12) is a non-linear fractional differential equation, it is linearized as: where = 1+ 1 2 will be considered as a fixed quantity. At the beginning of the stability analysis procedure, the difference equation for solving Equation (36) is constructed first, which can be expressed as: According to von Neumann stability analysis, the following transformations are considered: where = √−1.
Substituting transformations (38) into Equation (37) results in: Dividing both sides of Equation (39) by , we obtain the following: We rewrite Equation (40) in the form: Since the scheme is constructed at three-time levels, one more equation is formed. For this reason, we write as: We rewrite Equations (42) and (43) in a matrix-vector equation as: The stability conditions will be obtained by imposing a condition on the amplification matrix of Equation (44). These conditions are imposed on the eigenvalues of the matrix given in Equation (44). Let 1 and 2 be the two eigenvalues of the amplification matrix. The stability conditions are expressed as: . Therefore, if satisfies condition (45), the scheme will be stable. This can be verified by using maximum or minimum values of and some value of . By doing so, two functions in terms of can be obtained. These two functions can be plotted over different ranges of parameter , and two-dimensional graphs of the eigenvalues can be obtained. On these graphs, the eigenvalues should lie within the closed interval [−1,1] for any range of .

Consistency of Fractional Proposed Scheme
The consistency of the proposed scheme will be checked for the model of fractional differential equations. Equation (19) is the discretized equation for Equation (12), so to check the consistency of the scheme, we expand +1 and −1 in a fractional Taylor series form, which already expanded in Equations (20) and (21). The next step is to expand the terms +1 and −1 in the Taylor series as: Adding +1 and −1 , the following is obtained: Substituting expansions for +1 and −1 from (20) and (21) and (48) into Equation (19) results in: We rewrite Equation (49) as: By substituting Equation (25) in Equation (50), we obtain: Equation (51) can be simplified as: Applying ( ) → 0, → 0, the original Equation (12) is obtained, evaluated at the ℎ grid point and the th time level. Thus, the proposed scheme is consistent with Equation (12). Similarly, it can be proved that the proposed fractional scheme is consistent for the remaining Equations (13)-(16).

Issue of Accuracy in Fractional NSFD
The non-standard finite difference method (conventional-NSFD) is one of the existing numerical methods providing unconditionally positivity-preserving solutions and unconditionally stable solutions for epidemic models. However, it has the drawback of lacking first-order accuracy. In this work, it is pointed out that fractional NSFD is not even first-order accurate in fractional time. In order to provide theoretical evidence, the fractional Taylor series is implemented on a difference equation, obtained by applying fractional NSFD on any of the equations in the fractional diffusive epidemic models (12)- (16). To determine its lack of accuracy, we consider the difference equation obtained by discretizing Equation (12) using the fractional non-standard finite difference method as: Since expansions for fractional Taylor series of +1 are already given in this work, we substitute the expansion for +1 from (20) into Equation (53) and obtain the following: Substituting Equation (12) into Equation (54) evaluated at the ℎ grid point and ℎ time level gives: Comparing coefficients of both sides of Equation (55) results in: From Equation (56), it can be concluded that the coefficients of first-order fractional time derivative terms do not vanish, so fractional NSFD is not first-order accurate in time, according to the fractional Taylor series approach considered in this work. Similarly, it can be proved that fractional NSFD is not first-order accurate in time for Equations (13)-(16).

Convergence of Proposed Fractional Scheme
For finding the convergence condition(s) of the proposed scheme for the system of fractional Equations (12)-(16), a matrix-vector equation is formed. Since both Equations (12) and (13) are non-linear, the Jacobian evaluated at disease-free equilibrium point (7) is calculated. Since the Jacobian contains both positive and negative entries, it is divided into two matrices that contain all positive entries. The Jacobian can be expressed as: We discretize Equation (61) using the proposed scheme and obtain: The unknown will be found similarly as in the construction of the fractional proposed scheme, which is given as: .

Results and Discussions
The numerical experiments show the effectiveness of the proposed fractional scheme. Three numerical schemes are employed to find numerical solutions of considered diffusive fractional models. The MATLAB built-in solver pdepe, which can be used to find solutions to parabolic and elliptic equations, is employed to check the accuracy of the schemes with integer-order time derivatives. For the integer-order time derivative, the considered diffusive equations become parabolic equations with source terms. Thus, these types of linear and non-linear equations can also be solved by the MATLAB built-in solver pdepe. The solver can be implemented for the cases that correspond to the positive solutions. If the solutions of the considered model are positive for particular values of the contained parameters, any high-order scheme can be applied to obtain solutions to the problem(s). The MATLAB built-in solver is used for integer-order derivatives, and for fractional derivatives, the first-order standard Euler scheme is considered. Since the nonstandard finite difference method considered in this work is explicit, the proposed and first-order forward Euler schemes are also explicit, so these schemes do not require the use of an iterative scheme to solve the discretized or difference equations. However, to treat the boundary conditions, an additional iterative scheme is required. To handle boundary conditions, first-order forward and backward Euler formulas are employed. The boundary condition for susceptible individuals at the left endpoint of the domain is discretized as: Since the information of susceptible individuals is required at the first grid point and at ( + 1) ℎ time level, Equation (76) can be expressed as: However, the information of 2 +1 is not given, so 1 +1 cannot be found explicitly. Therefore, an additional iterative scheme is employed at the first grid point as: In Equation (78), in the superscript is used for the iteration number. Since an initial guess is chosen at the first iteration, the information of the chosen initial guess will be used to find susceptible individuals at the first grid point. These iterations will continue until the stopping criteria are met. The stopping criteria are based on the norm of the difference between solutions computed on two consecutive iterations. When the norms of each equation are less than some given smaller value close to zero, the iterations will be stopped, and the final solution will be obtained. Figures 1-3 show the comparison of the proposed and non-standard finite difference methods. The scheme will be more accurate if the absolute error in Figures 1-3 is near zero. The proposed scheme produced a smaller error as compared to NSFD over time with an integer-order time derivative. The error is determined by finding the absolute difference between the solutions obtained by NSFD/proposed scheme and MATLAB built-in solver pdepe. For non-integer time derivatives, the comparison is made with the first-order fractional forward Euler method. The three-dimensional error plots are given in Figures 4-8. Again, the scheme will be more accurate if the absolute error is less than the absolute error produced by the other scheme. From Figures 4-8, it can be seen that the proposed fractional numerical scheme is more accurate than the fractional non-standard finite difference method. Moreover, these results are numerical evidence that the non-standard finite difference method considered in this work is not first-order accurate.        Figure 9 compares the impact of the bilinear incidence rate parameter 1 on susceptible, exposed, asymptomatic, and infectious individuals. This figure shows that susceptible individuals de-escalate by increasing values of bilinear incidence rate, and exposed, asymptomatic, and infectious people have both increasing and decreasing behavior. These three categories of people have increasing behavior for the first interval of time, and then this behavior is reversed. Figure 10 shows the contour plots for recovered people by varying bilinear incidence rate parameters. The difference between both contour plots can be seen in this figure. Figure 11 shows the impact of fractional parameter on four categories of people. Clearly, Figure 11 shows the increasing behavior of susceptible, exposed, asymptomatic, and infectious people by enhancing the values of the fractional parameter. By increasing the order of fractional time derivatives, a growth in four categories of people can be seen. Figure 12 shows the impact of fractional parameters on recovered people in the form of contours. Figure 9. Impact of the bilinear incidence rate parameter 1 on susceptible, exposed, asymptomatic, and infected individuals using = 50, = 300, = 0.9, = 0.7, = 0.9, = 0.9, 1 = 2, 2 = 2.5, 3 = 2.5, 4 = 2, 5 = 2, = 0.9, = 30.  Figure 10. Impact of the bilinear incidence rate parameter 1 on recovered individuals using = 50, = 300, = 0.9, = 0.7, = 0.9, = 0.9, 1 = 2, 2 = 2.5, 3 = 2.5, 4 = 2, 5 = 2, = 0.9. Figure 11. Impact of fractional order parameter on susceptible, exposed, asymptomatic, and infected individuals using = 50, = 300, 1 = 0.9, = 0.9, = 0.7, = 0.9, = 0.9, 1 = 2, 2 = 2.5, 3 = 2.5, 4 = 2, 5 = 2, = 30.  Table 1 shows the comparison of the proposed scheme with the nonstandard finite difference method for the integer-order time derivative. The errors produced by both schemes can be seen in this table. The error is found when the norm of the difference of two solutions is computed. Both schemes are checked when their comparison is made with MATLAB solver pdepe. The higher norm of error shows the higher error in the solution. In Figures 9 and 11, denotes the length of the boundary. Figures 9 and 11 are drawn at a fixed value of = 14.6939 and in Table 1, Nx = 50, and Nt = 300 denote the number of grid points and the number of time levels respectively. Table 1. Comparison of proposed and NSFD schemes over spatial coordinate x using d1 = 2, d2 = 2.5, d3 = 2.5, d4 = 2, d5 = 2, Nx = 50, and Nt = 300. It is to be noted that the values of parameters are chosen randomly and are somewhat suitable for the comparison of the two schemes.

Conclusions
In this work, a fractional numerical scheme is proposed to solve fractional timedependent problems. A modified model of the diffusive epidemic model was constructed with a non-linear incidence rate. The proposed scheme was second-order accurate in space. The stability of the fractional numerical scheme has been given, and convergence conditions have also been found. It was also pointed out that the fractional non-standard finite difference method is not accurate enough to use in its present form. The performance of the proposed scheme was also checked by comparing the results with the MATLAB built-in solver pdepe for the classical case. This comparison showed the accuracy of the proposed scheme and described the inaccuracy of the considered form of the nonstandard finite difference method for problems with diffusion effects. The proposed fractional numerical scheme can be used in situations where the conditions of obtaining a positive solution are desired with unconditional stability for numerous problems in science and engineering. As a result of the completion of this study, it is feasible to propose further applications for the currently available approach [34][35][36].

Data Availability Statement:
The manuscript includes all required data and implementing information.