# Stability Analysis for a Fractional HIV Infection Model with Nonlinear Incidence

Authors: Linli Zhang; Gang Huang; Anping Liu; Ruili Fan
Journal: Discrete Dynamics in Nature and Society
 (2015)

---

## Abstract

We introduce the fractional-order derivatives into an
HIV infection model with nonlinear incidence and show that the established
model in this paper possesses nonnegative solution, as desired in any population
dynamics. We also deal with the stability of the infection-free equilibrium,
the immune-absence equilibrium, and the immune-presence equilibrium. Numerical
simulations are carried out to illustrate the results.


---

## Body

## 1. Introduction

### 1.1. The History of Fractional Calculus

Fractional calculus which is a branch of mathematical analysis extends derivatives and integrals to an arbitrary order (real, even, or complex order). The study on fractional calculus started at the end of seventeenth century and the first reference was proposed by Leibniz and L’Hospital in 1695, in which half-order derivative was mentioned. The emergence of fractional calculus aroused interest of some famous mathematicians such as Lacroix, Lagrange, Fourier, and Laplace. J. L. Lagrange developed the law of exponents for differential operators in 1772, which promoted indirectly the development of fractional calculus. P. S. Laplace defined a fractional derivative by means of integral in 1812 and S. F. Lacroix and J. B. J. Fourier in succession mentioned the arbitrary order derivative, respectively, in 1819 and in 1822. However, the basic theory of fractional calculus was established with the studies of Liouville, Grunwald, Letnikov, and Riemann until the end of the nineteenth century.At the initial stage, the development of fractional calculus had been restricted to the pure mathematics research. Until the nineteen nineties, the theory of fractional calculus was applied to nature and social sciences. In the field of anomalous diffusion, the researchers applied fractional calculus to describe the diffusion processes which do not conform to the brown motion and proposed the fractional anomalous diffusion model. According to the theoretical analysis and the experimental data or results, they found that the fractional model is more reasonable to describe these processes [1]. Afterwards, many mathematicians and applied researchers also have tried to demonstrate applications of fractional differentials in the areas of non-Newtonian fluids [2], signal processing [3–5], viscoelasticity [6, 7], fluid-dynamic traffic model [8], colored noise [9], bioengineering [10], solid mechanics [11], continuum and statistical mechanics [12], and economics [13] and brought new research view for those fields.


### 1.2. Mathematical Modeling

In biology, it has been deduced that the membranes of cells of biological organism have fractional-order electrical conductance [14]; hence some mathematical models which describe cells behavior are classified into groups of non-integer-order models. In the field of rheology, fractional-order derivatives embody essential features of cell rheological behavior and have enjoyed greatest success [15]. Some researchers also found that fractional ordinary differential equations are naturally related to systems with memory which exists in most biological systems. And fractional-order equations are also closely related to fractals, which are abundant in biological systems.Mathematical models have been proven valuable in understanding the dynamics of viral infection. Although a large number of works on modeling the dynamics of viral infection have been done [16–20], it has been restricted to integer-order (delay) differential equations. In recent years, it has turned out that many phenomena in virus infection can be described very successfully by the models using fractional-order differential equations [21, 22]. Motivated by these references, in this paper, we will consider a fractional-order HIV model.Among a large number of virus infection models described with integer-order or delay differential equations, a typical model is given as the following simple three-dimensional system:(1)x′t=λ-dxt-βxtvt,y′t=βxtvt-pyt,v′t=kyt-uvt.Here x(t) represents the concentration of uninfected cells at time t; y(t) represents the concentration of infected cells that produce virus at time t; v(t) represents the concentration of viruses at time t. d(d>0) and β(β>0) are the death rate of uninfected cells and the rate constant characterizing infection of the cells, respectively. p(p>0) is the death rate of infected cells due to either the virus or the immune system. Free virus is produced from infected cells at the rate ky(t) and removed at rate uv(t).The infection rate in model (1) is usually assumed to be bilinear with respect to virus v(t) and uninfected target cells x(t). However, there are some evidences which show that a bilinear infection might not be an effective assumption when the number of target cells is large enough [23–25]. And Huang et al. [17] incorporated Holling type II functional response infection rate g(x(t),v(t)) into basic virus dynamics (1), which was expressed as (2)gxt,vt=βxtvt1+qxt,where q (q≥0) is a constant. Hence, system (1) can be modified into the following system:(3)x′t=λ-dxt-βxtvt1+qxt,y′t=βxtvt1+qxt-pyt,v′t=kyt-uvt.The immune response following viral infection is universal and necessary to eliminate or control the disease. Antibodies, cytokines, natural killer cells, B cells, and T cells are all essential components of a normal immune response to viral infection. Since cytotoxic T lymphocytes (CTLs) play a critical role in antiviral defense by attacking the virus-infected cells, Huang et al. [17] considered the immune response to model (3) with z′(t)=cy(t)z(t)-bz(t).System (3) with z′(t)=cy(t)z(t)-bz(t) can be further simplified if we take into consideration the fact that an average life span of viral particles is usually significantly shorter than one of the infected cells. It can be assumed, therefore, that, compared with a “slow” variation of the infected cells level, the virus load v(t) relatively quickly reaches a quasiequilibrium level. The equality v˙(t)=0 holds in the quasiequilibrium state and hence v(t)=ky(t)/u. This assumption is referred to as “separation of time scales” and is in common use in the virus dynamics [18]. We have to stress that this assumption does not imply that the virus concentration v(t) remains constant; on the contrary, it is assumed to be proportional to the varying concentration of infected cells y(t). Accordingly, system (3) with z′(t)=cy(t)z(t)-bz(t) can now be reformulated as a system of three differential equations:(4)x′(t)=λ-dx(t)-βx(t)y(t)1+qx(t),y′(t)=βx(t)y(t)1+qx(t)-ay(t)-py(t)z(t),z′t=cytzt-bzt.Further we introduce the fractional order into model (4) and a new system can be described by the following set of FODEs of order α>0:(5)xα(t)=λ-dx-βxy1+qx,yα(t)=βxy1+qx-ay-pyz,zα(t)=cyz-bz,(6)x(0)=x0,y(0)=y0,z(0)=z0,where the immune response is assumed to get stronger at a rate cyz, which is proportional to the number of infected cells and their current concentration; the immune response decays exponentially at a rate bz, which is proportional to their current concentration; the parameter q expresses the efficacy of nonlytic component.The paper is organized as follows. In the next section, we give two definitions and two lemmas about fractional calculus. In Section3, we show the nonnegativeness of the solutions of system (5) with initial condition (6). In Section 4, we give a detailed stability analysis for three equilibria. Finally numerical simulations are presented to illustrate the obtained results in Section 5.


## 1.1. The History of Fractional Calculus

Fractional calculus which is a branch of mathematical analysis extends derivatives and integrals to an arbitrary order (real, even, or complex order). The study on fractional calculus started at the end of seventeenth century and the first reference was proposed by Leibniz and L’Hospital in 1695, in which half-order derivative was mentioned. The emergence of fractional calculus aroused interest of some famous mathematicians such as Lacroix, Lagrange, Fourier, and Laplace. J. L. Lagrange developed the law of exponents for differential operators in 1772, which promoted indirectly the development of fractional calculus. P. S. Laplace defined a fractional derivative by means of integral in 1812 and S. F. Lacroix and J. B. J. Fourier in succession mentioned the arbitrary order derivative, respectively, in 1819 and in 1822. However, the basic theory of fractional calculus was established with the studies of Liouville, Grunwald, Letnikov, and Riemann until the end of the nineteenth century.At the initial stage, the development of fractional calculus had been restricted to the pure mathematics research. Until the nineteen nineties, the theory of fractional calculus was applied to nature and social sciences. In the field of anomalous diffusion, the researchers applied fractional calculus to describe the diffusion processes which do not conform to the brown motion and proposed the fractional anomalous diffusion model. According to the theoretical analysis and the experimental data or results, they found that the fractional model is more reasonable to describe these processes [1]. Afterwards, many mathematicians and applied researchers also have tried to demonstrate applications of fractional differentials in the areas of non-Newtonian fluids [2], signal processing [3–5], viscoelasticity [6, 7], fluid-dynamic traffic model [8], colored noise [9], bioengineering [10], solid mechanics [11], continuum and statistical mechanics [12], and economics [13] and brought new research view for those fields.


## 1.2. Mathematical Modeling

In biology, it has been deduced that the membranes of cells of biological organism have fractional-order electrical conductance [14]; hence some mathematical models which describe cells behavior are classified into groups of non-integer-order models. In the field of rheology, fractional-order derivatives embody essential features of cell rheological behavior and have enjoyed greatest success [15]. Some researchers also found that fractional ordinary differential equations are naturally related to systems with memory which exists in most biological systems. And fractional-order equations are also closely related to fractals, which are abundant in biological systems.Mathematical models have been proven valuable in understanding the dynamics of viral infection. Although a large number of works on modeling the dynamics of viral infection have been done [16–20], it has been restricted to integer-order (delay) differential equations. In recent years, it has turned out that many phenomena in virus infection can be described very successfully by the models using fractional-order differential equations [21, 22]. Motivated by these references, in this paper, we will consider a fractional-order HIV model.Among a large number of virus infection models described with integer-order or delay differential equations, a typical model is given as the following simple three-dimensional system:(1)x′t=λ-dxt-βxtvt,y′t=βxtvt-pyt,v′t=kyt-uvt.Here x(t) represents the concentration of uninfected cells at time t; y(t) represents the concentration of infected cells that produce virus at time t; v(t) represents the concentration of viruses at time t. d(d>0) and β(β>0) are the death rate of uninfected cells and the rate constant characterizing infection of the cells, respectively. p(p>0) is the death rate of infected cells due to either the virus or the immune system. Free virus is produced from infected cells at the rate ky(t) and removed at rate uv(t).The infection rate in model (1) is usually assumed to be bilinear with respect to virus v(t) and uninfected target cells x(t). However, there are some evidences which show that a bilinear infection might not be an effective assumption when the number of target cells is large enough [23–25]. And Huang et al. [17] incorporated Holling type II functional response infection rate g(x(t),v(t)) into basic virus dynamics (1), which was expressed as (2)gxt,vt=βxtvt1+qxt,where q (q≥0) is a constant. Hence, system (1) can be modified into the following system:(3)x′t=λ-dxt-βxtvt1+qxt,y′t=βxtvt1+qxt-pyt,v′t=kyt-uvt.The immune response following viral infection is universal and necessary to eliminate or control the disease. Antibodies, cytokines, natural killer cells, B cells, and T cells are all essential components of a normal immune response to viral infection. Since cytotoxic T lymphocytes (CTLs) play a critical role in antiviral defense by attacking the virus-infected cells, Huang et al. [17] considered the immune response to model (3) with z′(t)=cy(t)z(t)-bz(t).System (3) with z′(t)=cy(t)z(t)-bz(t) can be further simplified if we take into consideration the fact that an average life span of viral particles is usually significantly shorter than one of the infected cells. It can be assumed, therefore, that, compared with a “slow” variation of the infected cells level, the virus load v(t) relatively quickly reaches a quasiequilibrium level. The equality v˙(t)=0 holds in the quasiequilibrium state and hence v(t)=ky(t)/u. This assumption is referred to as “separation of time scales” and is in common use in the virus dynamics [18]. We have to stress that this assumption does not imply that the virus concentration v(t) remains constant; on the contrary, it is assumed to be proportional to the varying concentration of infected cells y(t). Accordingly, system (3) with z′(t)=cy(t)z(t)-bz(t) can now be reformulated as a system of three differential equations:(4)x′(t)=λ-dx(t)-βx(t)y(t)1+qx(t),y′(t)=βx(t)y(t)1+qx(t)-ay(t)-py(t)z(t),z′t=cytzt-bzt.Further we introduce the fractional order into model (4) and a new system can be described by the following set of FODEs of order α>0:(5)xα(t)=λ-dx-βxy1+qx,yα(t)=βxy1+qx-ay-pyz,zα(t)=cyz-bz,(6)x(0)=x0,y(0)=y0,z(0)=z0,where the immune response is assumed to get stronger at a rate cyz, which is proportional to the number of infected cells and their current concentration; the immune response decays exponentially at a rate bz, which is proportional to their current concentration; the parameter q expresses the efficacy of nonlytic component.The paper is organized as follows. In the next section, we give two definitions and two lemmas about fractional calculus. In Section3, we show the nonnegativeness of the solutions of system (5) with initial condition (6). In Section 4, we give a detailed stability analysis for three equilibria. Finally numerical simulations are presented to illustrate the obtained results in Section 5.


## 2. Fractional Calculus

Definition 1 (see [26]).
The Riemann-Liouville (R-L) fractional integral operator of orderα>0 of a function f:R+→R is defined as (7)Iαf(x)=1Γ(α)∫0x‍x-tα-1f(t)dt.Here Γ(·) is the Euler gamma function which is defined as (8)Γ(n)=∫0∞‍tn-1e-tdt.This function is generalization of a factorial in the following form: (9)Γ(n)=(n-1)!.Definition 2 (see [26]).
The Caputo (C) fractional derivative of orderα>0, n-1<α<n, n∈N, is defined as (10)Dαf(t)=In-αDnf(t)=1Γ(n-α)∫0t‍f(n)(s)t-sα+1-nds,where the function f(t) has absolutely continuous derivatives up to order (n-1). In particular, when 0<α<1, one has (11)Dαf(t)=1Γ(1-α)∫0t‍f′(s)t-sαds.In this paper we use Caputo fractional derivative definition. The main advantage of Caputo’s definition is that the initial conditions for fractional differential equations with Caputo derivatives take the same form as that for integer-order differential equations.Lemma 3 (see [26]).
Consider the following commensurate fractional-order system:(12)Dαx=f(x),x(0)=x0,with 0<α≤1 and x∈Rn. The equilibrium points of system (12) are calculated by solving the following equation: f(x)=0. These points are locally asymptotically stable if all eigenvalues ri of Jacobian matrix J=∂f/∂x evaluated at the equilibrium points satisfy(13)argri>θπ2.Definition 4 (see [27]).
The discriminantD(f) of a polynomial (14)fx=xn+c1xn-1+c2xn-2+⋯+cnis defined by D(f)=(-1)nn-1/2R(f,f′), where f′ is the derivative of f. If g(x)=xl+d1xl-1+d2xl-2+⋯+dn, R(f,g) is the determinant of the corresponding Sylvester (n+l)⊗(n+l) matrix. The Sylvester matrix is formed by filling the matrix beginning with the upper left corner with the coefficients of f(x) and then shifting down one row and one column to the right and filling in the coefficients starting there until they hit the right side. The process is then repeated for the coefficients of g(x).Lemma 5 (see [27]).
For the polynomial equation,(15)P(λ)=λn+h1λn-1+h2λn-2+⋯+hn=0,the conditions which make all the roots of (15) satisfy (13) are displayed as follows: (i)
forn=1, the condition is h1>0;(ii)
forn=2, the conditions are either Routh-Hurwitz conditions or (16)h1<0,4h2>h12,tan-14h2-h12h1>απ2;(iii)
forn=3,
(a)
if the discriminant ofP(λ), D(P) is positive, then Routh-Hurwitz conditions are the necessary and sufficient conditions; that is, h1>0, h3>0, and h1h2>h3 if D(P)>0;(b)
ifD(P)<0, h1≥0, h2≥0, and h3>0, then (13) for (15) holds when α<2/3;(c)
ifD(P)<0, h1<0, and h2<0, then (13) for (15) holds when α>2/3;(d)
ifD(P)<0, h1>0, h2>0, and h1h2=h3, then (13) for (15) holds for all α∈[0,1).


## 3. Nonnegative Solutions

LetR+3={W∈R3:W≥0} and W(t)=(x(t),y(t),z(t))T. For the proof of the theorem about nonnegative solutions, we would need the following lemma.Lemma 6 (see [28] generalized mean value theorem).
Letf(x)∈C[a,b] and Dαf(x)∈C[a,b] for 0<α≤1. Then one has (17)fx=fa+1ΓαDαfξx-aα,with 0≤ξ≤x, ∀x∈(a,b], where Γ(x)=∫0+∞‍tx-1e-tdt.Remark 7.
Suppose thatf(x)∈C[a,b] and Dαf(x)∈C[a,b], for 0<α≤1. It is clear from Lemma 6 that if Dαf(x)≥0, ∀x∈(a,b), then f(x) is nondecreasing for each x∈[a,b]. If Dαf(x)≤0, ∀x∈(a,b), then f(x) is nonincreasing for each x∈[a,b].We now prove the main theorem.Theorem 8.
There is a unique solution for the initial value problem (5) with (6) and the solution remains in R+3.Proof.
From Theorem 3.1 and Remark 3.2 of [29], we obtain the solution on (0,∞) solving the initial value problem (5) with (6) which is not only existent but also unique. Next, we will show the nonnegative orthant R+3 is a positively invariant region. What is needed for this is to show that, on each hyperplane bounding the nonnegative orthant, the vector field points into R+3. From (5), we find(18)xα(t)x=0=λ≥0,yα(t)y=0=0,zα(t)z=0=0.By Remark 7, the solution will remain in R+3.


## 4. Equilibrium States and Their Stability

In this section, we investigate the stability of the fractional-order model of HIV infection of CD4 + T cells, that is, system (5) with (6). Consider the initial value problem (5) with (6) with α satisfying 0<α≤1.In order to obtain the equilibria of system (5), we set xα(t)=0, yαt=0, and zαt=0 and we have(19)λ-dx-βxy1+qx=0,βxy1+qx-ay-pyz=0,cyz-bz=0.System (5) has three types of relevant nonnegative equilibrium states. System (5) always has an infection-free equilibrium E0=(x0,0,0), where (20)x0=λd,which means that the infected cells are cleared. The basic reproductive number of the viruses for system (5) is given by(21)R0=λβa(d+qλ).This number describes the average number of newly generated infected cells from one infected cell at the beginning of the infection process. When R0>1, in addition to the infection-free equilibrium E0, there is an immune-absence equilibrium E1=(x1,y1,z1), where (22)x1=aβ-aq,y1=d+qλβ-aq(R0-1),z1=0.Note that R0>1 implies β>aq. When the HIV infection is in the immune-absence equilibrium E1, the infected cells and virus exist but the immune response is not activated yet. Further, we denote(23)R1=R0-bβ-aqcd+qλ.Note that R1=((β-aq)/c(d+qλ))(cy1-b)+1, which implies that R1>1 is equivalent to cy1/b>1. The latter cy1/b is seen as the immune reproductive number, which expresses the average number of activated CTLs generated from one CTL during its life time 1/b through the stimulation of the infected cells y1. It is reasonable that immune response is activated in the case where R1>1. When R1>1, in addition to the infection-free equilibrium E0 and the immune-absence equilibrium E1, there is an interior immune-presence equilibrium E∗=(x∗,y∗,z∗), where (24)x∗=-(bβ+cd-cqλ)+bβ+cd-cqλ2+4c2dqλ2cdq,y∗=bc,z∗=cbp(λ-dx∗)-ap.E∗ expresses the state where CTLs immune response is present.Next, we establish the local asymptotic stability of model (5) by the characteristic equation.Theorem 9.
Consider system (5). (1)
The infection-free equilibriumE0 is locally asymptotically stable if R0<1.(2)
IfR0>1, the equilibrium E0 is unstable, and if R0=1, it is a critical case.Proof.
The characteristic equation for the infection-free equilibriumE0 is given as follows: (25)det⁡-d-r-βx01+qx000βx01+qx0-a-r000-b-r=0.
It is reduced to(26)r+dr+br+a-βx01+qx0=0.
It is clear that (26) has the characteristic roots r1=-d<0 which means argr1=π>α(π/2), r2=-b<0 which means argr2=π>α(π/2), and r3=(βx0/(1+qx0))-a=a(R0-1). Since the imaginary part of characteristic root r3 is zero, R0<1 which means argr3=π>α(π/2) is necessary and sufficient to ensure the local asymptotic stability of E0. If R0>1, argr3=0<α(π/2); hence E0 is unstable. If R0=1, r3=0, which is a critical case.Theorem 10.
Consider system (5). (1)
The immune-free equilibriumE1 is locally asymptotically stable if R1<1.(2)
IfR1>1, the immune-free equilibrium E1 is unstable; if R1=1, it is a critical state.Proof.
The characteristic equation for the immune-free equilibriumE1 is given as follows: (27)det⁡-d-βy11+qx12-r-βx11+qx10βy11+qx12βx11+qx1-a-r-py100cy1-b-r=0.Using βx1/(1+qx1)=a, it is reduced to(28)r+b-cy1r2+d+βy11+qx12r+aβy11+qx12=0.The root of (28)  r1=cy1-b is negative and argr1=π>α(π/2) when R1<1, positive and argr1=0<α(π/2) when R1>1, and zero when R1=1, which is a critical case.
Now, we consider the equation(29)r2+d+βy11+qx12r+aβy11+qx12=0.Denote(30)B=βy11+qx12,C=aβy11+qx12.Since B>0 and C>0, (29) has two negative real roots and we denote them by r2 and r3. It is easy to see argr2=π>α(π/2) and |argr3|=π>α(π/2). Hence, when R1<1, the immune-free equilibrium E1 is locally asymptotically stable, when R1>1, E1 is unstable, and when R1=1, it is a critical case.To discuss the local stability of the immune-present equilibrium stateE∗, we consider the linearized system of (5) at E∗. The Jacobian matrix at E∗ is given by (31)J(E∗)=-d-βy∗1+qx∗2-βx∗1+qx∗0βy∗1+qx∗20-py∗0cz∗0.The characteristic equation of the linearized system is(32)P(r)=r3+a1r2+a2r+a3=0,where(33)a1=d+βy∗1+qx∗2>0,a2=cpy∗z∗+β2x∗y∗1+qx∗3>0,a3=cpy∗z∗d+βy∗1+qx∗2>0,a1a2-a3=β2x∗y∗1+qx∗3d+βy∗1+qx∗2>0.Based on Definition 4, we obtain the discriminant of (32)(34)DP=-1a1a2a3001a1a2a332a1a200032a1a200032a1a2=18a1a2a3+a1a22-4a13a3-4a23-27a32.Using the result (iii) of Lemmas 5 and 3, we have the following theorem.Theorem 11.
Consider system (5). Under the condition of R1>1, (1)
if the discriminant ofP(r), D(P) is positive, namely, D(P)>0, then the immune-present equilibrium E∗ is locally asymptotically stable for 0<α≤1.(2)
ifD(P)<0, then the immune-present equilibrium E∗ is locally asymptotically stable for 0<α<2/3.


## 5. Numerical Method

Atanackovic and Stankovic introduced a numerical method to solve the single linear FDE in 2004 [1]. A few years later, they developed again a method to solve the nonlinear FDE [30]. It was shown that the fractional derivative of a function f(t) with order α satisfying 0<α<1 may be expressed as(35)Dαft=1Γ2-α×f1ttα-11+∑m=1∞‍Γm-1+αΓα-1m!-α-1tαft+∑m=2∞‍Γ(m-1+α)Γ(α-1)(m-1)!∑m=2∞‍Γm-1+αΓα-1m-1!fttα+Vmfttm-1+αfttα+Vmfttm-1+α×∑m=2∞‍Γm-1+αΓα-1m-1!fttα+Vmfttm-1+αf(t)tα+Vm(f)(t)tm-1+α,where(36)Vmft=-m-1∫0t‍τm-2fτdτ,m=2,3,…,with the following properties:(37)ddtVm(f)=-(m-1)tm-2f(t),m=2,3,….We approximate Dαf(t) by using M terms in sums appearing in (35) as follows:(38)Dαft≃1Γ2-α×f1ttα-11+∑m=1M‍Γm-1+αΓα-1m!-α-1tαft+∑m=2M‍Γm-1+αΓα-1m-1!fttα+Vmfttm-1+α×fttα+Vmfttm-1+α.We can rewrite (38) as follows:(39)Dαft≃Ωθ,t,Mf1t+Φθ,t,Mft+∑m=2M‍A(θ,t,m)Vm(f)(t)tm-1+α,where (40)Ω(α,t,M)=1+∑m=1M(Γ(m-1+α)/Γ(α-1)m!)Γ(2-α)tα-1,R(a,t)=1-αtαΓ2-α,A(α,t,m)=-Γ(m-1+α)Γ(2-α)Γ(α-1)(m-1)!,Φ(α,t,M)=R(a,t)+∑m=2M‍A(α,t,m)tα.We set (41)Θ1(t)=x(t),Θm(t)=Vm(x)(t),ΘM+1(t)=y(t),ΘM+m(t)=Vm(y)(t),Θ2M+1(t)=z(t),Θ2M+m(t)=Vm(z)(t),for m=2,3,…. We can rewrite system (19) as the following form:(42)Ωα,t,MΘ1′t+Φα,t,MΘ1t+∑m=2M‍Aα,t,mΘmttm-1+α=λ-dΘ1(t)-βΘ1(t)ΘM+1(t)1+qΘ1(t),Ωα,t,MΘM+1′t+Φα,t,MΘM+1t+∑m=2M‍Aα,t,mΘM+mttm-1+α=βΘ1tΘM+1t1+qΘ1t-aΘM+1(t)-pΘM+1(t)Θ2M+1(t),Ωα,t,MΘ2M+1′t+Φα,t,MΘ2M+1t+∑m=2M‍Aα,t,mΘ2M+mttm-1+α=cΘM+1(t)Θ2M+1(t)-bΘ2M+1(t),where(43)Θmt=-m-1∫0t‍τm-2Θ1τdτ,ΘM+mt=-m-1∫0t‍τm-2ΘM+1τdτ,Θ2M+mt=-m-1∫0t‍τm-2Θ2M+1τdτ,m=2,3,…,M.Now we can rewrite (39) and (42) as the following form:(44)Θ1′t=1Ωα,t,Mλ-d+Φα,t,MΘ1t∑m=2M‍Aα,t,mΘmttm-1+α-βΘ1tΘM+1t1+qΘ1t-∑m=2M‍A(α,t,m)Θm(t)tm-1+α,Θm′(t)=-(m-1)tm-2Θ1(t),m=2,3,…,M,ΘM+1′t=1Ωα,t,MβΘ1tΘM+1t1+qΘ1t∑m=2M‍Aα,t,mΘM+mttm-1+α-a+Φα,t,MΘM+1t-pΘM+1tΘ2M+1t-∑m=2M‍A(α,t,m)ΘM+m(t)tm-1+α,ΘM+m′(t)=-(m-1)tm-2ΘM+1(t),m=2,3,…,M,Θ2M+1′t=1Ωα,t,McΘM+1tΘ2M+1t∑m=2M‍Aα,t,mΘ2M+mttm-1+α-b+Φα,t,MΘ2M+1t-∑m=2M‍A(α,t,m)Θ2M+m(t)tm-1+α,Θ2M+m′(t)=-(m-1)tm-2Θ2M+1(t),m=2,3,…,M,with the following initial conditions:(45)Θ1(δ)=x0,Θm(δ)=0,m=2,3,…,M,ΘM+1(δ)=y0,ΘM+m(δ)=0,m=2,3,…,M,Θ2M+1(δ)=z0,Θ2M+m(δ)=0,m=2,3,…,M.Now we consider the numerical solution of system of ordinary differential equations (44) with the initial conditions (45) by using the well-known Runge-Kutta method of fourth order.


## 6. Numerical Simulation and Discussion

Firstly, by using GEM (generalized Euler method) [21], we simulate system (5) with the parameter values as shown in Table 1.Table 1
Parameter
Value
Parameter
Value
λ
23.3
c
0.0031
β
0.5
d
0.09
a
0.02
p
10
b
0.15
q
0.79By direct calculation, we haveR0=31.4916, R1=30.2250, and D(P)=0.00076142>0 and the simulations display that the immune-present equilibrium E∗ is asymptotically stable for α=0.85, α=0.95, and α=1 (see Figure 1). From Figure 1, we can clearly see that, compared with the case of order α=0.85, the trajectory of model with order α=0.95 is closer to the trajectory of the model with the integer-order 1. That is, the farther from α to 1, the bigger of the trajectory difference of them.Figure 1
The diagrams show the time evolution of the trajectory of system (5) for healthy cells x(t), infected cells y(t), and the CTL z(t), respectively, with respect to α=0.85, α=0.95, and α=1, in condition of D(P)>0. The red lines correspond to x(t) in (a), y(t) in (b), and z(t) in (c) for α=0.85; the green lines correspond to x(t) in (a), y(t) in (b), and z(t) in (c) for α=0.95 and the blue lines for α=1.
(a)
Time evolution of the state variables
(b)
Time evolution of the state variables
(c)
Time evolution of the state variablesWe choose the parameter values as shown in Table2.Table 2
Parameter
Value
Parameter
Value
λ
10
c
0.051
β
0.01
d
0.18
a
0.1
p
0.6
b
0.15
q
0.01Based on the parameter values in Table2, we have R0=3.5714, R1=3.4769, and D(P)=-0.00069<0. With the same simulated method, it is shown that the immune-present equilibrium E∗ is asymptotically stable for α=0.56 and α=0.65 (see Figure 2).Figure 2
The diagrams show the time evolution of the trajectory of system (5) for healthy cells x(t), infected cells y(t), and the CTL z(t), respectively, with respect to α=0.56 and α=0.65, in condition of D(P)<0. The red lines correspond to x(t) in (a), y(t) in (b), and z(t) in (c) for α=0.56; the green lines correspond to x(t) in (a), y(t) in (b), and z(t) in (c) for α=0.65.
(a)
Time evolution of the state variables
(b)
Time evolution of the state variables
(c)
Time evolution of the state variablesNext, we use the method that is shown in the previous section to simulate system (5) by transforming system (5) to one order ordinary differential equation. Here we set the parameters values as shown in Table 3 with M=10, Δt=0.005. We have R0=9.0445, R1=8.7438, and D(P)=-10.0374<0. The approximate solutions x(t), y(t), and z(t) for α=0.6, α=0.7, α=0.8, and α=0.9 are displayed in Figure 3. It shows that the immune-present equilibrium E∗ is asymptotically stable for α=0.6<2/3 and there exists the limit circle for α>2/3.Table 3
Parameter
Value
Parameter
Value
λ
29.46
c
0.09
β
0.286
d
0.9862
a
0.0384
p
0.85
b
2.568
q
0.79Figure 3
The three-dimensional diagrams show the approximate solutions ofx(t), y(t), and z(t) for α=0.6, α=0.7, α=0.8, and α=0.9 in condition of D(P)<0.We also simulate the situation of system (5) by the method in the previous section when α=1 with the parameter values of Table 3 in Figure 4. Both Figures 1 and 4 show that the immune-present equilibrium E∗ is asymptotically stable for α=1 and that is nothing to do with D(P)>0 or D(P)<0. Hence the results of fractional-order system when α=1 are consistent with the result of integer-order HIV model (4).Figure 4
The diagrams show the approximate solutions ofx(t), y(t), and z(t) for α=1 in condition of R1>1 and D(P)<0. The panel (a) corresponds to healthy cells, (b) to infected cells, and (c) to antigen-specific CTL. The panel (d) shows the three-dimensional trajectory of system (5) with time evolution.
(a)
(b)
(c)
(d)


## 7. Conclusion

Fractional differential equations have garnered a lot of attention and appreciation due to their ability to provide an exact description of different nonlinear phenomena. The advantage of fractional-order systems is that they allow greater degrees of freedom in the model. Nowadays, more and more investigators begin to study the qualitative properties and numerical solutions of fractional-order virus infection models. In this paper, we introduced a fractional-order HIV infection model with nonlinear incidence and dealt with the mathematical behaviors of the model.We showed that system (5) possesses nonnegative solutions and studied the stability behavior of the infection-free equilibrium, the immune-absence equilibrium, and the immune-presence equilibrium. We found that the stability of the infection-free equilibrium and the immune-absence equilibrium of system (5) is the same as that of system (4). When the basic reproduction number of viruses (R0) is less than one, the infection-free equilibrium is stable; however, when R0 is more than one, the infection-free equilibrium is unstable and when the immune reproduction number (R1) is less than one, the immune-absence equilibrium is stable; however, when R1 is more than one, the immune-absence equilibrium is unstable. However, the results for the immune-presence equilibrium of system (5) are different to those of system (4). In system (4), the immune-presence equilibrium is stable when R1 is more than one, while, in system (5), when R1 is more than one, the immune-presence equilibrium is not always stable. In the condition of R1>1, when D(P)>0, the immune-presence equilibrium is stable for 0<α≤1, while when D(P)<0, the immune-presence equilibrium is stable only for 0<α<2/3. But using the simulation, we found when D(P)<0, the immune-presence equilibrium is stable for α=1. From the simulation, we also found the farther from α to 1, the bigger of their trajectory difference. These results show that the integer-order model can be viewed as a special case from the more general fractional-order model. Although a large part of results is illustrated by both theory analysis and numerical simulation, the result for the immune-presence equilibrium when D(P)<0 and α=1 can be just verified by the simulation in this paper.In this paper, we introduce the fractional calculus into the HIV infection model with nonlinear incidence and, from the theory analysis and numerical simulations, it is illustrated that the integer-order HIV infection model can be viewed as a special case of fractional-order model. We hope that this work can create interest and further do research effort in this field, since the fractional modeling might provide more insight into understanding the dynamical behaviors of such systems.



---

*Source: 563127-2015-05-14.xml*