A Fractional Order Model for Viral Infection with Cure of Infected Cells and Humoral Immunity

In this paper, we study the dynamics of a viral infectionmodel formulated by five fractional differential equations (FDEs) to describe the interactions between host cells, virus, and humoral immunity presented by antibodies. The infection transmission process is modeled by Hattaf-Yousfi functional response which covers several forms of incidence rate existing in the literature.We first show that the model is mathematically and biologically well-posed. By constructing suitable Lyapunov functionals, the global stability of equilibria is established and characterized by two threshold parameters. Finally, some numerical simulations are presented to illustrate our theoretical analysis.


Introduction
The immune response plays an important role to control the dynamics of viral infections such as human immunodeficiency virus (HIV), hepatitis B virus (HBV), hepatitis C virus (HCV), and human T-cell leukemia virus (HTLV).Therefore, many mathematical models have been developed to incorporate the role of immune response in viral infections.Some of these models considered the cellular immune response mediated by cytotoxic T lymphocytes (CTL) cells that attack and kill the infected cells [1][2][3][4][5] and the others considered the humoral immune response based on the antibodies which are produced by the B-cells and are programmed to neutralize the viruses [6][7][8][9][10][11].However, all these models have been formulated by using ordinary differential equations (ODEs) in which the memory effect is neglected while the immune response involves memory [12,13].
Fractional derivative is a generalization of integer derivative and it is a suitable tool to model real phenomena with memory which exists in most biological systems [14][15][16].The fractional derivative is a nonlocal operator in contrast to integer derivative.This means that if we want to compute the fractional derivative at some point  =  1 , it is necessary to take into account the entire history from the starting point  =  0 up to the point  =  1 .For these reasons, modeling some real process by using fractional derivative has drawn attention of several authors in various fields [17][18][19][20][21][22].In biology, it has been shown that the fractional derivative is useful to analyse the rheological proprieties of cells [23].Furthermore, it has been deduced that the membranes of cells of biological organism have fractional order electrical conductance [24].Recently, much works have been done on modeling the dynamics of viral infections with FDEs [25][26][27][28][29][30][31].These works ignored the impact of the immune response and the majority of them deal only with the local stability.
In some viral infections, the humoral immune response is more effective than cellular immune response [32].For this reason, we improve the above ODE and FDE models by proposing a new fractional order model that describes the interactions between susceptible host cells, viral particles, and the humoral immune response mediated by the antibodies; that is,    () =  −  −  (, V) V + ,    () =  (, V) V − ( +  + ) , 2 International Journal of Differential Equations where (), (), (), V(), and () are the concentrations of susceptible host cells, latently infected cells (infected cells which are not yet able to produce virions), productive infected cells, free virus particles, and antibodies at time , respectively.Susceptible host cells are assumed to be produced at a constant rate , die at the rate , and become infected by virus at the rate (, V)V.Latently infected cells die at the rate  and return to the uninfected state by loss of all covalently closed circular DNA (cccDNA) from their nucleus at the rate .Productive infected cells are produced from latently infected cells at the rate  and die at the rate .Free virus particles are produced from productive infected cells at the rate , cleared at the rate V, and are neutralized by antibodies at the rate V.Antibodies are activated against virus at the rate V and die at the rate ℎ.
In system (1),   represents the Caputo fractional derivative of order  defined for an arbitrary function  by with 0 <  ≤ 1 [33].Further, the infection transmission process in (1) is modeled by Hattaf-Yousfi functional response [34] which was recently used in [35,36] and has the form (, V) = /( 0 +  1  +  2 V +  3 V), where  0 ,  1 ,  2 ,  3 ≥ 0 are the saturation factors measuring the psychological or inhibitory effect and  > 0 is the infection rate.In addition, this functional response generalizes many common types existing in the literature such as the specific functional response proposed by Hattaf et al. in [37] and used in [2,31] when  0 = 1; the Crowley-Martin functional response introduced in [38] and used in [39] when  0 = 1 and  3 =  1  2 ; and the Beddington-DeAngelis functional response proposed in [40,41] and used in [3,4,10] when  0 = 1 and  3 = 0. Also, the Hattaf-Yousfi functional response is reduced to the saturated incidence rate used in [9] when  0 = 1 and  1 =  3 = 0 and the standard incidence function used in [27] when  0 =  3 = 0 and  1 =  2 = 1, and it was simplified to the bilinear incidence rate used in [5,6] when  0 = 1 and On the other hand, system (1) becomes a model with ODEs when  = 1, which improves and generalizes the ODE model with bilinear incidence rate [42], the ODE model with saturated incidence rate [43], and the ODE model with specific functional response [44].
The rest of the paper is organized as follows.The next section deals with some basic proprieties of the solutions and the existence of equilibria.The global stability of equilibria is established in Section 3. To verify our theoretical results, we provide some numerical simulations in Section 4, and we conclude in Section 5.

Basic Properties and Equilibria
In this section, we will show that our model is well-posed and we discuss the existence of equilibria.
Since system (1) describes the evolution of cells, then we need to prove that the cell numbers should remain nonnegative and bounded.For biological considerations, we assume that the initial conditions of (1) satisfy (3) Then we have the following result.
Theorem 1. Assume that the initial conditions satisfy (3).Then there exists a unique solution of system (1) defined on [0, +∞).Moreover, this solution remains nonnegative and bounded for all  ≥ 0.
Proof.First, system (1) can be written as follows: where ) and  () = ( ( ( ) . ( It is important to note that when  = 1, (4) becomes a system with ODEs.In this case, we refer the reader to [45] for the existence of solutions and to the works [46][47][48][49][50] for the stability of equilibria.In the case of FDEs, we will use Lemma 2.4 in [31] to prove the existence and uniqueness of solutions.Hence, we put We discuss four cases: (i) If  0 ̸ = 0, () can be formulated as follows: where Hence, (ii) If  1 ̸ = 0, we can write () in the form where Moreover, we get where Further, we obtain (iv) If  3 ̸ = 0, we have where Then Hence, the conditions of Lemma 2.4 in [31] are verified.Then system (1) has a unique solution on [0, +∞).Now, we show the nonnegativity of solutions.By (1), we have As in [31, Theorem 2.7], we deduce that the solution of ( 1) is nonnegative.
Finally, we prove the boundedness of solutions.We define the function Then, we have where  = min{, , /2, , ℎ}.Thus, we obtain Since 0 ≤   (−  ) ≤ 1, we get This completes the proof.Now, we discuss the existence of equilibria.It is clear that system (1) has always an infection-free equilibrium  0 (/, 0, 0, 0, 0).Then the basic reproduction number of ( 1) is as follows: To find the other equilibria, we solve the following system: From ( 29), we get  = 0 or V = ℎ/.Then we discuss two cases.
Let us introduce the reproduction number for humoral immunity as follows: which 1/ℎ denotes the average life expectancy of antibodies and V 1 is the number of free viruses at  1 .For the biological significance,  1 represents the average number of the antibodies activated by virus.
In this case, (32)  We summarize the above discussions in the following theorem.
(ii) If  0 > 1, then system (1) has an infection equilibrium without humoral immunity of the form

Global Stability of Equilibria
In this section, we focus on the global stability of equilibria.
Theorem 3. If  0 ≤ 1, then the infection-free equilibrium  0 is globally asymptotically stable and it becomes unstable if  0 > 1.
Proof.The proof of the first part of this theorem is based on the construction of a suitable Lyapunov functional that satisfies the conditions given in [51, Lemma 4.6].Hence, we define a Lyapunov functional as follows: where Φ() =  − 1 − ln() for  > 0. It is not hard to show that the functional  0 is nonnegative.In fact, the function Φ has a global minimum at  = 1.Consequently, Φ() ≥ 0 for all  > 0.
The proof of the instability of  0 is based on the computation of the Jacobean matrix of system (1) and the results presented in [53][54][55].The Jacobean matrix of (1) at any equilibrium (, , , V, ) is given by

Theorem 4.
(i) The infection equilibrium without humoral immunity  1 is globally asymptotically stable if  0 > 1,  1 ≤ 1, and Proof.Define a Lyapunov functional as follows: Calculating the fractional derivative of  1 (), we get Using , we obtain Hence, Using the arithmetic-geometric inequality, we have Since It is easy to see that this condition is equivalent to (44).Furthermore,
At  1 , the characteristic equation of ( 40) is given as follows: where We can easily see that (50) has the root  1 = V 1 − ℎ.Then, when  1 > 1, we have  1 > 0. In this case,  1 is unstable.
It is important to note that when  is sufficiently small or  is sufficiently large, the two conditions ( 44) and ( 52) are satisfied.Then, we have the following corollary.Corollary 6. Assume that  0 > 1.When  is sufficiently small or  is sufficiently large, then we have the following: (i) The infection equilibrium without humoral immunity  1 is globally asymptotically stable if  1 ≤ 1.
(ii) The infection equilibrium with humoral immunity  2 is globally asymptotically stable if  1 > 1.

Numerical Simulations
In this section, we validate our theoretical results to HIV infection.Firstly, we take the parameter values as shown in Table 1.
Next, we take  = 0.0004 and do not change the other parameter values.In this case, we have  1 = 3.3338, ℎ = 0.0000024, and ( +  + )( 0  +  2 ℎ) + ( 1  +  3 ℎ) = 0.000006.Hence, condition ( 52) is satisfied.Consequently, system (1) has an infection equilibrium with humoral immunity  2 (423.4261,92.0442, 3.4090, 500, 245.4473) which is globally asymptotically stable.Figure 3 illustrates this result.We can observe that the activation of the humoral immune response increases the healthy cells and decreases the productive infected cells and viral load to a lower levels but it is not able to eradicate the infection.

Conclusion
In the present paper, we have studied the dynamics of a viral infection model by taking into account the memory effect represented by the Caputo fractional derivative and the humoral immunity.We have proved that the solutions of the model are nonnegative and bounded which assure the well-posedness.We have shown that the proposed model has three infection equilibriums, namely, the infection-free equilibrium  0 , the infection equilibrium without humoral immunity  1 , and the infection equilibrium with humoral immunity  2 .By constructing suitable Lyapunov functionals, the global stability of these equilibria is fully determined by two threshold parameters  0 and  1 .More precisely, when  0 ≤ 1,  0 is globally asymptotically stable, whereas if  0 > 1, it becomes unstable and another equilibrium point appears, that is,  1 , which is globally asymptotically stable whenever  1 ≤ 1 and condition (44) is satisfied.In the case that  1 > 1,  1 becomes unstable and there exists another equilibrium point  2 which is globally asymptotically stable when condition (52) is satisfied.In addition, we remarked that when  is sufficiently small or  is sufficiently large, conditions (44) and ( 52) are verified, and then the global stability of  1 and  2 is characterized only by  0 and  1 .
From our theoretical and numerical results, we deduce that the order of the fractional derivative  has no effect on the dynamics of the model.However, when the value of  decreases (long memory), the solutions of our model converge rapidly to the steady states (see Figures 1-3).This behavior can be explained by the memory term 1/Γ(1 − )( − )  included in the fractional derivative which represents