A SARS-CoV-2 Fractional-Order Mathematical Model via the Modiﬁed Euler Method

: This article develops a within-host viral kinetics model of SARS-CoV-2 under the Caputo fractional-order operator. We prove the results of the solution’s existence and uniqueness by using the Banach mapping contraction principle. Using the next-generation matrix method, we obtain the basic reproduction number. We analyze the model’s endemic and disease-free equilibrium points for local and global stability. Furthermore, we ﬁnd approximate solutions for the non-linear fractional model using the Modiﬁed Euler Method (MEM). To support analytical ﬁndings, numerical simulations are carried out.


Introduction
COVID-19 caused by SARS-CoV-2 is a highly transmissible and pathogenic coronavirus that appeared in late 2019 and has posed a danger to both human health and public safety [1,2]. The World Health Organization (WHO) declared it an endemic because it killed and affected many worldwide, particularly in the USA and Europe, in a short period of time [3]. Clinical research has not yet produced a remedy that can completely eradicate the virus from the human body. However, the research remains ongoing. Many therapies (such as vaccination, monoclonal antibody therapy, and plasma therapy) have been developed by researchers who successfully treat early-stage diseases, including MERS-CoV, SARS-CoV, Ebola, HIV, and influenza-like viral diseases. Additionally, our body's immune system responds well to infections or disorders. Throughout this uncontrollable condition, nearly all virus-infected regions were locked down, social gatherings were banned, and strict social distancing measures in all situations were implemented to control virus spread. Research worked to control the spread of this virus from different points of analysis such as microbiology [4,5], pathology [6,7], and applied mathematics [8].
In addition to biological and medical research, theoretical studies based on mathematical models may also be crucial to this anti-pandemic effort in analyzing the behavior of the outbreak that made it an epidemic, making decisions about how to stop its spread, and in comprehending patterns of virus transmission within hosts. Several mathematical models for COVID- 19 have been formulated at an epidemiological level [6,[9][10][11][12][13][14][15][16][17][18][19][20][21]. The replication cycle of the SARS-CoV-2 virus and its interactions with the innate and adaptive immune systems are only limited by research done at the within-host level [5,22,23]. Susceptible, infected, and removed (SIR) models are generally used for measles, rubella, and other infectious diseases to investigate the quick spread behavior of these infectious diseases [24]. A susceptible, exposed, infected, and removed (SEIR) model is comparable to the SIR model. The classes S, I, and R represent the number of populations in each partition at a specific point in time [25,26]. However, the incubation period is used in the SEIR model, so a specific incubation period for infectious diseases is more applicable [26,27]. The incubation period of SARS-CoV-2 is 1-14 days [28][29][30][31]. Different models have been developed to study the mechanism behind the spread of COVID-19 [32][33][34].
At the time of SARS-CoV-2 infection, macrophages are first targeted by SARS-CoV-2, and SARS-CoV-2 propagates to T cells afterwards. At this stage of the virus route, T cells are activated and further involve differentiation. In addition, T cells produce cytokines (INF-α), IL-6, and IL-10) associated with the different types of a T cell. A large amount of cytokines provides a greater activation of the immune response to fight against the virus. Particularly, T cells, CD4 + T cells, and CD8 + T cells have been playing a significant antiviral role in the fight against pathogens. There is also a risk of mounting autoimmunity or devastating inflammation. CD4 + T cells help the immune system of the body by generating virusspecific antibodies with the activation of T-dependent B cells. However, CD8 + T cells can kill virally infected cells, as they are cytotoxic. In general, many CD8 + T cells in the infected SARS-CoV-2 body are found in nearly 80% of the total infiltrative inflammatory cells in the interstitial pulmonary tract, which play a significant role in clearing CoVs. The loss of CD4 + T cells is correlated with the reduced conscription of lymphocytes and neutralizing the production of antibodies and cytokines, resulting in severe immune-mediated interstitial pneumonitis and delayed SARS-CoV-2 lung clearance [34][35][36][37].
Researchers have shown that a long-lasting and persistent response of T cells to S and other structural proteins (including the proteins M and N), which provides the sufficient knowledge to draft the SARS vaccine by combining viral structural proteins. These vaccines may provide a robust, efficient, and long-term response to the virus by memory cells [36]. Clinical trials also show that a monoclinic antibody therapy is an effective treatment tool that responds well to SARS-CoV-2 [34].
Fractional calculus is often employed for epidemic models [38][39][40][41][42][43][44]. It has been shown that the fractional-order model outperforms the integer-order model in handling the modeling process since it has many other desirable characteristics, such as excellent fitting with real data. Furthermore, memory and heredity characteristics make it more effective in modeling and analyzing real-world problems. Numerous definitions or operators in the fractional calculus, such as Caputo, Riemann-Liouville, Atangana-Baleanu, and Caputo-Fabrizio derivative [45,46], are helpful in modeling epidemic diseases.
In our paper, we study a model of in-host viral kinetics that describes the response of SARS-CoV-2 to epithelial cells. Section 2 presents the formulation of the model, Section 3 describes preliminaries, Section 4 presents the uniqueness and existence of the solution to the proposed model, Section 5 presents the steady state points and the derivation of the primary reproduction ratio, Sections 6 and 7 present local and global stability analyses, respectively, Section 8 presents the numerical solution algorithm, Section 9 shows graphical representations to support the logical result, and Section 10 concludes.

Mathematical Model
We consider the mathematical model [5,47] used to analyze the fractional model of the within-host viral kinetics of SARS-CoV-2.
Here, the whole population is divided into three compartments: virus-free pulmonary epithelial cells, represented by E p (t), virus-free cells, represented by F(t), and virusinfected pulmonary epithelial cells, denoted by I p (t). Here δ, β, and ω denote the death rate of uninfected pulmonary epithelial cells, infected pulmonary epithelial cells, and the virus, respectively. The production rate of the virus is denoted by ρ. In the model, the term E p (0) denotes the initial number of uninfected cells at time zero, and δE p (0) shows a constant regeneration of uninfected epithelial cells.
Several researchers have taken into account fractional-order derivative (FOD) for modelling the infectious diseases [48][49][50]. We used FOD in the Caputo sense, because it gives better results than the classical order. Several FOD operators were developed, but the Caputo and the Riemann-Liouville FD operators are the most generally used due to their simplicity and similarities [38,39,51]. Other derivatives include Katugampola, Hadamard, Atangana-Baleanu, and Caputo-Fabrizio [52].

Existences and Uniqueness
Firstly, we discuss the uniqueness and existence of the solution to the proposed model.
Applying the fractional integral as Using the initial conditions, (5) becomes We define the kernels in System (6) as Therefore, by using (7), the system (3) becomes Proof. By considering the kernels given in (3), we have Hence, the proposed system satisfies the Lipschitz continuity. This means that the system is continuous and bounded.
and it is also a closed set. We now define the operator T in F as We then have Hence, Similarly, we easily show that Finally, we need to show that T satisfies the contraction mapping theorem.
Assuming that E p (0) = E * p (0), which implies that By using the value of |( Hence, Since A 1 C 1 M 1 < 1, one can conclude that T satisfies the contraction mapping theorem and it has a unique fixed point.
In the same procedure, we obtain Thus, System (3) has a unique solution.

Steady State and Derivation of Reproduction Number R 0
Let E 0 be the disease-free equilibrium point, and from the equations in (3), we set the right-hand sides as equal to zero and solve for variables.
If I p and F are equal to zero, we obtain the disease-free equilibrium as To find the endemic equilibrium E * = E * p , I * p , F * , we consider E p (t) = 0, I p (t) = 0, in Equations (13)-(15) as Next, we want to derive the reproduction number R 0 . The reproduction number is a threshold quantity representing the total number of secondary illnesses caused by an infected individual in a fully susceptible population throughout the infectious period [53].
For this, we develop the next generation matrix approach. Let M(S ) represent the rate of new infections, and let X(S ) represent the transfer rate of individual. We then have The Jacobian matrices J M(S ) and JX(S ) at the non-infected steady state (16) are The reproduction number is given by the spectral radius of the next generation matrix mx −1 . We then obtain λ 1 = 0, λ 2 = ρ α γ α E p (0) λ α β α , from the absolute highest eigenvalues from mx −1 . Thus, from the model (3), we obtain the expression for R 0 as The following graphical representation describes the potential impact of introducing a new SARS-CoV-2 variant with high transmission rates. Figure 1 represents the effects of the basic reproduction number with respect to different parameters.

Local Stability
In this part of the paper, we discuss the local stability of the steady-state point (16). For this need, we write the Jacobian matrix for the system (3) as Theorem 4. The non-infected steady state E 0 is locally asymptotically stable when R 0 < 1.
Proof. The Jacobian matrix for the steady state (16) as According to the Routh-Hurwitz criteria [54], if all the real parts of the value of the real eigenvalues of the Jacobian matrix are negative, then E 0 is locally asymptotically stable. In order to show that, we have The characteristics equation of the above matrix is with a 1 = ω α + β α , and a 2 = β α ω α − ρ α γ α E p (0). Based on an equation, the eigenvalues λ 1 = δ α > 0, and a 1 = ω α + β α > 0 were obtained. According to Routh-Hurwitz stability criteria, β α ω α ⇒ R 0 < 1. Therefore, if R 0 < 1, then all the conditions of Routh-Hurwitz criteria are satisfied. Hence, the non-infected steady-state E 0 becomes locally asymptotically stable if R 0 < 1 and unstable if R 0 > 1.

Global Stability
In this part of the paper, we use the Lyapunov function approach to investigate the result for the model's globally asymptotically stability (GAS) in the disease-free state. For this, we have the following theorem.
Theorem 5. When R 0 < 1, the disease-free equilibrium E 0 is globally asymptotically stable; otherwise, it is unstable.
Proof. For this, we derive the Lyapunov candidate function (LCF) for the fractional-order as in [38]. We then take the family of the Lyapunov function L(e 1 , e 2 , ..., e n ) = LCF is defined as Applying the Caputo fractional derivative and using its linearity property yields Applying Lemma 2 of [38] to (23), we obtain Substituting the disease-free equilibrium, we obtain where Hence, by Theorem 2 in [39], the disease-free equilibrium is globally asymptotically stable.
The existence of equilibria of the model is shown for different values of R 0 . In plotting Figures 2 and 3, we have varied the value of infection rate γ. It is observed that the diseasefree equilibrium E 0 is stable for low infection rates, corresponding to R 0 < 1. It becomes unstable for high infection rates, corresponding to R 0 > 1.

Algorithm for the Solution
In this part of the paper, we derive the general procedure for the approximate solution of the considered fractional-order model (3). For this, we extend the numerical Euler method [55]. Now, using (9), we write (3) as Let [0, T] be the interval where we want to find the solution of (27). For this, we divide the given interval into an i sub-interval [t i , t i+1 ] of equal width h = T i , using the nodes t i = ih for i = 0, 1, 2, 3, ..., k − 1. Here, we assume that E p , c 0 D α t [E p ], and c 0 D 2α t [E p ] are continuous up to a higher order on [0, T].
We take into account the modified Euler method (MEM) for t = t 0 = 0 to (27). The expression for t 1 is In this procedure, we neglect the second-order involving t 2α because we the step size h is very small.
The subsequent terms are In the above fashion, a general formula at t i+1 = t i + h is developed: where i = 0, 1, 2, ..., k − 1.

Graphical Representations
In this section, we discuss the obtain approximate solutions of (31) graphically using the values of the parameters shown in Table 1 with different α values. We conduct graphical representation to represent the solution effects of our model with different fractional-order α. By using MATLAB software, we set up an algorithm to simulate the results in Figure 4-9. For this, we consider some appropriate values used in the model for the parameters in Table 1. The parameter values of the model (3) have been estimated based on chest radiograph score data [5] using the Monte Carlo Markov Chain method. We utilize that parameter; more specifically, E p (0) = 50, I p (0) = 50, and F(0) = 100.
In Figure 4, we take the fractional-order α = 0.9, which is close to the integer-order α = 1. For this model, the results obtained from the fractional-order model of SARS-CoV-2 are very similar to those of the integer-order model [5]. Further, Figure 5-9 show series plots of the approximate solutions (31) with the fractional orders slowly decreasing as follows: α = 0.90, 0.80, 0.70, 0.60, 0.50. We observe from these graphs that the virus-infected, virusfree epithelial and the virus-free curves have approximately the same trend for different values of α. However, their convergence toward stability of the virus-infected epithelial is slightly changed. There is a notable difference in the graph at order 0.90 and 0.50. We also note that the equilibrium points are the same for the different fractional orders in Figure 4-9, but the equilibrium points of the virus-infected epithelial go to the fixed point over a longer period when the fractional-order increases.      From the above plotting, we have concluded that the infection rate γ can be low. This is because, if we take the greater infection rate in such a way that R 0 > 1, the infection prevails in the body since the virus cannot be controlled. However, if the infection rate γ is lower than R 0 < 1, this leads to being controlled by the virus-infected epithelial cells rate in the body. Otherwise, the number of virus-free epithelial cells will decrease. Such infected virus epithelial cells and the virus-free cells will not stabilize very quickly in a body with time. In this case, the patient's condition will vitiate with time. We also conclude that the fractional derivative is the generalization of the integer-order derivative. We use fractional-order to obtain good accuracy and reduce the possible errors arising from the mistreated parameters in the distributed modeling system and the system with memory.

Conclusions
This article presents a Caputo fractional-order model for the within-host dynamics of SARS-CoV-2. First, we found the existence and uniqueness of the model's solution by using the fixed point theory. We obtained the disease-free and endemic equilibrium points. The critical parameter (reproduction number) R 0 was determined using the next-generation matrix approach. Moreover, we developed the considered model's local and global stability conditions. Further, we calculated the approximate solution to the model (3) via a powerful method due to Euler. Finally, graphical representations of the obtained approximate solutions were performed to verify the theoretical analysis using MATLAB. By biology, an infection will be cleaned from a human body when the basic reproduction number R 0 < 1 is approached; otherwise, therapy must be used to minimize and eliminate the infection from the body. Additionally, it was shown that an infection will continue for any level of viral load in the host's body if the basic reproduction number is greater than one. For the eradication of the virus from an infected human body, our future work will include using various treatment strategies and control parameters in this model. Funding: This research received no specific grant from any funding agency in the public, commercial, or not-for-profit sectors.