Nonlinear Dynamics and Chaos in a Fractional-Order HIV Model

Abstract

We introduce fractional order into an HIV model. We consider the effect of viral diversity on the human immune system with frequency dependent rate of proliferation of cytotoxic T-lymphocytes (CTLs) and rate of elimination of infected cells by CTLs, based on a fractional-order differential equation model. For the one-virus model, our analysis shows that the interior equilibrium which is unstable in the classical integer-order model can become asymptotically stable in our fractional-order model and numerical simulations confirm this. We also present simulation results of the chaotic behaviors produced from the fractional-order HIV model with viral diversity by using an Adams-type predictor-corrector method.

1. Introduction

An important part of the human immune response against viral infections is cytotoxic T lymphocytes (CTLs) [1]. They recognize and kill cells which are infected by virus. There are many immune models describing the virus dynamics with CTL immune response. Nowak and Bangham [2, 3] proposed an ODE model which explores the relation among CTL immune responses, virus load, and virus diversity. In [2], a rate of specific CTL (𝑍_𝑗) proliferation in response to the corresponding specific infected cells (𝐼_𝑗) depends on the mass action law 𝑐𝐼_𝑗𝑍_𝑗. This model has been important in the field of mathematical modelling of HIV infection. In their model, there is no interaction among different types of CTL (𝑍_𝑗). Iwami et al. [4] assumed that the correlation is incorporated as a function of the frequency that the specific CTLs (𝑍_𝑗) encounter in the specific infected cells (𝐼_𝑗). In a similar manner, they considered the rate of elimination of specific infected cells (𝐼_𝑗) by the specific CTLs (𝑍_𝑗) to be proportional to this frequency. However, these models do not take into account the fractional order derivatives that have been extensively applied in many fields (e.g., [5–17] and the reference cited therein). Recently many mathematicians and applied researchers have tried to model real processes using the fractional order differential equations (FODE) [16]. In biology, it has been deduced that the membranes of cells of biological organism have fractional order electrical conductance [13] and then, they are classified into group of noninteger order models. Also, it has been shown that modelling the behavior of brainstem vestibule-oculumotor neurons by FODE has more advantages than classical integer order modelling [8].

Particular emphasis is that a major difference between fractional order models and integer order models is that fractional order models possess memory [5, 12], while the main features of immune response involve memory [18]. Hence, we attempt to model HIV infection with immune response using a fractional order system. Our presentation is based on the immune model of HIV infection which is developed by Iwami et al. [4]. For the one-virus model, we carry out a detailed analysis on stability of equilibrium. Our analysis shows that the interior equilibrium which is unstable in the classical integer order model can become asymptotically stable in our fractional order model. We also find that chaos does exist in the fractional order HIV model with viral diversity.

2. Model Derivation

We first give the definition of fractional order integration and fractional order differentiation [14, 16]. For the concept of fractional derivative we will adopt Caputo's definition which is a modification of the Riemann-Liouville definition and has the advantage of dealing properly with initial value problems.

Definition 2.1. The fractional integral of order 𝛼>0 of a function 𝑓∶𝑅⁺→𝑅 is given by 𝐼^𝛼1𝑓(𝑥)=Γ(𝛼)^𝑥₀(𝑥−𝑡)^𝛼−1𝑓(𝑡)𝑑𝑡 (2.1) provided the right side is pointwise defined on 𝑅⁺.

Definition 2.2. The Caputo fractional derivative of order 𝛼∈(𝑛−1,𝑛) of a continuous function 𝑓∶𝑅⁺→𝑅 is given by 𝐷^𝛼𝑓(𝑥)=𝐼^𝑛−𝛼𝐷^𝑛𝑑𝑓(𝑥),𝐷=𝑑𝑡. (2.2)

Now we introduce fractional order into the ODE model by Iwami et al. [4]. The new system is described by the following set of FODE: 𝐷^𝑞1𝑇=𝜆−𝑑𝑇−^𝑛_𝑙=1𝛽^′_𝑙𝑇𝑉_𝑙,𝐷^𝑞2𝑗𝐼_𝑗=𝛽^′_𝑗𝑇𝑉_𝑗−𝑎𝐼_𝑗−𝑏𝑍_𝑗𝐼_𝑗𝑇+∑^𝑛_𝑙=1𝐼_𝑙,𝐷^𝑞3𝑗𝑍_𝑗=𝑐𝑍_𝑗𝐼_𝑗𝑇+∑^𝑛_𝑙=1𝐼_𝑙−𝛿𝑍_𝑗,𝐷^𝑞4𝑗𝑉_𝑗=𝑘𝑎𝐼_𝑗−𝑢𝑉_𝑗(𝑗=1,2,…,𝑛), (2.3) where 𝑇(𝑡) represents the concentration of uninfected cells at time 𝑡, 𝐼_𝑗(𝑡) represents the concentration of infected cells with a virus particle of type 𝑗, 𝑉_𝑗(𝑡) the concentration of free virus particle of type 𝑗, and 𝑍_𝑗(𝑡) denotes the magnitude of the specific CTL response against variant 𝑗. Here, 0.95≤𝑞₁,𝑞_2𝑗,𝑞_3𝑗,𝑞_4𝑗≤1(𝑗=1,2,…,𝑛) are restricted such that fractional derivative can be approximately described the rate of change in number.

Following [4], uninfected cells are assumed to be generated at a constant rate 𝜆. Uninfected cells, infected cells, free viruses, and CTLs decline at rates 𝑑,𝑎,𝑢, and 𝛿, respectively. The total number of virus particles produced from one cell is 𝑘. The rate of CTL proliferation in response to antigen is given by 𝑐𝑍_𝑗𝐼_𝑗/(𝑇+∑^𝑛_𝑙=1𝐼_𝑙) and the specific infected cells are killed by specific CTLs at rate 𝑏𝑍_𝑗𝐼_𝑗/(𝑇+∑^𝑛_𝑙=1𝐼_𝑙), while infected cells are produced from uninfected cells and free virus at rate 𝛽^′_𝑙𝑇𝑉_𝑙. That is, a rate of specific CTL(𝑍_𝑗) proliferation in response to the corresponding specific infected cells (𝐼_𝑗) depends on the frequency, instead of the mass action law.

To simplify the model, it is reasonable to assume that the decay rate of free virus, 𝑢, is much larger than that of the infected cells, 𝑎, and this system describes the qualitative dynamics of the asymptomatic phase of HIV infection. Thus, we may introduce as a good approximation that the virus is in steady state (i.e., 𝐷^𝑞4𝑗𝑉_𝑗=0) and hence 𝑉_𝑗=𝑘𝑎𝐼_𝑗/𝑢 (see [4, 19]). This leads to the following simplified system of FODE: 𝐷^𝑞1𝑇=𝜆−𝑑𝑇−^𝑛_𝑙=1𝛽_𝑙𝑇𝐼_𝑙,𝐷^𝑞2𝑗𝐼_𝑗=𝛽_𝑗𝑇𝐼_𝑗−𝑎𝐼_𝑗−𝑏𝑍_𝑗𝐼_𝑗𝑇+∑^𝑛_𝑙=1𝐼_𝑙,𝐷^𝑞3𝑗𝑍_𝑗=𝑐𝑍_𝑗𝐼_𝑗𝑇+∑^𝑛_𝑙=1𝐼_𝑙−𝛿𝑍_𝑗(𝑗=1,2,…,𝑛), (2.4) where 𝛽_𝑗=𝑘𝑎𝛽^′_𝑗/𝑢.

3. One-Virus Model

In this section, we discuss in detail an important special case of model (2.4) and perform an equilibrium and stability analysis for this special case. We consider the one-virus model (𝑛=1) and assume that 𝑞₁=𝑞₂₁=𝑞₃₁=𝛼(0.95≤𝛼≤1). This one-virus model is described by the following system of FODE: 𝐷^𝛼𝑇=𝜆−𝑑𝑇−𝛽₁𝑇𝐼₁,𝐷^𝛼𝐼₁=𝛽₁𝑇𝐼₁−𝑎𝐼₁−𝑏𝑍₁𝐼₁𝑇+𝐼₁,𝐷^𝛼𝑍₁=𝑐𝑍₁𝐼₁𝑇+𝐼₁−𝛿𝑍₁. (3.1)

To evaluate the equilibria, let 𝐷^𝛼𝐷𝑇=0,_𝛼𝐼₁𝐷=0,_𝛼𝑍₁=0. (3.2) Then system (3.1) has three equilibria: the uninfected equilibrium 𝐸_𝐻=(𝜆/𝑑,0,0), the boundary equilibrium 𝐸_𝐼=(𝑇^∗,𝐼^∗₁,0), where 𝑇^∗=𝑎𝛽₁,𝐼^∗₁=𝜆𝑎−𝑑𝛽₁, (3.3) and the interior equilibrium 𝐸_𝑐=(𝑇,𝐼₁̂,𝑍₁), where 𝑇=−𝑑+𝑑²+4𝜆𝛽2𝛽,𝐼₁̂=𝛿𝑍𝑐−𝛿𝑇,₁=𝑐𝑇𝑏(𝑐−𝛿)(𝛽₁𝑇−𝑎),𝛽=𝛿𝛽₁.𝑐−𝛿 (3.4) Following the analysis in [4], we introduce a basic reproduction number which is defined by 𝑅₀=𝜆𝛽₁𝑎𝑑. (3.5)

Denote 𝑅³₊={(𝑇,𝐼₁,𝑍₁)∈𝑅³∣𝑇≥0,𝐼₁≥0,𝑍₁≥0} and we always assume that 𝑐>𝛿. Note that 𝐷^𝛼𝑍₁<0 always holds true if 𝑐≤𝛿. By generalized mean value theorem [15], we get 𝑍₁(𝑡) is decreasing if 𝑐≤𝛿.

Next we will discuss the existence and stability of the equilibria of the model (3.1).

Theorem 3.1. (a) The uninfected equilibrium 𝐸_𝐻 is locally asymptotically stable (LAS) if 𝑅₀∈(0,1) and unstable if 𝑅₀>1.
(b) If 𝑅₀>1, then the boundary equilibrium 𝐸_𝐼 exists. This equilibrium is LAS if 1<𝑅₀<(𝑎𝛿/𝑑(𝑐−𝛿))+1 and unstable if 𝑅₀>(𝑎𝛿/𝑑(𝑐−𝛿))+1.(c) If 𝑅₀>(𝑎𝛿/𝑑(𝑐−𝛿))+1, then 𝐸_𝑐 exists in Int𝑅³₊, where Int𝑅³₊ is the interior of 𝑅³₊.

Proof. (a) The Jacobian matrix 𝐽(𝐸_𝐻) for system (3.1) evaluated at 𝐸_𝐻 is given by 𝐽(𝐸_𝐻𝛽)=−𝑑−¹𝜆𝑑00𝛽₁𝜆𝑑−𝑎000−𝛿. (3.6)𝐸_𝐻 is locally asymptotically stable if all of the eigenvalues 𝑝 of the Jacobian matrix 𝐽(𝐸_𝐻) satisfy the following condition [6, 17]: |arg(𝑝)|>𝛼𝜋2. (3.7) The eigenvalues of 𝐽(𝐸_𝐻) are −𝑑,𝛽₁𝜆/𝑑−𝑎,−𝛿. It is clear that 𝐸_𝐻 is LAS if 𝑅₀<1 and is unstable if 𝑅₀>1.
(b) If 𝑅₀>1, then the existence of 𝐸_𝐼 is obvious.
The Jacobian matrix 𝐽(𝐸_𝐼) for system (3.1) evaluated at 𝐸_𝐼 is given by 𝐽(𝐸_𝐼)=−𝑑−𝛽₁𝐼^∗₁−𝛽₁𝑇^∗0𝛽₁𝐼^∗₁0−𝑏𝐼^∗₁𝑇^∗+𝐼^∗₁00𝑐𝐼^∗₁𝑇^∗+𝐼^∗₁−𝛿. (3.8) For 𝐽(𝐸_𝐼) given by (3.8), the characteristic equation becomes [𝑝²+(𝑑+𝛽₁𝐼^∗₁)𝑝+𝛽²₁𝐼^∗₁𝑇^∗](𝑝−𝑐𝐼^∗₁𝑇^∗+𝐼^∗₁+𝛿)=0. (3.9) and hence all the eigenvalues are 𝑝_1,2=−(𝑑+𝛽₁𝐼^∗₁)±(𝑑+𝛽₁𝐼^∗₁)²−4𝛽²₁𝐼^∗₁𝑇^∗2,𝑝₃=𝑐𝐼^∗₁𝑇^∗+𝐼^∗₁−𝑑. (3.10) If 𝑅₀>1, then 𝑇^∗>0,𝐼^∗₁>0, and 𝑝_1,2 have negative real parts. Furthermore, if 1<𝑅₀<(𝑎𝛿/𝑑(𝑐−𝛿))+1, then𝑝₃<0 and 𝐸_𝐼 is LAS. If 𝑅₀>(𝑎𝛿/𝑑(𝑐−𝛿))+1, then 𝑝₃>0 and 𝐸_𝐼 is unstable.
(c) If 𝑅₀>(𝑎𝛿/𝑑(𝑐−𝛿))+1, then we obtain 𝑍₁>0. Thus, 𝐸_𝑐 exists in Int𝑅³₊. Therefore, the proof is complete.

To discuss the local stability of the interior equilibrium 𝐸_𝑐, we consider the linearized system of (3.1) at 𝐸_𝑐. The Jacobian matrix at 𝐸_𝑐 is given by 𝐽(𝐸_𝑐)=−𝑑−𝛽₁𝐼₁̂−𝛽₁𝛽𝑇0₁𝐼₁̂+𝑏𝑍₁𝐼₁̂(𝑇+𝐼₁̂)²𝛽₁𝑏𝑍𝑇−𝑎−¹𝑇(𝑇+𝐼₁̂)²−𝑏𝐼₁̂𝑇+𝐼₁̂−𝑐𝑍₁𝐼₁̂(𝑇+𝐼₁̂)²𝑐𝑍₁𝑇(𝑇+𝐼₁̂)²𝑐𝐼₁̂𝑇+𝐼₁̂−𝛿. (3.11) For convenience, we denote 𝑇=𝑇,𝐼₁̂=𝐼,𝑍₁=𝑍, and 𝛽₁=𝛽. In view of the above assumptions and using 𝛽𝑇−𝑎−𝑏𝑍𝑇+𝐼=0,𝜆−𝑑𝑇−𝛽𝑇𝐼=0,𝑐𝐼𝑇𝑇+𝐼=𝛿,=𝑇+𝐼𝑐−𝛿𝑐, (3.12)𝐽(𝐸_𝑐) can now be written as follows: 𝐽(𝐸_𝑐𝛿)=−𝑑−𝛽𝐼−𝛽𝑇0𝛽𝐼+𝑐𝛿(𝛽𝑇−𝑎)𝑐(𝛽𝑇−𝑎)−𝑏𝛿𝑐−𝛿𝑏(𝛽𝑇−𝑎)𝑐−𝛿𝑏(𝛽𝑇−𝑎)0. (3.13) Then the characteristic equation of the linearized system of (3.1) is Φ(𝑝)=𝑝³+𝑎₁𝑝²+𝑎₂𝑝+𝑎₃=0, (3.14) where 𝑎₁=𝑑+𝛿𝑎𝑐+𝛿²𝛽𝑇,𝑎𝑐(𝑐−𝛿)₂𝛿=𝛽𝑇(𝛽𝐼+𝑐𝛿(𝛽𝑇−𝑎))+𝑐𝛿(𝑐−𝛿)(𝛽𝑇−𝑎)−𝑐𝑎(𝑑+𝛽𝐼)(𝛽𝑇−𝑎),₃=𝛿𝑐(𝑐−𝛿)(𝛽𝑇−𝑎)(𝑑+𝛽𝐼)+𝛽𝛿²𝑐𝑇(𝛽𝑇−𝑎). (3.15)

Proposition 3.2. The interior equilibrium 𝐸_𝑐 is LAS if all of the eigenvalues 𝑝 of 𝐽(𝐸_𝑐) satisfy |arg(𝑝)|>𝛼𝜋/2.
Denote 𝐷(Φ)=−1𝑎₁𝑎₂𝑎₃001𝑎₁𝑎₂𝑎₃32𝑎₁𝑎₂00032𝑎₁𝑎₂00032𝑎₁𝑎₂=18𝑎₁𝑎₂𝑎₃+(𝑎₁𝑎₂)²−4𝑎₃𝑎³₁−4𝑎³₂−27𝑎²₃. (3.16)

Using the results of [5, 20], we have the following proposition.

Proposition 3.3. One assumes that 𝐸_𝑐 exists in Int𝑅³₊.
(i)If the discriminant of Φ(𝑝),𝐷(Φ) is positive and Routh-Hurwitz conditions are satisfied, that is, 𝐷(Φ)>0,𝑎₁>0,𝑎₃>0,𝑎₁𝑎₂>𝑎₃, (3.17) then the interior equilibrium 𝐸_𝑐 is LAS.(ii)If 𝐷(Φ)<0,𝑎₁>0,𝑎₂>0,𝑎₁𝑎₂=𝑎₃,𝛼∈[0,1), then the interior equilibrium 𝐸_𝑐 is LAS.(iii)If 𝐷(Φ)<0,𝑎₁<0,𝑎₂<0,𝛼>2/3, then the interior equilibrium 𝐸_𝑐 is unstable.

In our first example we set 𝜆=10,𝑑=0.02 which are chosen according to [21] and set 𝑎=𝛿=0.04,𝑐=𝑏=0.8,𝛽₁=4×10⁻⁴ which come from [4]. With these parameter values, 𝑅₀=5>(𝑎𝛿/𝑑(𝑐−𝛿))+1,𝐷(Φ)=−9.4073×10⁻⁷<0,𝑎₁𝑎₂<𝑎₃. By Proposition 3.2., we obtain the interior equilibrium 𝐸_𝑐=(362.0335,19.0544,49.9289) is LAS when 𝛼<0.9916. Numerical simulations show that trajectories of system (3.1) approach to the interior equilibrium (see Figures 1(a) and 1(b)). However, when 𝛼=1 (that is the case of classical integer order), 𝐸_𝑐 is unstable by the Routh-Hurwitz criterion(see Figures 2(a) and 2(b)).

(a)

(b)

Figure 1 

Numerical solutions of system (3.1). The plots show that trajectories of system (3.1) approach to the interior equilibrium for 𝛼=0.95.

(a)

(b)

Figure 2 

Numerical solutions of system (3.1). The plots show that the interior equilibrium is unstable for 𝛼=1.

4. Two-Virus Model

In this section, we consider viral diversity. We examine the two-virus model using numerical simulations. By examining the behavior of this simpler model we hope to get an idea as to how the more general models in system (2.4) may behave. The two-virus model is given by the following system of FODE: 𝐷^𝑞1𝑇=𝜆−𝑑𝑇−𝛽₁𝑇𝐼₁−𝛽₂𝑇𝐼₂,𝐷^𝑞21𝐼₁=𝛽₁𝑇𝐼₁−𝑎𝐼₁−𝑏𝑍₁𝐼₁𝑇+𝐼₁+𝐼₂,𝐷^𝑞22𝐼₂=𝛽₂𝑇𝐼₂−𝑎𝐼₂−𝑏𝑍₂𝐼₂𝑇+𝐼₁+𝐼₂,𝐷^𝑞31𝑍₁=𝑐𝑍₁𝐼₁𝑇+𝐼₁+𝐼₂−𝛿𝑍₁,𝐷^𝑞32𝑍₂=𝑐𝑍₂𝐼₂𝑇+𝐼₁+𝐼₂−𝛿𝑍₂, (4.1) with initial value condition 𝑇(0)=𝑇₀,𝐼_𝑖(0)=𝐼_𝑖0,𝑍_𝑖(0)=𝑍_𝑖0,𝑖=1,2, (4.2) where 0.95≤𝑞₁,𝑞_2𝑗,𝑞_3𝑗(𝑗=1,2)≤1.

To find numerical solution to (4.1) and (4.2) in the interval [0,𝑇], we reduce the systems (4.1) and (4.2) to a set of fractional integral equations, by using an equivalence (see [16, Theorem 3.24]) 𝐷^𝛼𝑋=𝑓(𝑋)⟺𝑋(𝑡)=𝑋(0)+𝐼^𝛼𝑓(𝑋). (4.3) Then we apply the generalized Adams-type predictor-corrector method or, more precisely, Predict, Evaluate, Correct, Evaluate (PECE) methods (see [22, 23]).

For notational convenience, we denote 𝛼=(𝑞₁,𝑞₂₁,𝑞₂₂,𝑞₃₁,𝑞₃₂). We carry out numerical simulations for system, (4.1) and (4.2) with parameters 𝜆=10,𝑏=𝑐=0.8,𝑎=0.031,𝛿=0.03,𝑑=0.02,𝛽₁=4×10⁻⁴, and 𝛽₂=2.8×10⁻⁴ for the step size 0.07. Numerical solutions of systems (4.1) and (4.2) support that the system exhibits a chaotic behavior and systems (4.1) and (4.2) have a strange attractor in Int𝑅⁵₊ for 𝛼=[0.95,0.95,1,1,1] (see Figures 3(a)–3(c)). It is clear that chaos does exist in our fractional order model with viral diversity as in the case of integer order model. The effect of viral diversity and the frequency dependence results in collapse of the immune system and make the behavior of the system dynamics complex [4]. However, as the value of some component or more components of the order 𝛼 further decreases, for example, 𝛼=[0.95,0.95,0.95,1,1], the chaotic motion disappears and the systems (4.1) and (4.2) stabilize to a fixed point (see Figures 4(a)–4(c)).

(a)

(b)

(c)

Figure 3 

Numerical solutions of system (4.1) for 𝛼=[0.95,0.95,1,1,1]. (a) A strange attractor in the 𝑍1-𝑍2-𝐼1 phase. (b) Infected cell 1. (c) CTL 2.

(a)

(b)

(c)

Figure 4 

Numerical solutions of system (4.1) for 𝛼=[0.95,0.95,0.95,1,1]. (a) 𝑍1-𝑍2-𝐼1 phase. (b) Infected cell 1. (c) CTL 2.

5. Conclusions

In this paper, we have proposed a fractional order HIV model, as a generalization of an integer order model, developed by Iwami et al. [4]. The premise of the proposed model is the fact that fractional order models possess memory while the main features of immune response involve memory. It is an attempt to incorporate fractional order into the mathematical model of HIV-immune system dynamics and it is still an interesting exercise to determine, mathematically, how the order of a fractional differential system affects the dynamics of system.

In the case of one-virus model, the fractional order system has an interior equilibrium under some restriction. By using stability analysis on fractional order system, we obtain sufficient condition on the parameters for the stability of the interior equilibrium. Our analysis shows that the interior equilibrium which is unstable in the classical integer order model can become asymptotically stable in our fractional order model. Note that the interior equilibrium is globally asymptotically stable (GAS) (see [24]) if the terms associated with immune reactions are given by 𝑐𝑍₁𝐼₁ and 𝑏𝑍₁𝐼₁ instead of 𝑐𝑍₁𝐼₁/(𝑇+𝐼₁) and 𝑏𝑍₁𝐼₁/(𝑇+𝐼₁) in (3.1). That is, the interior equilibrium of the one-virus model can become unstable because of the frequency dependence (see [4]). However, in our fractional order model with the frequency dependence, the interior equilibrium can also become asymptotically stable if the order 𝛼<0.9916.

We then consider viral diversity. If the terms associated with immune reaction depend on the mass action law instead of frequency, an interior equilibrium in [24] is GAS. Similar to the integer order model in [4], we find that strange chaotic attractors can be obtained under fractional order model with frequency dependence. That is, the effect of viral diversity and the frequency dependence results in collapse of the immune system and make the behavior of the system dynamics complex. However the chaotic motion may disappear and the fractional order system stabilizes to a fixed point if the value of the order 𝛼 decreases. The specific biological meaning is deserved to further study.

Acknowledgments

The authors are very grateful to the referees for their valuable suggestions, which helped to improve the paper significantly. This work was supported in part by Specialized Research Fund for the Doctoral Program of Higher Education from Ministry of Education of China (no. 20060255006), Project of the Shanghai Committee of Science and Technology (no. 08JC1400100), and the Open Fund from the Key Laboratory of MICCAI of Shanghai (06dz22103).