Analysis of a fractional model for HIV CD4 T-cells with treatment under generalized Caputo fractional derivative

Abstract: In this paper, a mathematical model of generalized fractional-order is constructed to study the problem of human immunodeficiency virus (HIV) infection of CD4 T-cells with treatment. The model consists of a system of four nonlinear differential equations under the generalized Caputo fractional derivative sense. The existence results for the fractional-order HIV model are investigated via Banach’s and Leray-Schauder nonlinear alternative fixed point theorems. Further, we also established different types of Ulam’s stability results for the proposed model. The effective numerical scheme so-called predictor-corrector algorithm has been employed to analyze the approximated solution and dynamical behavior of the model under consideration. It is worth noting that, unlike many discusses recently conducted, dimensional consistency has been taken into account during the fractionalization process of the classical model.


Introduction
The human immunodeficiency virus (HIV) is a kind of retrovirus and it is the cause of acquired immunodeficiency syndrome (AIDS) [1]. CD4 + T-cells are targeted of HIV, which are the most abundant white blood cells of the immune system and play a major role in protecting the body from infection. When HIV infection occurs, the virus attacks the CD4 + T-cells and tries to replicate itself. The process of infection completes is called the HIV life cycle [2], which carried out in five stages: Stage I: (binding and fusion) when HIV is entering into the CD4 + T-cell, it is detected and bound by receptors at the cell membrane and then it has been fused, which allows HIV to enter inside the CD4 + T-cell. Stage II: (reverse transcription) within the CD4 + T-cell, HIV releases and uses reverse transcriptase enzyme to convert its genetic material RNA to DNA, which allows HIV to enter the CD4 + T-cell nucleus. Stage III: (integration and transcription) within the nucleus of CD4 + T-cell, HIV releases and uses integrase enzyme to integrate its viral DNA into DNA of the CD4 + T-cell, this is called proviral DNA. Afterward, the cell is motivated, it will turn on the transcription of the proviral DNA to produce the viral RNA and move out of the nucleus to the cytoplasm. Stage IV: (replication and assembly) within the cytoplasm, the viral RNA begins to use the machinery of the CD4 + T-cell to make long chains of HIV proteins, which are the building blocks for HIV replication. After that HIV protein and RNA move to the surface of the cell membrane and assemble into immature HIV. Stage V: (budding) the immature HIV pushes itself out of the host cell and releases protease enzyme to split the long protein chains in the immature virus. Lastly, it becomes a mature virus which makes effects on other CD4 + T-cells. Antiretroviral therapy is the use of HIV medicines to treat HIV infection and to protect the immune system by blocking the virus from reproducing at different stages of the HIV life cycle for instances fusion inhibitors, reverse transcription inhibitors, integrase inhibitors, and protease inhibitors [3][4][5].
Mathematical modeling and simulations were used as an essential instrument to forecast the possibility and severity of diseases and allow information to know its dynamic behavior of the infection. These models are useful tools and play a significant role to explain the dynamics of the immunological response to HIV infection, we refer readers to some books, see. [6][7][8]. Presenting some important previous works about the HIV model, researchers commonly used mathematical expressions in the form of differential equations to establish the relation between HIV viruses and uninfected CD4 + cells, and the effect of drug therapy on infected cells is described by mathematical models. In 1989, Perelson [9] initiated a simple model for the HIV infection and analyzed some dynamics behavior of the HIV infection of CD4 + T-cells in 1993 with his team [10]. In 2004, Tuckwell et al. [11] studied the dynamics behavior of solutions in the two-component model of the HIV infection model as follows: where S is denoted the uninfected CD4 + T-cells, I is denoted the density of infected such cells, λ is constant rate, µ 1 is their per capita death rate, µ 2 is the rate of disappearance of infected cells, and k is infection rate. After three years, Rong et al. [12] developed a mathematical model of HIV model to predicate the initial constraints that shape the evolution of viral resistance to antiretroviral drugs in 2007 as: where they considered three compartments of CD4 + T-cells including concentration of uninfected (susceptible) cells S , infected cells I in eclipse phase, and free HIV virus particles in the blood cells V. Here denotes δ death rate of infected S cells and includes the possibility of death by bursting of infected S cells, hence δ ≥ µ 1 . Further, the parameter b is the rate at which infected cells return to the uninfected class, while c is the death rate of the virus and N is the average number of viral particles produced by an infected cell, where k is the rate of infection. In 2009, Srivastava et al. [13] presented a primary infection model with reverse transcriptase (RT) inhibitor which takes place before the infected cell starts producing the virus. They studied the dynamics of the uninfected CD4 + T-cells, the virus population, and the infected CD4 + T-cells. They also distinguished the infected class of CD4 + T-cells into two subclasses: Pre-RT class which is the phase of the infected cells in which the reverse transcription is not completed and post-RT class which is the phase of those infected cells in which the reverse transcription is completed. The model of integer order is shown of the form: with the initial conditions (S (0), I(0), V(0), L(0)) T = (S 0 , 0, 0, L 0 ) T , where (·) T represents the transpose operation, S , I, V and L are denoted the density of susceptible CD4 + T-cells, the density of infected CD4 + T-cells in pre-RT class, the density of infected CD4 + T-cells in post-RT class and the density of virus, respectively. The stability and numerical results are obtained in their work. The unknown parameters will be provide in the details in Section 3. The researchers used differential equations to applied in various areas to real-world problems. Since the theory of differential equations is well-known for applied researchers. It has popular in many fields and is working on different aspects of the mentioned area. In the last three centuries, various researchers developed fractional-order differential equations. It was extensively applied in various areas to real-world problems; see [14][15][16]. The well-known aspects that has been established recently are existence theory, stability analysis, numerical and optimization methods. The mathematician researchers introduced the different fractional derivatives that have been provided in fractional calculus. There are many types of research studies that have been created by using fractional derivatives and apply them to the HIV model. For examples; In 2009, Ding and Ye [17] studied the stability of equilibrium for a fractional HIV infection of CD4 + T-cells model and numerical simulations presented. In 2011, Ertürk et al. [18] discussed an approximated solution of a fractional HTLV-I infection of CD4 + T-cells model using a multi-step generalized differential transform method. In 2012, Arafa et al. [19] constructed and analyzed fractional modeling dynamics of HIV and CD4 + T-cells during primary infection. The generalized Euler method (GEM) used to find a numerical solution of this model. A few years later, Arafa et al. [20] studied and analyzed the behavior of solutions for a fractional-order model of HIV infection with CD4 + T-cells in Caputo sense by using the GEM. In 2017, Arshad et al. [21] analyzed the effects of the HIV infection model on CD4 + T-cell population based on a Caputo fractional derivative. By applying a finite difference approximation, the numerical scheme has been proposed. Khan et al. [22] established the stability results and numerical solutions of the fractional HIV/AIDS model in sense of an Atangana-Baleanu derivative operator. In 2019, Lichae et al. [23] studied an HIV-1 infection model of CD4 + T-cells with a description of the effect of antiviral drug treatment under a Caputo fractional derivative. The approximate solution has been gained by applying the Laplace transform and the Adomian decomposition method. Ferrari et al. [24] presented an HIV model under Caputo fractional derivative, which the presences of a reverse transcriptase inhibitor based on [13] is indicated. They proved the existence and uniqueness of the model and shown that the model is positive invariance. The equilibrium points and their stability of the model are also investigated. Recently, the qualitative theory has been interested and applied to mathematical models. For example, in 2020, Khan et al. [25] studied and analyzed the existence and stability results for the fractional order HIV-TB model with the Mittag-Leffler kernel. Further, numerical solutions are obtained. The HIV model under the Caputo-Fabrizio fractional derivative is studied by Nazir et al. [26]. They obtained the existence criterias of the solutions by using fixed point theory. The stability of the concerned solution is also showed in the Hyers-Ulam sense. The approximate solution of the model is verified by using the technique of Laplace transform with the Adomian decomposition method. Abdeljawad et al. [27] studied the existence of solutions of the q-fractional difference equations with help of Krasnoselskii's fixed point theorem and show the results by the Lotka-Volterra model.
Abdeljawad and co-workers [28] applied the Banach fixed point theorem to prove existence results for the discrete fractional differences in sense of the discrete Mittage-Leffler kernel and numerical results were provided for the presence of the results. Khan et al. [29] discussed the existence of solutions and stability for the fractional-order advection-diffusion model in sense of the Mittag-Leffler kernel and shown the results numerically. Khan et al. [30] studied the time-fractional COVID-19 model with the help of the corrector-predictor method in the sense of a generalized Caputo fractional derivative [31]. The existence results are investigated by using Schauders and Weissingers fixed point theorems. Many papers have been developed under various types of fractional derivatives, see [32][33][34][35][36][37][38][39][40][41][42][43][44][45] and references cited therein.
Inspired by the mentioned discussion previously works in the HIV infection model with CD4 + Tcells, we consider the model [13] under the generalized Caputo fractional derivative [31]. The main aim of this research is to examine the existence and uniqueness of the solutions for the fractional HIV model via the famous fixed point theorems such as Banach's and Larey-Schauder's nonlinear alternative. Moreover, stability analysis in the context of various Ulam's stability is discussed. Finally, we use the predictor-corrector algorithm proposed by Erturk et al. [30] to find the approximate solutions of the fractional HIV model for various fractional derivative orders. This paper is organized as follows: Some preliminary facts, viral definitions, theoretical results, and descriptions of the fractional HIV model are presented in Section 2. In Section 3, we derive the relation between the Volterra integral equation and the proposed model. Also, with the help of these models, the existence and uniqueness of solutions are discussed by using the concepts of Banachs fixed point theorem and Larey-Schuader's nonlinear alternative. In Section 4, Ulam's stability of the proposed model in the frame of Ulam-Hyers, generalized Ulam-Hyers, Ulam-Hyers-Rassias, and generalized Ulam-Hyers-Rassias is provided. As an application, the general numerical scheme and numerical simulations illustrate the theoretical results in Section 5. Finally, the details of the conclusion have been completed in Section 6.

Background materials
This section gives some notations, basic definitions, and fundamental results of generalized fractional derivative and integral operators which are useful throughout this paper, see more details in [31,46,47]. 46]). The generalized fractional integral of the function f , ρ I α a + , of order α > 0 where ρ > 0, is expressed (provided it exists) as: and Γ(·) represents the Gamma function. 47]). The generalized Riemann-type fractional derivative of the function f , ρ D α a + , of order α > 0 where ρ > 0, is defined by 31]). The new generalized Caputo-type fractional derivative of a function f , Cρ D α a + , of order α > 0 where ρ > 0, is defined by the following relation holds: In particular, if 0 < α ≤ 1, we have Lemma 2.6. (Banach fixed point theorem [48]). Let X be a Banach space, and D be a non-empty and closed subset of X. If the operator Q : D → D is a contraction, then, Q has a fixed point in D.
Lemma 2.7. (Leray-Schauder nonlinear alternative [49]). Let X be a Banach space, C a closed, convex subset of X, U an open subset of C and 0 ∈ U. Suppose that D : U → C is a continuous, compact (that is, D(U) is a relatively compact subset of C) map. Then either (i) D has a fixed point in U; or (ii) there is a x ∈ ∂U (the boundary of U in C) and ν ∈ (0, 1) with x = νD(x).

Model description
This work is based on the HIV models proposed by [13,24], which are considered the antiretroviral therapy of reverse transcriptase (RT) inhibitor. The unknown variables could be considered in the model represent the following: S (t) is the density of susceptible CD4 + T-cells; I(t) is the density of infected CD4 + T-cells before reverse transcription (pre-RT class); V(t) is the density of infected CD4 + T-cells in which reverse transcription is completed (post-RT class) and so are capable of producing virus; L(t) is the virus density. Besides, the following positive parameters have been considered: λ is inflows rate of CD4 + T-cells; k is an interaction infection rate of CD4 + T-cells; µ 1 is the natural death rate of CD4 + T-cells; η is the efficacy of RT inhibitor (η ∈ (0, 1)); σ is transition rate of pre-RT class infected CD4 + T-cells to pose-RT class; b is reverting rate of infected cells to uninfected class due to non-completion of reverse transcription; µ 2 is the death rate of infected CD4 + T-cells; ω is the death rate of actively infected CD4 + T-cells; N is the total number of viral particles produced by an infected CD4 + T-cells; c is the clearance rate of the virus.
We study the HIV model infection with CD4 + T-cells (1.3) under the generalized Caputo fractional derivative [31] and the meaning of all parameters is shown in the previous section: where Cρ D α 0 + (·) is the generalized Caputo-type fractional derivative of order α with 0 < α ≤ 1 and ρ > 0. We simplify the model (2.7) for easy description in the following setting: where the nonlinear functions g 1 -g 4 are given by with the initial conditions (S (0),

Existence results
This section discusses the existence and uniqueness of solutions of the proposed model by helping fixed point theorems. Let which is equivalent to the Volterra integral equation where ξ(t) = (S (t), I(t), V(t), L(t)) T with the initial conditions ξ(0) = (S 0 , I 0 , V 0 , L 0 ) T , and G(t, ξ(t)) = (g i (t, S , I, V, L)) T for i = 1, 2, 3, 4.
In view of the integral equation (3.2), an operator Q : E → E defined by It should be noticed that the problem (3.1) which is equivalent to the model (2.7) has solutions if and only if the operator Q has fixed points. (H 1 ) there exists a constant L G > 0 such that |G(t, ξ 1 (t)) − G(t, ξ 2 (t))| ≤ L G |ξ 1 (t) − ξ 2 (t)|, for any ξ 1 , ξ 2 ∈ E and for all t ∈ [0, T ].

If
L G T ρα ρ α Γ(α + 1) < 1, Proof. Early the proving, we convert the problem (3.1) into a fixed point problem, ξ = Qξ, where Q is given as in (3.3). Applying the Banach contraction principle, we show that Q has a unique fixed point, which is the unique solution of the problem (3.1). . (3.5) Note that B r 1 is a bounded, closed, and convex subset of E.
Step I. We show that QB r 1 ⊂ B r 1 . For any ξ ∈ B r 1 , we have which implies that QB r 1 ⊂ B r 1 .

Existence result via Leray-Schauder's nonlinear alternative
The Leray-Schauder's nonlinear alternative (Lemma 2.7) is presented for our first existence result.
Secondly, we show that the operator Q maps bounded sets into equicontinuous sets of E. Let τ 1 , τ 2 ∈ [0, T ] with τ 1 < τ 2 and ξ ∈ B r 2 . Then, we obtain Clearly, which independent of ξ ∈ B r 2 , the right hand side of the inequality (3.7) tends to zero as τ 2 → τ 1 . Therefore, by the Arzelá-Ascoli theorem, Q : E → E is completely continuous.
Finally, we show that the boundedness of the set of all solutions to ξ = κQξ for κ ∈ (0, 1). Let ξ ∈ E be a solution. Then, for t ∈ [0, T ], and following computation similar process to the first step, we get Taking the norm into the both sides of the above inequality, for t ∈ [0, T ], which implies that In view of (H 3 ), there exists a constant M 2 > 0 such that ξ M 2 . Let us set K := {ξ ∈ E : ξ < M 2 }. Note that the operator Q : K → C is continuous and completely continuous. From the choice of K, there is no ξ ∈ ∂K such that ξ = κQξ for some κ ∈ (0, 1). Hence, by the nonlinear alternative of Leray-Schauder type (Lemma 2.7), we deduce that the operator Q has a fixed point ξ ∈ K which implies that the model (2.7) has at least one solution on [0, T ]. The result is desired.

Ulam-Hyers stability
In this section, we are developing some sufficient conditions under which the model (2.7) will satisfy the assumptions of the different types of Ulam's stability such as Ulam-Hyers stable, generalized Ulam-Hyers stable, Ulam-Hyers-Rassias stable, and generalized Ulam-Hyers-Rassias stable.
Before we state the Ulam's stability theorem, the following definitions are needed. Let > 0 be a positive real number and Φ G : [0, T ] → R + be a continuous function. We consider the following inequalities: where = max( j ) T for j = 1, 2, 3, 4.
[50] The problem (3.1) is said to be generalized Ulam-Hyers stable if there exists a function Φ G ∈ C(R + , R + ), with Φ G (0) = 0 such that, for each solution z ∈ E of inequality (4.2), there exists a solution ξ ∈ E of the problem (3.1) such that

5)
where Φ G = max(Φ G j ) T for j = 1, 2, 3, 4.  Remark 4.6. A function z ∈ E is a solution of the inequality (4.1) if and only if there exists a function w ∈ E (which depends on z) such that the following properties:

The Ulam-Hyers stability and generalized Ulam-Hyers stability
We demonstrate an important role that will be used in the proofs of Ulam-Hyers stability and generalized Ulam-Hyers stability.
Lemma 4.8. Let α ∈ (0, 1] and ρ > 0. If z ∈ E is a solution of the inequality (4.1), then z is a solution of the following inequality .
Proof. Let z be a solution of the inequality (4.1). In view of Remark 4.6 (ii), we obtain Thus, the solution of the problem (4.9) can be written By applying Remark 4.6 (i), Therefore, the inequality (4.8) is achieved.
Next, we investigate the Ulam-Hyers stability and generalized Ulam-Hyers stability results. Proof. Suppose that > 0 and let z ∈ E be any solution of the inequality (4.1). Let ξ ∈ E be the unique solution of the problem (3.1), By using (3.2) and Lemma 4.8 with the property of triangle inequality, we get This implies that |z(t) − ξ(t)| ≤ C G , where . Hence, the model (2.7) is Ulam-Hyers stable. Now, by setting Φ G ( ) = C G such that Φ G (0) = 0 yields that the model (2.7) is generalized Ulam-Hyers stable. The proof is completed.

The Ulam-Hyers-Rassias stability and generalized Ulam-Hyers-Rassias stability
We assume that the following assumption: (H 4 ) There exists an increasing function Φ G ∈ E and there exists λ Φ G > 0, such that for any t ∈ [0, T ], the following integral inequality (4.10) Let us present an important lemma that will be used in the proofs of the Ulam-Hyers-Rassias stability and generalized Ulam-Hyers-Rassias stability results. Lemma 4.10. Let α ∈ (0, 1] and ρ > 0. If z ∈ E is a solution of the inequality (4.2), then z is a solution of the following inequality Proof. Let z be a solution of the inequality (4.2). In view of Remark 4.7 (ii), we obtain Thus, the solution of the problem (4.12) can be written . By applying Remark 4.7 (i), . Hence, the inequality (4.8) is achieved.
In the last result, we are ready to prove the Ulam-Hyers-Rassias stability and generalized Ulam-Hyers-Rassias stability results. Proof. Let > 0 and z ∈ E be the solution of the inequality (4.3). Let ξ ∈ E be the unique solution of the problem (3.1). By applying (3.2) and Lemma 4.10, we have , we get the following inequality Therefore, the model (2.7) is Ulam-Hyers-Rassias stable. Moreover, if we set = 1, in (4.13), with Φ G (0) = 0, then the model (2.7) is generalized Ulam-Hyers-Rassias stable. The proof is completed.

Simulation results and discussion
In this section, the fractional variant of the HIV model under consideration via the generalized Caputo fractional derivative is numerically simulated by using the predictor-corrector algorithm as proposed in [30].
The first step of the algorithm, the interval [0, T ] is divided into n unequal subintervals {[t k , t k+1 ], k = 0, 1, . . . , n − 1} using the mesh points where h = T ρ n . The approximate solutions of the HIV model (2.7) can be written as According to the biological works of literature [13,24], we will consider the model's constant parameters and initial states as shown in Table 1.  It can be easily seen that when the values of fractional orders α approach to 1, the solutions S α (t) decrease to the integer order soluion S 1 (t) = S (t) of the classical model [13], see in Figure 1, but the solutions I α (t), V α (t) and L α (t) also increase to the interger order solution I 1 (t) = I(t), V 1 (t) = V(t) and L 1 (t) = L(t), respectively, of the classical model [13] as see in

Conclusions
The analysis of the mathematical model of HIV infections with antiretroviral therapy involving the generalized Caputo fractional derivative established in this paper. We showed that the solution of the considered model (2.7) exists and has at least one solution with the help of the Larey-Schauder nonlinear alternative and the unique solution proved by using Banach's fixed point theorem. The stability of the solution was examined by employing the various Ulam's stability such as Ulam-Hyers stable, generalized Ulam-Hyers stable, Ulam-Hyers-Rassias stable, and generalized Ulam-Hyers-Rassias stable. Further on using the predictor-corrector algorithm, we illustrated the approximate solutions for the different fractional-order α via the Matlab program and we concluded that the solutions of the considered model (2.7) tend to the standard solution when the fractional order α → 1. For further work, one can be developed and applied to the different types of fractional derivatives with the similar techniques in many real-world problems. We hope that this work would be a good alternative way for various researches.