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 ( [subscript] R0[/subscript] ) by the next generation matrix and then analyzed stability of two equilibriums by constructing Lyapunov function. When [subscript] R0[/subscript] <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.

(ProQuest: ... denotes non-US-ASCII text omitted.)

Qianqian Li 1 and Shengshan Cao 1 and Xiao Chen 1 and Guiquan Sun 2 and Yunxi Liu 3 and Zhongwei Jia 4, 5

Recommended by Youssef Raffoul

1, School of Mathematical Sciences, Ocean University of China, Qingdao 266100, China
2, Department of Mathematics, North University of China, Taiyuan 030051, China
3, Department of Nosocomial Infection Management and Disease Control, Chinese PLA General Hospital, 28 Fuxing Road, Haidian District, Beijing 100853, China
4, National Institute of Drug Dependence, Peking University, Beijing 100191, China
5, Takemi Program in International Health, Department of Global Health and Population, Harvard School of Public Health, 665 Huntington Avenue, Boston, MA 02115, USA

Received 11 August 2012; Accepted 22 October 2012

This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

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 S ( t ) , I ( t ) , J ( t ) , A ( t ) , and T (t ) 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 R (t ) 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 d 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.

[figure omitted; refer to PDF]

The model was collected as the following differential equations: [figure omitted; refer to PDF] where Λ is the recruitment rate of the susceptible population, _{β 1} is the probability of transmission by an infection in the first stage, _{β 2} is the probability of transmission by an infection in the second phase, _{β 3} =l _{β 2} (l <1 ) is the probability of transmission by a patient being treated, _{β 4} is the probability of transmission by a drug resistance individual, α is the transfer rate constant from the asymptomatic phase I to the symptomatic phase J , σ is the proportion of treatment, γ is the exit treatment rate without developing drug resistance, _{k 1} is probability constant of infection by transmission of drug-resistant strains, and _{k 2} is the rate of acquiring drug resistance during treatment; _{ρ i} (i =1,2 ,3 ) denote transfer rate constant by an infection from phase J , T , R to the AIDS cases A , respectively.

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

Since the A of system ( 2.1) does not appear in the equations, in the following analysis, we only consider the system as follows: [figure omitted; refer to PDF]

Theorem 2.1.

Let the initial data be S ( 0 ) = _{S 0} >0 , I ( 0 ) = _{I 0} >0 , J ( 0 ) = _{J 0} >0 , T ( 0 ) = _{T 0} >0 and R ( 0 ) = _{R 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 Ω = { ( S ,I ,J ,T ,R ) ∈ ^{... + 5} :S +I +J +T +R ...4; Λ / μ , 0 <S ...4; Λ / μ } , and Ω for system ( 2.2) is positively invariant.

Proof.

From the first equation of ( 2.2) [figure omitted; refer to PDF] consider the following two categories.

(1) When _{t 0} >0 and S ( _{t 0} ) =0 .

: Equation ( 2.3) becomes ^{S [variant prime]} = Λ ( t = _{t 0} ) , due to Λ >0 ; we have ^{S [variant prime]} ( _{t 0} ) >0 , that is, when t > _{t 0} , S ( t ) is an increasing function about t . Therefore, we can conclude that when t is the neighborhood of _{t 0} and t > _{t 0} , S ( t ) ...5;S ( _{t 0} ) >0 .

(2) When _{t 0} >0 and S ( _{t 0} ) >0 .

: Equation ( 2.3) can be written as [figure omitted; refer to PDF] that is [figure omitted; refer to PDF] thus, [figure omitted; refer to PDF]

Similarly for the other equations of system ( 2.2) we can easily show that I ( t ) , J ( t ) , T ( t ) , R (t ) are increasing functions about t when _{t 0} >0 and I ( _{t 0} ) =0 , _{t 0} >0 and J ( _{t 0} ) =0 , _{t 0} >0 and T ( _{t 0} ) =0 , _{t 0} >0 and R ( _{t 0} ) =0 , respectively; then when t > _{t 0} , X ( t ) ...5;X ( _{t 0} ) >0 ( X =I ,J ,T ,R ) . Otherwise, when _{t 0} >0 and X ( _{t 0} ) >0 ( X =I ,J ,T ,R ) , we have the following results corresponding to the respective hypothesis: [figure omitted; refer to PDF] Thus, [figure omitted; refer to PDF] Thus, [figure omitted; refer to PDF] Thus, [figure omitted; refer to PDF] Thus, [figure omitted; refer to PDF] Then, we have that I ( t ) , J ( t ) , T ( t ) , and R (t ) are all strictly positive for t >0 . Thus we can conclude that all solutions of system ( 2.2) remain positive for all t >0 .

Next, add all the the equations of system ( 2.2); we have [figure omitted; refer to PDF]

Then, [figure omitted; refer to PDF]

In a similar fashion we have ^{S [variant prime]} ...4; Λ - μS from the first equation of ( 2.2); then _{lim t [arrow right] ∞} supS ...4; Λ / μ .

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), _{P 0} = ( Λ / μ ,0,0 ,0,0 ) . Following the paper [ 17], we obtain the basic reproduction number by using the next generation operator approach.

Let y = ^{(I ,J ,T ,R ,S ) T} ; thus we have [figure omitted; refer to PDF] where [figure omitted; refer to PDF] Then the derivatives of ... (y ) and ...B1; (y ) at the DFE _{P 0} - = (0,0 ,0,0 , Λ / μ ) are partitioned as [figure omitted; refer to PDF] F and V are the 4 ×4 matrices as follows: [figure omitted; refer to PDF] where [figure omitted; refer to PDF]

Hence the reproduction number, denoted by _{R 0} , is the spectral radius of the next generation matrix F ^{V -1} : [figure omitted; refer to PDF] Here [figure omitted; refer to PDF]

3.2. DFE and Stability

Theorem 3.1.

The disease-free equilibrium _{P 0} of system is globally asymptotically stable for _{R 0} <1 and unstable for _{R 0} >1 .

Proof.

(1) The Jacobian matrices of system ( 2.2) at the DFE are [figure omitted; refer to PDF] where [figure omitted; refer to PDF]

Obviously, -D = [ _{d ij} ] is a 4 ×4 matrix with _{d ij} ...4;0 , for i ...0;j ,i ,j =1 , ... ,4 and _{d ii} >0 for i =2 , ... ,4 . If _{R 0} <1 , we have _{R 1} = ( _{β 1} / _{Q 1} ) · ( Λ / μ ) < _{R 0} <1 from the expression of ( 3.8); thus _{d 11} = _{Q 1 - β 1} · ( Λ / μ ) = _{Q 1} ( 1 - _{R 1} ) >0 .

Define the positive vector subsequently: [figure omitted; refer to PDF] If _{R 0} <1 , -D ·x = ^{[ _{Q 1} ( 1 - _{R 0} ) ,0,0 ,0 ] T} ...5;0 .

Then begin to show that all the eigenvalues of -D are nonzero: [figure omitted; refer to PDF] Simplify the above expression through substituting formula ( 3.8) into ( 3.12): [figure omitted; refer to PDF] Here _{Q 2 Q 3} - σ γ = ( _{ρ 1} + σ + _{k 1} + μ ) ( _{ρ 2} + γ + _{k 2} + μ ) - σ γ >0 ; thus, det ( -D ) >0 ; namely, -D has non zero eigenvalue for _{R 0} <1 .

In conclusion, if _{R 0} <1 , -D is an irreducible matrix with _{d ii} >0 and _{d ij} ...4;0 (i ...0;j ) , there exists a positive vector x such that -D ·x ...5;0 . Hence, the real part of each nonzero eigenvalue of -D is positive according to the M -matrix theory; that is, each eigenvalue of D has negative real part.

Through the structure of the Jacobian matrix J ( _{P 0} ) , it can be seen that the eigenvalues of J ( _{P 0} ) consist of - μ and all eigenvalues of D . Hence, all eigenvalues of J ( _{P 0} ) have negative real part for _{R 0} <1 ; thus, disease-free equilibrium _{P 0} is locally asymptotically stable.

(2) Let _{U 1} = _{m 1} I + _{m 2} J + _{m 3} T + _{m 4} R , where _{m i} ( i =1 ...4 ) are positive constants as follows: [figure omitted; refer to PDF] When _{R 0} <1 , the time derivative of _{U 1} is [figure omitted; refer to PDF] Let M =max ( _{β 1} / _{m 1} , _{β 2} / _{m 2} , _{β 3} / _{m 3} , _{β 4} / _{m 4} ) ; due to _{lim t [arrow right] ∞} supS ...4; Λ / μ , M is the finite number, then [figure omitted; refer to PDF] Solve ^{U 1 [variant prime]} ...4;M ( _{R 0} -1 ) _{U 1} , have _{U 1} ...4; _{U 1} (0 ) ^{e M ( _{R 0} -1 )t} , that is, lim t [arrow right] ∞ _{U 1} ( t ) =0 , thus when t [arrow right] ∞ , ( I ,J ,T ,R ) [arrow right] (0,0 ,0,0 ) . When _{R 0} <1 , ^{U 1 [variant prime]} ...4;0 , and equalities hold if and only if S = Λ / μ ,I =J =T =R =0 , that is ^{U 1 [variant prime]} =0 if and only if t [arrow right] ∞ . We can conclude that the solutions of system ( 2.2) are all in Ψ = { (S ,I ,J ,T ,R ) :S = Λ / μ ,I =J =T =R =0 } and the only invariant set in Ψ is _{P 0} by the LaSalle's invariance principle. Thus the solutions of system ( 2.2) are limits to the endemic equilibrium _{P 0} when _{R 0} <1 . Combine locally asymptotically stable of _{P 0} with convergence properties of the _{P 0} , we conclude that _{P 0} of system is globally asymptotically stable for _{R 0} <1 .

(3) If _{R 0} >1 , detD = _{Q 1 Q 4} ( _{Q 2 Q 3} - σ γ ) ( 1 - _{R 0} ) <0 ; thus D has eigenvalue with positive real part, otherwise, detD >0 ; this is contradiction. Hence, _{P 0} is unstable for _{R 0} >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 ^{P *} = ( ^{S *} , ^{I *} , ^{J *} , ^{T *} , ^{R *} ) for _{R 0} >1 , here [figure omitted; refer to PDF]

Theorem 4.1.

Endemic equilibrium ^{P *} of system is globally asymptotically stable for _{R 0} >1 .

Proof.

Equating each equation in system ( 2.2) to zero, the equilibrium equations as follows are useful: [figure omitted; refer to PDF] where _{Q i} is defined in ( 3.6).

Setting x = ( S ,I ,J ,T ,R ) ∈ Ω ⊂ ^{... 5 +} , construct a Lyapunov function [figure omitted; refer to PDF] where ^{x *} = ^{P *} = ( ^{S *} , ^{I *} , ^{J *} , ^{T *} , ^{R *} ) and _{A i} >0 is constant. There, _{U 2} (x ) ...5;0 for x ∈Int Ω , and _{U 2} ( x ) =0 ...x = ^{x *} .

Computing the time derivative of _{U 2} , we have [figure omitted; refer to PDF] Using ( 2.2) and ( 4.10), we obtain [figure omitted; refer to PDF] Similarly, we obtain [figure omitted; refer to PDF] Substituting formula ( 4.5) and ( 4.6) into ( 4.4) and arranging the equation we have [figure omitted; refer to PDF] where [figure omitted; refer to PDF] In ( 4.7), _{Q 2} consists of all constant terms, _{Q 3} contains all linear terms of I , J , T , R , and _{Q 4} contains all negative nonlinear.

In order to determine the coefficient _{A i} of _{U 2} , let _{Q 3} ...1;0 (in Ω ); then the coefficients of state variables I , J , T , R are equal to zero, that is: [figure omitted; refer to PDF] Solving ( 4.9), and using the expression of _{R 0} and ^{S *} , we have [figure omitted; refer to PDF] then, _{Q 1} =0 .

Let S / ^{S *} =x , I / ^{I *} =y , J / ^{J *} =z , T / ^{T *} =u , R / ^{R *} =v ; substituting these expressions into ( 4.9), and then substituting the changing expression into ( 4.7) [figure omitted; refer to PDF] Using the arithmetic mean geometric to get that 2 -x - 1 / x is less than or equal to zero, substituting _{A i} , and simplifying the other expressions in ( 4.11) after tedious algebraic manipulations, we can get the other expressions such as ( 2 - xz / y - 1 / x ) , ( 2 - xu / y - 1 / x ) , and ( 2 - xv / y - 1 / x ) is less than or equal to zero; this indicates that ^{U 2 [variant prime]} ...4;0 , and equalities hold if and only if x =1 ,y =z =u =v . Furthermore, ^{S =S *} , I / ^{I *} = J / ^{J *} = T / ^{T *} = R / ^{R *} =a ; then substituting ^{S =S *} , I =a ^{I *} , J =a ^{J *} , T =a ^{T *} , R =a ^{R *} into the first equation of system ( 2.2), and in contrast to the first equality of ( 4.2), we have a =1 .

In conclusion, the limit sets of solutions in Ω are all in Γ = { ( S ,I ,J ,T ,R ) : ^{S =S *} , ^{I =I *} , ^{J =J *} , ^{T =T *} , ^{R =R *} } , and the only invariant set in Γ is ^{P *} by the LaSalle's invariance principle. Thus the solutions of system ( 2.2) in Ω are limits to the endemic equilibrium ^{P *} , and ^{P *} is globally asymptotically stable for _{R 0} >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 _{R 0} . If control _{R 0} <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 _{R 0} in the HIV procession. However, we cannot yet give this positive result based on the formula of _{R 0} 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 _{k 1} , _{k 2} and treatment exiting γ are the feasible measures to reduce _{R 0} . Improving treatment standard and patients' compliance are the feasible and effective measures to reduce _{k 1} , _{k 2} . Certainly, new effective antiretroviral drugs might be the real determinants.

_{R 0} is also linked with drug resistance by the transmission rate _{β 4} of drug resistance individual and removing rate _{ρ 3} from the population R . When the other parameters keep constant, _{R 0} is positive with _{β 4} and negative with _{ρ 3} . The value of _{R 0} 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 _{k 1} , _{k 2} can be prevented by improving treatment quality, and the transmission coefficient parameter _{β 4} 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).

References

[1] WHO, "UNAIDS World AIDS Day Report," 2011

[2] Ministry of Health, "2011 Estimates for the HIV/AIDS Epidemic in China," 2011

[3] R. M. Anderson, "The role of mathematical models in the study of HIV transmission and the epidemiology of AIDS," Journal of Acquired Immune Deficiency Syndromes , vol. 1, no. 3, pp. 241-256, 1988.

[4] R. M. Anderson, G. F. Medley, R. M. May, A. M. Johnson, "A preliminary study of the transmission dynamics of the human immunodeficiency virus (HIV), the causative agent of AIDS," IMA Journal of Mathematics Applied in Medicine and Biology , vol. 3, no. 4, pp. 229-263, 1986.

[5] R. M. May, R. M. Anderson, "Transmission dynamics of HIV infection," Nature , vol. 326, no. 6109, pp. 137-142, 1987.

[6] C. C. McCluskey, "A model of HIV/AIDS with staged progression and amelioration," Mathematical Biosciences , vol. 181, no. 1, pp. 1-16, 2003.

[7] H. Guo, M. Y. Li, "Global dynamics of a staged-progression model for HIV/AIDS with amelioration," Nonlinear Analysis: Real World Applications , vol. 12, no. 5, pp. 2529-2540, 2011.

[8] Y. Wang, Y.-c. Zhou, "Mathematical modeling and dynamics of HIV progression and treatment," Chinese Journal of Engineering Mathematics , vol. 27, no. 3, pp. 534-548, 2010.

[9] S. Blower, "Calculating the consequences: HAART and risky sex," AIDS , vol. 15, no. 10, pp. 1309-1310, 2001.

[10] M. C. Boily, F. I. Bastos, K. Desai, B. Mâsse, "Changes in the transmission dynamics of the HIV epidemic after the wide-scale use of antiretroviral therapy could explain increases in sexually transmitted infections--results from mathematical models," Sexually Transmitted Diseases , vol. 31, no. 2, pp. 100-112, 2004.

[11] M. Bachar, A. Dorfmayr, "HIV treatment models with time delay," Comptes Rendus--Biologies , vol. 327, no. 11, pp. 983-994, 2004.

[12] S. M. Blower, H. B. Gershengorn, R. M. Grant, "A tale of two futures: HIV and antiretroviral therapy in San Francisco," Science , vol. 287, no. 5453, pp. 650-654, 2000.

[13] O. Sharomi, A. B. Gumel, "Dynamical analysis of a multi-strain model of HIV in the presence of anti-retroviral drugs," Journal of Biological Dynamics , vol. 2, no. 3, pp. 323-345, 2008.

[14] N. J. D. Nagelkerke, P. Jha, S. J. De Vlas, E. L. Korenromp, S. Moses, J. F. Blanchard, F. A. Plummer, "Modelling HIV/AIDS epidemics in Botswana and India: impact of interventions to prevent transmission," Bulletin of the World Health Organization , vol. 80, no. 2, pp. 89-96, 2002.

[15] L. Cai, X. Li, M. Ghosh, B. Guo, "Stability analysis of an HIV/AIDS epidemic model with treatment," Journal of Computational and Applied Mathematics , vol. 229, no. 1, pp. 313-323, 2009.

[16] J. Q. Li, Y. L. Yang, W. Wang, "Qualitative analysis of an HIV transmission model with therapy," Chinese Journal of Engineering Mathematics , vol. 26, no. 2, pp. 226-232, 2009.

[17] P. van den Driessche, J. Watmough, "Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission," Mathematical Biosciences , vol. 180, pp. 29-48, 2002.

[18] M. S. Cohen, Y. Q. Chen, M. McCauley, "Prevention of HIV-1 infection with early antiretroviral therapy," The New England Journal of Medicine , vol. 365, no. 6, pp. 493-405, 2011.

[]