Stability Analysis of an HIV/AIDS Dynamics Model with Drug Resistance

Abstract

A mathematical model of HIV/AIDS transmission incorporating treatment and drug resistance was built in this study. We firstly calculated the threshold value of the basic reproductive number (𝑅₀) by the next generation matrix and then analyzed stability of two equilibriums by constructing Lyapunov function. When 𝑅₀<1, the system was globally asymptotically stable and converged to the disease-free equilibrium. Otherwise, the system had a unique endemic equilibrium which was also globally asymptotically stable. While an antiretroviral drug tried to reduce the infection rate and prolong the patients’ survival, drug resistance was neutralizing the effects of treatment in fact.

1. Introduction

It was reported that 2.7 million people were newly infected by HIV/AIDS virus and 1.8 million patients died of AIDS-related causes in 2010 worldwide. By the end of 2010, about 34 million people were living with HIV/AIDS in the world [1]. China estimated that 2.8 million died of AIDS-related causes in 2011, and there were about 7.8 million HIV-infected people by the end of 2011 [2].

Since the initial infectious diseases model was presented by Anderson et al. in 1986 [3–5], various mathematical models have been developed among which the treatment has been addressed [6–16]. For example, Wang and Zhou’s model tried to address HIV treatment and progression by CD4 ~+ T-cells and virus particles in microcosmic [8]; Blower, Boily et al., and Bachar and Dorfmayr’s works tried to investigate the effect of treatment on sexual behaviors [9–11]; Blower et al., Sharomi and Gumel and Nagelkerke et al.studied the epidemic contagion transmission in some specific regions or groups considering drug resistance [12–14].

In this work, we established a model by adding effect of drug resistance into the similar models in the literatures [6–11, 15, 16]. The model in [12–14] are include both virus drug resistance and drug sensitive on the treatment. However the model in [12] did not distinguish the stage of HIV and AIDS, our model aware of the different between HIV infections and AIDS patients. Moreover, compare to those models in [13, 14], we carefully consider some infections would exit treatment group without developing drug-resistance due to other reasons such as migration. Theoretical analysis on global stability of endemic equilibrium has then been implemented.

2. Dynamic Model

According to the progression of disease, the total populations were separated into five groups: susceptible population, early-stage HIV population, symptomatic population, AIDS patients, and those who are accepting ART; we marked them with 𝑆(𝑡), 𝐼(𝑡), 𝐽(𝑡), 𝐴(𝑡), and 𝑇(𝑡) separately. Treatment has three outcomes: (1) a patient can respond to treatment and remain the ART; (2) exit treatment due to clinical failure, migration, or other reasons without developing drug resistance; (3) virologically fail and develop drug resistance. We use 𝑅(𝑡) to denote patients in situation (3).

We made some assumption as below.(1)Infection occurred when susceptible and infected contact with each other took place.(2)Only people in the period of AIDS may die of AIDS disease-related, then use 𝑑 denote the disease-related death rate of the AIDS. And use 𝜇 denote the mortality rate in the total population.(3)Although AIDS patients have the higher viral load, we assume they will not infect others because they have the obvious clinical symptoms and were accepting ART.(4)When in the situation (2) of treatment, we assume these people transformed into asymptomatic individuals.

According to the assumptions, the flow diagram of the six subpopulations was shown in Figure 1.

Figure 1 

Relationship between different populations.

The model was collected as the following differential equations: 𝑆^′=Λ−𝜇𝑆−𝛽₁𝐼+𝛽₂𝐽+𝛽₃𝑇+𝛽₄𝑅𝐼𝑆,_′=𝛽₁𝐼+𝛽₂𝐽+𝛽₃𝑇+𝛽₄𝑅𝐽𝑆−𝛼𝐼−𝜇𝐼,_′=𝛼𝐼−𝜌₁+𝜎+𝑘₁𝑇𝐽−𝜇𝐽+𝛾𝑇,_′=𝜎𝐽−𝜌₂+𝛾+𝑘₂𝑅𝑇−𝜇𝑇,_′=𝑘₁𝐽+𝑘₂𝑇−𝜌₃𝐴𝑅−𝜇𝑅,_′=𝜌₁𝐽+𝜌₂𝑇+𝜌₃𝑅−𝑑+𝜇𝐴, (2.1) where Λ is the recruitment rate of the susceptible population, 𝛽₁ is the probability of transmission by an infection in the first stage, 𝛽₂ is the probability of transmission by an infection in the second phase, 𝛽₃=𝑙𝛽₂(𝑙<1) is the probability of transmission by a patient being treated, 𝛽₄ is the probability of transmission by a drug resistance individual, 𝛼 is the transfer rate constant from the asymptomatic phase 𝐼 to the symptomatic phase 𝐽, 𝜎 is the proportion of treatment, 𝛾 is the exit treatment rate without developing drug resistance, 𝑘₁ is probability constant of infection by transmission of drug-resistant strains, and 𝑘₂ is the rate of acquiring drug resistance during treatment; 𝜌_𝑖(𝑖=1,2,3) denote transfer rate constant by an infection from phase 𝐽, 𝑇, 𝑅 to the AIDS cases 𝐴, respectively.

The model is established in practice; thus we assume all parameters are nonnegative.

Since the 𝐴 of system (2.1) does not appear in the equations, in the following analysis, we only consider the system as follows: 𝑆^′=Λ−𝜇𝑆−𝛽₁𝐼+𝛽₂𝐽+𝛽₃𝑇+𝛽₄𝑅𝐼𝑆,_′=𝛽₁𝐼+𝛽₂𝐽+𝛽₃𝑇+𝛽₄𝑅𝐽𝑆−𝛼+𝜇𝐼,_′=𝛼𝐼−𝜌₁+𝜎+𝑘₁𝑇+𝜇𝐽+𝛾𝑇,_′=𝜎𝐽−𝜌₂+𝛾+𝑘₂𝑅+𝜇𝑇,_′=𝑘₁𝐽+𝑘₂𝑇−𝜌₃+𝜇𝑅. (2.2)

Theorem 2.1. Let the initial data be 𝑆(0)=𝑆₀>0, 𝐼(0)=𝐼₀>0, 𝐽(0)=𝐽₀>0, 𝑇(0)=𝑇₀>0 and 𝑅(0)=𝑅₀>0; then the solutions of system (2.2) are all positive for all t > 0. For the model, the feasible region of system (2.2) is Ω={(𝑆,𝐼,𝐽,𝑇,𝑅)∈ℝ⁵₊∶𝑆+𝐼+𝐽+𝑇+𝑅≤Λ/𝜇,0<𝑆≤Λ/𝜇}, and Ω for system (2.2) is positively invariant.

Proof. From the first equation of (2.2) 𝑆^′=Λ−𝜇𝑆−𝛽₁𝐼+𝛽₂𝐽+𝛽₃𝑇+𝛽₄𝑅𝑆, (2.3) consider the following two categories.(1) When 𝑡₀>0and𝑆(𝑡₀)=0.   Equation (2.3) becomes 𝑆^′=Λ(𝑡=𝑡₀), due to Λ>0; we have 𝑆^′(𝑡₀)>0, that is, when 𝑡>𝑡₀,𝑆(𝑡) is an increasing function about 𝑡. Therefore, we can conclude that when 𝑡 is the neighborhood of 𝑡₀ and 𝑡>𝑡₀,𝑆(𝑡)≥𝑆(𝑡₀)>0.(2)When 𝑡₀>0and𝑆(𝑡₀)>0.  Equation (2.3) can be written as 𝑆^′𝑆=Λ𝑆−𝜇+𝛽₁𝐼+𝛽₂𝐽+𝛽₃𝑇+𝛽₄𝑅, (2.4) that is 𝑆^′𝑆≥−𝜇+𝛽₁𝐼+𝛽₂𝐽+𝛽₃𝑇+𝛽₄𝑅, (2.5) thus, 𝑆(𝑡)≥𝑆(0)exp−^𝑡₀𝜇+𝛽₁𝐼(𝑠)+𝛽₂𝐽(𝑠)+𝛽₃𝑇(𝑠)+𝛽₄𝑅(𝑠)𝑑𝑠>0. (2.6)
Similarly for the other equations of system (2.2) we can easily show that 𝐼(𝑡), 𝐽(𝑡), 𝑇(𝑡), 𝑅(𝑡) are increasing functions about 𝑡 when 𝑡₀>0 and 𝐼(𝑡₀)=0, 𝑡₀>0 and 𝐽(𝑡₀)=0, 𝑡₀>0 and 𝑇(𝑡₀)=0,𝑡₀>0 and 𝑅(𝑡₀)=0, respectively; then when 𝑡>𝑡₀, 𝑋(𝑡)≥𝑋(𝑡₀)>0(𝑋=𝐼,𝐽,𝑇,𝑅). Otherwise, when 𝑡₀>0 and 𝑋(𝑡₀)>0  (𝑋=𝐼,𝐽,𝑇,𝑅), we have the following results corresponding to the respective hypothesis: 𝐼^′𝐼≥−𝛼+𝜇. (2.7) Thus, 𝐼(𝑡)(0)−𝑡𝐽≥𝐼exp𝛼+𝜇>0,_′𝐽≥−𝜌₁+𝜎+𝑘₁+𝜇. (2.8) Thus, 𝐽(𝑡)≥𝐽(0)exp−𝜌₁+𝜎+𝑘₁𝑡𝑇+𝜇>0,_′𝑇≥−𝜌₂+𝛾+𝑘₂+𝜇. (2.9) Thus, 𝑇(𝑡)≥𝑇(0)exp−𝜌₂+𝛾+𝑘₂𝑡𝑅+𝜇>0,_′𝑅≥−𝜌₃+𝜇. (2.10) Thus, 𝑅(𝑡)≥𝑅(0)exp−𝜌₃+𝜇𝑡>0. (2.11) Then, we have that 𝐼(𝑡), 𝐽(𝑡), 𝑇(𝑡), and 𝑅(𝑡) are all strictly positive for 𝑡>0. Thus we can conclude that all solutions of system (2.2) remain positive for all 𝑡>0.
Next, add all the the equations of system (2.2); we have (𝑆+𝐼+𝐽+𝑇+𝑅)^′=Λ−𝜇(𝑆+𝐼+𝐽+𝑇+𝑅)−𝜌₁𝐽+𝜌₂𝑇+𝜌₃𝑅≤Λ−𝜇(𝑆+𝐼+𝐽+𝑇+𝑅). (2.12)
Then, lim_𝑡→∞()≤Λsup𝑆+𝐼+𝐽+𝑇+𝑅𝜇. (2.13)
In a similar fashion we have 𝑆^′≤Λ−𝜇𝑆 from the first equation of (2.2); then lim_𝑡→∞sup𝑆≤Λ/𝜇.
Thus, the feasible solution of system remains in the region Ω, and Ω as the feasible region for system is positively invariant. In the following, the dynamics of system (2.2) will be considered in Ω.

3. The Basic Reproduction Number and the Disease-Free Equilibrium

3.1. The Basic Reproduction Number

It is easy to see that the model has a disease-free equilibrium (DFE),  𝑃₀=(Λ/𝜇,0,0,0,0). Following the paper [17], we obtain the basic reproduction number by using the next generation operator approach.

Let 𝑦=(𝐼,𝐽,𝑇,𝑅,𝑆)^𝑇; thus we have 𝑦^′=ℱ𝑦−𝒱𝑦, (3.1) where ℱ𝑦=𝛽₁𝐼+𝛽₂𝐽+𝛽₃𝑇+𝛽₄𝑅𝑆0000𝒱𝑦=𝐼𝜌,𝛼+𝜇−𝛼𝐼+₁+𝜎+𝑘₁𝜌+𝜇𝐽−𝛾𝑇−𝜎𝐽+₂+𝛾+𝑘₂+𝜇𝑇−𝑘₁𝐽−𝑘₂𝑇+𝜌₃𝑅𝛽+𝜇−Λ+𝜇𝑆+₁𝐼+𝛽₂𝐽+𝛽₃𝑇+𝛽₄𝑅𝑆. (3.2)(3.3) Then the derivatives of ℱ(𝑦) and 𝒱(𝑦) at the DFE 𝑃₀=(0,0,0,0,Λ/𝜇) are partitioned as 𝐷ℱ𝑃₀=,𝐹000𝐷𝒱𝑃₀=𝑉0000𝛽₁⋅Λ𝜇𝛽₂⋅Λ𝜇𝛽₃⋅Λ𝜇𝛽₄⋅Λ𝜇𝜇. (3.4)𝐹 and 𝑉 are the 4×4 matrices as follows: 𝛽𝐹=¹⋅Λ𝜇𝛽₂⋅Λ𝜇𝛽₃⋅Λ𝜇𝛽₄⋅Λ𝜇0000,𝑄00000000𝑉=¹000−𝛼𝑄₂−𝛾00−𝜎𝑄₃00−𝑘₁−𝑘₂𝑄₄, (3.5) where 𝑄₁=𝛼+𝜇,𝑄₂=𝜌₁+𝜎+𝑘₁,𝑄+𝜇₃=𝜌₂+𝛾+𝑘₂+𝜇,𝑄₄=𝜌₃+𝜇. (3.6)

Hence the reproduction number, denoted by 𝑅₀, is the spectral radius of the next generation matrix 𝐹𝑉⁻¹: 𝑅₀=𝜌𝐹𝑉⁻¹=𝑅₁+𝑅₂+𝑅₃+𝑅₄. (3.7) Here 𝑅₁=𝛽₁𝑄₁⋅Λ𝜇,𝑅₂=𝛽₂𝛼𝑄₃𝑄₁𝑄₂𝑄₃⋅Λ−𝜎𝛾𝜇,𝑅₃=𝛽₃𝛼𝜎𝑄₁𝑄₂𝑄₃⋅Λ−𝜎𝛾𝜇,𝑅₄=𝛽₄𝛼𝑘₁𝑄₃+𝑘₂𝜎𝑄₁𝑄₂𝑄₃−𝜎𝛾𝑄₄⋅Λ𝜇. (3.8)

3.2. DFE and Stability

Theorem 3.1. The disease-free equilibrium 𝑃₀ of system is globally asymptotically stable for 𝑅₀<1 and unstable for 𝑅₀>1.

Proof. (1) The Jacobian matrices of system (2.2) at the DFE are 𝐽𝑃₀=−𝜇−𝛽₁⋅Λ𝜇−𝛽₂⋅Λ𝜇−𝛽₃⋅Λ𝜇−𝛽₄⋅Λ𝜇0000𝐷, (3.9) where 𝛽𝐷=¹⋅Λ𝜇−𝑄₁𝛽₂⋅Λ𝜇𝛽₃⋅Λ𝜇𝛽₄⋅Λ𝜇𝛼−𝑄₂𝛾00𝜎−𝑄₃00𝑘₁𝑘₂−𝑄₄. (3.10)
Obviously, −𝐷=[𝑑_𝑖𝑗] is a 4×4 matrix with 𝑑_𝑖𝑗≤0, for 𝑖≠𝑗,𝑖,𝑗=1,…,4 and 𝑑_𝑖𝑖>0 for 𝑖=2,…,4. If 𝑅₀<1, we have 𝑅₁=(𝛽₁/𝑄₁)⋅(Λ/𝜇)<𝑅₀<1 from the expression of (3.8); thus 𝑑₁₁=𝑄₁−𝛽₁⋅(Λ/𝜇)=𝑄₁(1−𝑅₁)>0.
Define the positive vector subsequently: 𝑥=1,𝛼𝑄₃𝑄₂𝑄₃,−𝜎𝛾𝛼𝜎𝑄₂𝑄₃,𝛼𝑘−𝜎𝛾¹𝑄₃+𝑘₂𝜎𝑄₂𝑄₃−𝜎𝛾𝑄₄^𝑇. (3.11) If 𝑅₀<1, −𝐷⋅𝑥=[𝑄₁(1−𝑅₀),0,0,0]^𝑇≥0.
Then begin to show that all the eigenvalues of −𝐷 are nonzero: det(−𝐷)=𝑄₁𝑄₄𝑄₂𝑄₃−𝜎𝛾−𝛽₁𝑄₂𝑄₃−𝜎𝛾𝑄₄−𝛼𝛽₂𝑄₃+𝛽₃𝜎𝑄₄−𝛼𝛽₂𝑄₃+𝛽₃𝜎𝑄₄. (3.12) Simplify the above expression through substituting formula (3.8) into (3.12): det(−𝐷)=𝑄₁𝑄₄𝑄₂𝑄₃−𝜎𝛾1−𝑅₀. (3.13) Here 𝑄₂𝑄₃−𝜎𝛾=(𝜌₁+𝜎+𝑘₁+𝜇)(𝜌₂+𝛾+𝑘₂+𝜇)−𝜎𝛾>0; thus, det(−𝐷)>0; namely, −𝐷 has non zero eigenvalue for 𝑅₀<1.
In conclusion, if 𝑅₀<1,−𝐷 is an irreducible matrix with 𝑑_𝑖𝑖>0 and 𝑑_𝑖𝑗≤0(𝑖≠𝑗), there exists a positive vector 𝑥 such that −𝐷⋅𝑥≥0. Hence, the real part of each nonzero eigenvalue of −𝐷 is positive according to the M-matrix theory; that is, each eigenvalue of 𝐷 has negative real part.
Through the structure of the Jacobian matrix 𝐽(𝑃₀), it can be seen that the eigenvalues of 𝐽(𝑃₀) consist of −𝜇 and all eigenvalues of 𝐷. Hence, all eigenvalues of 𝐽(𝑃₀) have negative real part for 𝑅₀<1; thus, disease-free equilibrium 𝑃₀ is locally asymptotically stable.
(2) Let 𝑈₁=𝑚₁𝐼+𝑚₂𝐽+𝑚₃𝑇+𝑚₄𝑅, where 𝑚_𝑖(𝑖=1⋯4) are positive constants as follows: 𝑚₁=𝛽₁𝑄₁+𝛽₂𝛼𝑄₃𝑄₁𝑄₂𝑄₃+𝛽−𝜎𝛾³𝛼𝜎𝑄₁𝑄₂𝑄₃+𝛽−𝜎𝛾⁴𝛼𝑘₁𝑄₃+𝑘₂𝜎𝑄₁𝑄₂𝑄₃−𝜎𝛾𝑄₄,𝑚₂=𝛽₂𝑄₃𝑄₂𝑄₃+𝛽−𝜎𝛾³𝜎𝑄₂𝑄₃+𝛽−𝜎𝛾⁴𝑘₁𝑄₃+𝑘₂𝜎𝑄₂𝑄₃−𝜎𝛾𝑄₄,𝑚₃=𝛽₂𝛾𝑄₂𝑄₃+𝛽−𝜎𝛾³𝑄₂𝑄₂𝑄₃+𝛽−𝜎𝛾⁴𝑘₁𝛾+𝑘₂𝑄₂𝑄₂𝑄₃−𝜎𝛾𝑄₄,𝑚₄=𝛽₄𝑄₄. (3.14) When 𝑅₀<1, the time derivative of 𝑈₁ is 𝑈^′_1(2.2)=𝑚₁𝐼^′+𝑚₂𝐽^′+𝑚₃𝑇^′+𝑚₄𝑅^′=𝑚₁𝛽₁𝑆−𝑚₁𝑄₁+𝑚₂𝛼𝐼+𝑚₁𝛽₂𝑆−𝑚₂𝑄₂+𝑚₃𝜎+𝑚₄𝑘₁𝐽+𝑚₁𝛽₃𝑆+𝑚₂𝛾−𝑚₃𝑄₃+𝑚₄𝑘₂𝑇+𝑚₁𝛽₄𝑆−𝑚₄𝑄₄𝑅=𝑚₁𝑆−1𝛽₁𝐼+𝛽₂𝐽+𝛽₃𝑇+𝛽₄𝑅. (3.15) Let 𝑀=max(𝛽₁/𝑚₁,𝛽₂/𝑚₂,𝛽₃/𝑚₃,𝛽₄/𝑚₄); due to lim_𝑡→∞sup𝑆≤Λ/𝜇, 𝑀is the finite number, then 𝑈^′₁≤𝑚₁Λ𝜇−1𝛽₁𝐼+𝛽₂𝐽+𝛽₃𝑇+𝛽₄𝑅=𝑅₀−1𝛽₁𝐼+𝛽₂𝐽+𝛽₃𝑇+𝛽₄𝑅≤𝑀𝑅₀−1𝑚₁𝐼+𝑚₂𝐽+𝑚₃𝑇+𝑚₄𝑅. (3.16) Solve 𝑈^′₁≤𝑀(𝑅₀−1)𝑈₁, have 𝑈₁≤𝑈₁(0)𝑒^{𝑀(𝑅0−1)𝑡}, that is, lim_𝑡→∞𝑈₁(𝑡)=0, thus when 𝑡→∞, (𝐼,𝐽,𝑇,𝑅)→(0,0,0,0). When 𝑅₀<1, 𝑈^′₁≤0, and equalities hold if and only if 𝑆=Λ/𝜇,𝐼=𝐽=𝑇=𝑅=0, that is 𝑈^′₁=0 if and only if 𝑡→∞. We can conclude that the solutions of system (2.2) are all in Ψ={(𝑆,𝐼,𝐽,𝑇,𝑅)∶𝑆=Λ/𝜇,𝐼=𝐽=𝑇=𝑅=0} and the only invariant set in Ψ is 𝑃₀ by the LaSalle’s invariance principle. Thus the solutions of system (2.2) are limits to the endemic equilibrium 𝑃₀ when 𝑅₀<1. Combine locally asymptotically stable of 𝑃₀ with convergence properties of the 𝑃₀, we conclude that 𝑃₀ of system is globally asymptotically stable for 𝑅₀<1.
(3) If 𝑅₀>1, det𝐷=𝑄₁𝑄₄(𝑄₂𝑄₃−𝜎𝛾)(1−𝑅₀)<0; thus 𝐷 has eigenvalue with positive real part, otherwise, det𝐷>0; this is contradiction. Hence, 𝑃₀ is unstable for 𝑅₀>1.

4. Endemic Equilibrium and Stability

Equating each equation in system (2.2) to zero and solving this equilibrium equations, system has the unique positive equilibrium 𝑃^∗=(𝑆^∗,𝐼^∗,𝐽^∗,𝑇^∗,𝑅^∗) for 𝑅₀>1, here 𝑆^∗=Λ𝜇⋅1𝑅₀,𝐼^∗=Λ𝑄₁11−𝑅₀,𝐽^∗=𝛼𝑄₃𝑄₂𝑄₃𝐼−𝜎𝛾^∗𝑇^∗=𝛼𝜎𝑄₂𝑄₃𝐼−𝜎𝛾^∗,𝑅^∗=𝛼𝑘₁𝑄₃+𝑘₂𝜎𝑄₂𝑄₃−𝜎𝛾𝑄₄𝐼^∗. (4.1)

Theorem 4.1. Endemic equilibrium 𝑃^∗ of system is globally asymptotically stable for 𝑅₀>1.

Proof. Equating each equation in system (2.2) to zero, the equilibrium equations as follows are useful: Λ=𝜇𝑆^∗+𝛽₁𝐼^∗+𝛽₂𝐽^∗+𝛽₃𝑇^∗+𝛽₄𝑅^∗𝑆^∗,𝑄₁𝐼^∗=𝛽₁𝐼^∗+𝛽₂𝐽^∗+𝛽₃𝑇^∗+𝛽₄𝑅^∗𝑆^∗,𝑄₂𝐽^∗=𝛼𝐼^∗+𝛾𝑇^∗,𝑄₃𝑇^∗=𝜎𝐽^∗,𝑄₄𝑅^∗=𝑘₁𝐽^∗+𝑘₂𝑇^∗, (4.2) where 𝑄_𝑖 is defined in (3.6).
Setting 𝑥=(𝑆,𝐼,𝐽,𝑇,𝑅)∈Ω⊂ℝ⁺₅, construct a Lyapunov function 𝑈₂=𝑈₂(𝑥)=𝑆−𝑆^∗−𝑆^∗𝑆ln𝑆^∗+𝐴₁𝐼−𝐼^∗−𝐼^∗𝐼ln𝐼^∗+𝐴₂𝐽−𝐽^∗−𝐽^∗𝐽ln𝐽^∗+𝐴₃𝑇−𝑇^∗−𝑇^∗𝑇ln𝑇^∗+𝐴₄𝑅−𝑅^∗−𝑅^∗𝑅ln𝑅^∗, (4.3) where𝑥^∗=𝑃^∗=(𝑆^∗,𝐼^∗,𝐽^∗,𝑇^∗,𝑅^∗) and 𝐴_𝑖>0 is constant. There,𝑈₂(𝑥)≥0 for 𝑥∈IntΩ, and 𝑈₂(𝑥)=0⇔𝑥=𝑥^∗.
Computing the time derivative of 𝑈₂, we have 𝑈^′₂=𝑆^′𝑆1−^∗𝑆+𝐴₁𝐼^′𝐼1−^∗𝐼+𝐴₂𝐽^′𝐽1−^∗𝐽+𝐴₃𝑇^′𝑇1−^∗𝑇+𝐴₄𝑅^′𝑅1−^∗𝑅. (4.4) Using (2.2) and (4.10), we obtain 𝑆^′𝑆1−^∗𝑆=𝑆1−^∗𝑆Λ−𝜇𝑆−𝛽₁𝐼+𝛽₂𝐽+𝛽₃𝑇+𝛽₄𝑅𝑆=𝑆1−^∗𝑆𝜇𝑆^∗+𝛽₁𝐼^∗+𝛽₂𝐽^∗+𝛽₃𝑇^∗+𝛽₄𝑅^∗𝑆^∗−𝜇𝑆−𝛽₁𝐼+𝛽₂𝐽+𝛽₃𝑇+𝛽₄𝑅𝑆=𝜇𝑆^∗𝑆2−𝑆^∗−𝑆^∗𝑆−𝛽₁𝐼+𝛽₂𝐽+𝛽₃𝑇+𝛽₄𝑅𝑆+𝛽₁𝐼+𝛽₂𝐽+𝛽₃𝑇+𝛽₄𝑅𝑆^∗+𝛽₁𝐼^∗+𝛽₂𝐽^∗+𝛽₃𝑇^∗+𝛽₄𝑅^∗𝑆^∗−𝛽₁𝐼^∗+𝛽₂𝐽^∗+𝛽₃𝑇^∗+𝛽₄𝑅^∗𝑆^∗2𝑆. (4.5) Similarly, we obtain 𝐴₁𝐼^′𝐼1−^∗𝐼=𝐴₁𝛽₁𝐼+𝛽₂𝐽+𝛽₃𝑇+𝛽₄𝑅𝑆+𝛽₁𝐼^∗+𝛽₂𝐽^∗+𝛽₃𝑇^∗+𝛽₄𝑅^∗𝑆^∗−𝑄₁𝐼−𝛽₁𝐼+𝛽₂𝐽+𝛽₃𝑇+𝛽₄𝑅𝑆𝐼^∗𝐼,𝐴₂𝐽^′𝐽1−^∗𝐽=𝐴₂𝛼𝐼−𝑄₂𝐽+𝛾𝑇+𝛼𝐼^∗+𝛾𝑇^∗𝐽−𝛼𝐼^∗𝐽𝐽−𝛾𝑇^∗𝐽,𝐴₃𝑇^′𝑇1−^∗𝑇=𝐴₃𝜎𝐽−𝑄₃𝑇+𝜎𝐽^∗𝑇−𝜎𝐽^∗𝑇,𝐴₄𝑅^′𝑅1−^∗𝑅=𝐴₄𝑘₁𝐽+𝑘₂𝑇−𝑄₄𝑅+𝑄₄𝑅^∗+𝑘₁𝐽^∗+𝑘₂𝑇^∗−𝑘₁𝐽𝑅^∗𝑅−𝑘₂𝑇𝑅^∗𝑅. (4.6) Substituting formula (4.5) and (4.6) into (4.4) and arranging the equation we have 𝑈^′₂=𝑄₀+𝑄₁+𝑄₂+𝑄₃+𝑄₄, (4.7) where 𝑄₀=𝜇𝑆^∗𝑆2−𝑆^∗−𝑆^∗𝑆,𝑄₁=𝐴₁−1𝛽₁𝐼+𝛽₂𝐽+𝛽₃𝑇+𝛽₄𝑅𝑄𝑆,₂=𝐴₁+1𝛽₁𝐼^∗+𝛽₂𝐽^∗+𝛽₃𝑇^∗+𝛽₄𝑅^∗𝑆^∗+𝐴₂𝛼𝐼^∗+𝐴₂𝛾𝑇^∗+𝐴₃𝜎𝐽^∗+𝐴₄𝑘₁𝐽^∗+𝐴₄𝑘₂𝑇^∗,𝑄₃=𝛽₁𝐼+𝛽₂𝐽+𝛽₃𝑇+𝛽₄𝑅𝑆^∗+𝐴₂𝛼𝐼+𝐴₂𝛾𝑇+𝐴₃𝜎𝐽+𝐴₄𝑘₁𝐽+𝐴₄𝑘₂𝑇−𝐴₁𝑄₁𝐼−𝐴₂𝑄₂𝐽−𝐴₃𝑄₃𝑇−𝐴₄𝑄₄𝑄𝑅,₄=−𝛽₁𝐼^∗+𝛽₂𝐽^∗+𝛽₃𝑇^∗+𝛽₄𝑅^∗𝑆^∗2𝑆−𝐴₁𝛽₁𝐼+𝛽₂𝐽+𝛽₃𝑇+𝛽₄𝑅𝐼^∗𝐼−𝐴₂𝐽𝛼𝐼^∗𝐽−𝐴₂𝐽𝛾𝑇^∗𝐽−𝐴₃𝑇𝜎𝐽^∗𝑇−𝐴₄𝑘₁𝐽𝑅^∗𝑅−𝐴₄𝑘₂𝑇𝑅^∗𝑅. (4.8) In (4.7), 𝑄₂ consists of all constant terms, 𝑄₃ contains all linear terms of 𝐼, 𝐽, 𝑇, 𝑅, and 𝑄₄ contains all negative nonlinear.
In order to determine the coefficient 𝐴_𝑖 of 𝑈₂, let 𝑄₃≡0 (in Ω); then the coefficients of state variables 𝐼, 𝐽, 𝑇, 𝑅 are equal to zero, that is: 𝛽₁𝑆^∗−𝐴₁𝑄₁+𝐴₂𝛽𝛼=0₂𝑆^∗−𝐴₂𝑄₂+𝐴₃𝜎+𝐴₄𝑘₁𝛽=0₃𝑆^∗+𝐴₂𝛾−𝐴₃𝑄₃+𝐴₄𝑘₂𝛽=0₄𝑆^∗−𝐴₄𝑄₄=0. (4.9) Solving (4.9), and using the expression of 𝑅₀ and 𝑆^∗, we have 𝐴₁=1,𝐴₂=1𝑄₂𝑄₃𝛽−𝜎𝛾²𝑄₃+𝛽₃𝛽𝜎+⁴𝑘₁𝑄₃+𝑘₂𝜎𝑄₄𝑆^∗,𝐴₃=1𝑄₂𝑄₃𝛽−𝜎𝛾²𝛾+𝛽₃𝑄₂+𝛽₄𝑘₁𝛾+𝑘₂𝑄₂𝑄₄𝑆^∗,𝐴₄=𝛽₄𝑄₄𝑆^∗, (4.10) then,𝑄₁=0.
Let𝑆/𝑆^∗=𝑥, 𝐼/𝐼^∗=𝑦, 𝐽/𝐽^∗=𝑧, 𝑇/𝑇^∗=𝑢, 𝑅/𝑅^∗=𝑣; substituting these expressions into (4.9), and then substituting the changing expression into (4.7) 𝑈^′₂=𝜇𝑆^∗12−𝑥−𝑥+2𝛽₁𝑆^∗𝐼^∗+2𝛽₂𝑆^∗𝐽^∗+2𝛽₃𝑆^∗𝑇^∗+2𝛽₄𝑆^∗𝑅^∗+𝐴₂𝛼𝐼^∗+𝐴₂𝛾𝑇^∗+𝐴₃𝜎𝐽^∗+𝐴₄𝑘₁𝐽^∗+𝐴₄𝑘₂𝑇^∗−𝛽₁𝑆^∗𝐼^∗1𝑥−𝛽₂𝑆^∗𝐽^∗1𝑥−𝛽₃𝑆^∗𝑇^∗1𝑥−𝛽₄𝑆^∗𝑅^∗1𝑥−𝛽₁𝑆^∗𝐼^∗𝑥−𝛽₂𝑆^∗𝐽^∗𝑥𝑧𝑦−𝛽₃𝑆^∗𝑇^∗𝑥𝑢𝑦−𝛽₄𝑆^∗𝑅^∗𝑥𝑣𝑦−𝐴₂𝛼𝐼^∗𝑦𝑧−𝐴₂𝛾𝑇^∗𝑢𝑧−𝐴₃𝜎𝐽^∗𝑧𝑢−𝐴₄𝑘₁𝐽^∗𝑧𝑣−𝐴₄𝑘₂𝑇^∗𝑢𝑣=𝜇𝑆^∗12−𝑥−𝑥+𝛽₁𝑆^∗𝐼^∗12−𝑥−𝑥+𝛽₂𝑆^∗𝐽^∗2−𝑥𝑧𝑦−1𝑥+𝛽₃𝑆^∗𝑇^∗2−𝑥𝑢𝑦−1𝑥+𝛽₄𝑆^∗𝑅^∗2−𝑥𝑣𝑦−1𝑥+𝐴₂𝛼𝐼^∗𝑦1−𝑧+𝐴₂𝛾𝑇^∗𝑢1−𝑧+𝐴₃𝜎𝐽^∗𝑧1−𝑢+𝐴₄𝑘₁𝐽^∗𝑧1−𝑣+𝐴₄𝑘₂𝑇^∗𝑢1−𝑣. (4.11) Using the arithmetic mean geometric to get that 2−𝑥−1/𝑥 is less than or equal to zero, substituting 𝐴_𝑖, and simplifying the other expressions in (4.11) after tedious algebraic manipulations, we can get the other expressions such as (2−𝑥𝑧/𝑦−1/𝑥), (2−𝑥𝑢/𝑦−1/𝑥), and (2−𝑥𝑣/𝑦−1/𝑥) is less than or equal to zero; this indicates that 𝑈^′₂≤0, and equalities hold if and only if 𝑥=1,𝑦=𝑧=𝑢=𝑣. Furthermore, 𝑆=𝑆^∗, 𝐼/𝐼^∗=𝐽/𝐽^∗=𝑇/𝑇^∗=𝑅/𝑅^∗=𝑎; then substituting 𝑆=𝑆^∗, 𝐼=𝑎𝐼^∗, 𝐽=𝑎𝐽^∗, 𝑇=𝑎𝑇^∗, 𝑅=𝑎𝑅^∗ into the first equation of system (2.2), and in contrast to the first equality of (4.2), we have 𝑎=1.
In conclusion, the limit sets of solutions in Ω are all in Γ={(𝑆,𝐼,𝐽,𝑇,𝑅)∶𝑆=𝑆^∗,𝐼=𝐼^∗,𝐽=𝐽^∗,𝑇=𝑇^∗,𝑅=𝑅^∗}, and the only invariant set in Γ is 𝑃^∗ by the LaSalle’s invariance principle. Thus the solutions of system (2.2) in Ω are limits to the endemic equilibrium 𝑃^∗, and 𝑃^∗ is globally asymptotically stable for 𝑅₀>1.

5. Discussion

This paper is an extended model about the works in [15, 16] by adding treatment and drug resistance in the whole transmission as well as considering the reasons of treatment exiting.

For public health view, to bring HIV/AIDS into control, the prerequisite is reducing the threshold value of basic reproductive number 𝑅₀. If control 𝑅₀<1, the disease can be eliminated from population. R052 RCT study indicated that treatment can prevent in HIV transmission [18]; this sounds that increasing proportion of treated population is helpful to control HIV epidemic overall. In our study, ART is clearly affecting 𝑅₀ in the HIV procession. However, we cannot yet give this positive result based on the formula of 𝑅₀ and 𝜎 in our work. That is because the treatment might also induce drug resistance which neutralizes the effect of treatment. ART might produce a more complicated HIV progress. However, decreasing acquiring drug-resistant rate 𝑘₁, 𝑘₂ and treatment exiting 𝛾 are the feasible measures to reduce 𝑅₀. Improving treatment standard and patients’ compliance are the feasible and effective measures to reduce 𝑘₁, 𝑘₂. Certainly, new effective antiretroviral drugs might be the real determinants.

𝑅₀ is also linked with drug resistance by the transmission rate 𝛽₄of drug resistance individual and removing rate 𝜌₃ from the population 𝑅. When the other parameters keep constant, 𝑅₀ is positive with 𝛽₄ and negative with𝜌₃. The value of 𝑅₀ will increase if more patients enter this kind of population; that means the drug resistance can fuel HIV epidemic. However, the acquiring drug-resistant rate 𝑘₁, 𝑘₂ can be prevented by improving treatment quality, and the transmission coefficient parameter 𝛽₄ can be reduced by decreasing contacts between these patients and other people at public health level. Generally, early finding by routine screening for drug resistance is an important way to find them.

Limitations in our model exist. We ignored the changing drug when treatment failed in practice. A refinement of the model can be done in future. Additionally, we did do simulation with actual data, which are ongoing under further study.

Acknowledgments

This study was supported by Chinese government grants administered under the Twelfth Five-Year Plan (2012ZX10001-002), the Beijing Municipal Commission of Education (KM201010025010), and National Nature Science Foundation (30973981).