# A Stochastic HIV Infection Model with Latent Infection and Antiretroviral Therapy

Authors: Jun Liu; Yan Wang; Luju Liu; Tingting Zhao
Journal: Discrete Dynamics in Nature and Society
 (2018)

---

## Abstract

Recent studies have demonstrated that the latent infection is a major obstacle to the viral elimination in HIV infection process. In this paper, we formulate a stochastic HIV infection model to include both latent infection and combination drug therapies. We derive that the model solution is unique and positive, and the solution is global. By constructing appropriate stochastic Lyapunov functions, the existence of an ergodic stationary distribution is obtained when the critical condition is greater than one. Furthermore, through rigorous analysis and deduction, the extinction of the virus is established under certain conditions. Numerical simulations are performed to show that small intensity of white noises can maintain the existence of a stationary distribution, while large intensity of white noises is beneficial to the extinction of the virus.


---

## Body

## 1. Introduction

More and more mathematical models have been developed to reflect the dynamics mechanism of HIV virus. The classical HIV virus model is based on three-dimensional ordinary differential equations, which contains the target T-cells population, infected T-cells population, and virions [1–3]. With the progress of HIV drug treatment, some researchers have investigated the effects of drug therapies on model behaviors [4–7]. The commonly used highly active antiretroviral therapy (HAART) is combined by reverse transcriptase inhibitors and protease inhibitors. Reverse transcriptase inhibitors can inhibit the activity of reverse transcriptase and prevent the formation of provirus during the virus infection; protease inhibitors can block the virus infection of new target T-cells. Unfortunately, the current drug treatment can not eradicate the virus thoroughly. With the development of cell and molecular biology, it has been discovered that the existence of latent infection is a major obstacle to clear the virus [8, 9].In recent years, mathematical models with the inclusion of latent infected T-cells have been developed to investigate the model behaviors [10–15]. However, all these models are deterministic model, and the effect of stochastic fluctuation factor is not considered. Actually, HIV transcription is an inherent random process and produces strong random fluctuations in the HIV gene products [16–18]. Thus, by experimental data, it has been proved that these random factors seriously affect the evolution of HIV virus during the protease inhibitor therapy [16]. Moreover, Weinberger et al. have demonstrated that stochastic fluctuation could play an important role in delaying HIV transactivation and contributing to latency [19]. Therefore, it is necessary to consider stochastic fluctuation in HIV virus model, which will be more appropriate to reflect the virus infection process.Recently, the stochastic differential equation models of infectious diseases [20–25] have been greatly developed, but the stochastic differential equation model of the virus dynamic model has just been developed. The classic three-dimensional stochastic model is studied in the earlier literatures [26–29], then the virus stochastic model with target T-cell logistic growth and CTL immune response is also developed in the last two years [30, 31]. However, all these models did not consider the latent infection mechanism. Conway and Coombs have formulated a stochastic model of latent infected T-cells [32], but they neither include the target T-cell logistic growth nor analyze the model dynamics theoretically. In this paper, we will formulate an HIV model with latent infection, healthy T-cell logistic growth, antiretroviral therapy, and random perturbations.The organization of this paper is as follows. In Section2, we construct our stochastic model and give the existence and uniqueness of the global positive solution. By constructing suitable stochastic Lyapunov functions, the existence of a unique ergodic stationary distribution is derived in Section 3. In Section 4, the sufficient conditions for the extinction of the virus are obtained, and numerical simulations are also employed to show how the white noises affect the model behaviors in Section 5. We summarize our results and give future work finally.


## 2. Model and Preliminaries

For deterministic model, Wang et al. [14] have studied the following HIV model with latent infection, T-cell logistic growth, and antiretroviral therapy,(1)ddtTt=λ-dTT+rT1-TTmax-1-nrtkVT,ddtLt=η1-nrtkVT-dLL-aL,ddtIt=1-η1-nrtkVT-δI+aL,ddtVt=1-npNδT∗-cV.Here, T(t), L(t), and I(t) represent the concentrations of healthy T-cells, latent infected T-cells, and actively infected T-cells at time t, respectively. V(t) denote the infectious viral particles at time t. Parameter λ is a constant input rate, r is the logistic growth rate, and Tmax is the maximum capacity of the number of healthy T-cells in the blood. k is the infection rate between healthy T-cells and infectious virus, η is the fraction of infections leading to latency, and a is the activated rate. Parameters dT, dL, δ, and c are the death rates of the healthy T-cells, latently infected T-cells, actively infected T-cells, and infectious virions, respectively. N is the burst size of one productively infected T-cell. nrt and np denote the drug efficacy of reverse transcriptase inhibitors and protease inhibitors (0≤nrt<1 and 0≤np<1), respectively.Wang et al. in literature [14] have analyzed the theoretical results of system (1). System (1) always has an uninfected equilibrium E0(T0,0,0,0), where (2)T0=Tmax2rr-dT+r-dT2+4rλTmax,and the basic reproduction number is(3)R0=kNT01-nrt1-npa+1-ηdLca+dL.When R0>1, only one infected equilibrium E¯(T¯,L¯,I¯,V¯) exists. Wang et al. have obtained the following theoretical results of system (1): (I)
IfR0<1, the uninfected equilibrium E0 is globally asymptotically stable.(II)
IfR0>1, then the infected equilibrium E¯ is locally asymptotically stable under certain condition, and the uninfected equilibrium E0 is unstable.In this paper, considering random fluctuation during HIV infection process [16–19], we follow the assumptions of previous literatures [26–31] and assume that the parameters -dT,  -dL,  -δ, and -c are perturbed by the white noise type. In other words, -dT→-dT+σ1B˙1(t),  -dL→-dL+σ2B˙2(t),  -δ→-δ+σ3B˙3(t), and -c→-c+σ4B˙4(t). Thus, (4)-dTdt⟶-dTdt+σ1dB1t,-dLdt⟶-dLdt+σ2dB2t,-δdt⟶-δdt+σ3dB3t,-cdt⟶-cdt+σ4dB4t,where Bi(t) are independent standard Brownian motions with Bi(0)=0 and σi2>0 denote the intensities of the white noise, i=1,2,3,4. Then, system (1) with white noises can be written as(5)dTt=λ-dTT+rT1-TTmax-1-nrtkVTdt+σ1TdB1t,dLt=η1-nrtkVT-dLL-aLdt+σ2LdB2t,dIt=1-η1-nrtkVT-δI+aLdt+σ3IdB3t,dVt=1-npNδI-cVdt+σ4VdB4t.In the following of this paper, we assume(Ω,F,{Ft}t≥0,P) be a complete probability space with a filtration {Ft}t≥0 satisfying the usual conditions (i.e., it is increasing and right continuous while F0 contains all P-null sets), and Bi(t)(i=1,2,3,4) are defined on this complete probability space. We also let R+4={x=(x1,x2,x3,x4)∈R4:xi>0,i=1,2,3,4}.For convenience, we use symbola∧b to represent min{a,b} and use symbol a∨b to represent max{a,b}. We further denote the following notations:(6)k¯≔1-nrtk,N¯≔1-npN.Based on the theory of Mao [33] and the article of Mao et al. [34], we show the existence and uniqueness of the solution of system (5) in the following theorem.Theorem 1.
For any initial valueT0,L0,I0,V0∈R+4 in system (5), there exists a unique positive solution (T(t),L(t),I(t),V(t)) and the solution will be contained in R+4 with probability one, namely, (T(t),L(t),I(t),V(t))∈R+4 for all t≥0 almost surely (a.s.).Proof.
All the coefficients of system (5) are locally Lipschitz continuous, by the theory of stochastic differential equation [33], so we obtain that there exists a unique maximum local solution Tt,Lt,It,Vt on t∈[0,τe), where τe is the explosion time.
We show the solution is global in the following; namely, we should derive thatτe=∞ a.s. Otherwise, it is supposed that there is a finite time so that the solution Tt,Lt,It,Vt can not explode to infinity. Let n0>0 be big enough so that T(0),L(0),I(0) and V(0) are all on the interval 1/n0,n0. For each integer n≥n0, the stopping time is defined as follows: (7)τn=inf⁡t∈0,τe:min⁡Tt,Lt,It,Vt≤1normax⁡Tt,Lt,It,Vt≥n, and we let inf⁡∅=∞ (commonly, ∅= the empty set). It is obvious that τn increases with n→∞. We denote τ∞=limn→∞⁡τn, then τ∞≤τe a.s.
(i) Ifτ∞=∞ is satisfied a.s., then τe=∞ a.s., so (T(t),L(t),I(t),V(t))∈R+4 a.s. for all t≥0.
(ii) Otherwise, we assume that there exists a pair of constantst~>0 and ϵ∈(0,1) such that P{τ∞≤t~}>ϵ. So, there exists an integer n1≥n0 that makes P{τ∞≤t~}≥ϵforalln≥n1 valid.
Define aC2-function W: R+4→R+, (8)WT,L,I,V=T-b-bln⁡Tb+L-1-ln⁡L+I-1-ln⁡I+12N¯V-1-ln⁡V, where b is a positive constant to be determined. Employing Ito^’s formula [33] to W yields (9)dWT,L,I,V=LWT,L,I,Vdt+σ1T-bdB1t+σ2L-1dB2t+σ3I-1dB3t+σ42N¯V-1dB4t, and LW: R+4→R+ is defined by (10)LWT,L,I,V=λ-dTT+rT1-TTmax-dLL-δ2N¯I-c2N¯V-bλT-dT+r1-TTmax-k¯V-ηk¯VTL-dL-a-1-ηk¯VTI-δ+aLI-12N¯N¯δIV-c+12bσ12+12σ22+12σ32+14N¯σ42=-rT2Tmax+r-dT+brTmaxT-dLL-δ2N¯I-bλT-ηk¯VTL-1-ηk¯VTI-aLI+bk¯-c2N¯V+K≤bk¯-c2N¯V+K. Here, we use the notation K=λ+a+dL+δ+c/2N¯+bdT+1/2bσ12+1/2σ22+1/2σ32+1/4N¯σ42, and we choose the constant 0<b≤c/2k¯N¯, such that bk¯-c/2N¯≤0, so(11)LWT,L,I,V≤K.The rest of the proof is similar to Theorem 2.2 in literature [34]. At the end, we deduce that τ∞=∞. This completes the proof.


## 3. Stationary Distribution

Our main concern of this paper is to investigate whether the random fluctuations could affect the model dynamic behavior or not. According to the theory of Hasminskii [35] in the Appendix, we prove the existence of a unique ergodic stationary distribution, which means the survival of the virus population in the future.Define the critical condition(12)R0s=λk¯N¯δa+1-ηdL+1/2σ22λ/T0+1/2σ12dL+a+1/2σ22δ+1/2σ32c+1/2σ42,and denote σ2=σ12∨σ22∨σ32∨σ42.Theorem 2.
WhenR0s>1, system (5) has a unique ergodic stationary distribution π(·).Proof.
According to LemmaA.1 in the Appendix, we divide it into two steps to prove our conclusion. Now, we prove the validity of condition (A1) in Lemma A.1. The diffusion matrix of system (5) is (13)A=σ12T20000σ22L20000σ32I20000σ42V2.Choose R=min(T,L,I,V)∈D¯l⊂R+4{σ12T2,σ22L2,σ32I2,σ42V2}; we have (14)∑i,j=14aijT,L,I,Vξiξj=σ12T2ξ12+σ22L2ξ22+σ32I2ξ32+σ42V2ξ42≥σ2ξ2,T,L,I,V∈D¯s,ξ=ξ1,ξ2,ξ3,ξ4∈R+4, where D¯l=[1/l,l]×[1/l,l]×[1/l,l]×[1/l,l] and l is a sufficiently large integer, then condition (A1) in Lemma A.1 is satisfied.
In the following, we illustrate the establishment of condition (A2) in LemmaA.1 in the Appendix. Denote(15)fT=λ+r-dTT-rTmaxT2, and f is a quadratic function with respect to T. Suppose T0 and -T1 are the two roots of f(T)=0, then(16)fT=-rTmaxT-T0T+T1. Since -(T-T0)2=-(T-T0)(T+T1)+(T-T0)(T0+T1)≤0, then (T-T0)(T+T1)≥(T-T0)(T0+T1), and(17)fT=-rTmaxT-T0T+T1≤-rTmaxT-T0T0+T1.
We define aC2-function W1: R+4→R,(18)W1=TT0+T1-ln⁡T, and then,(19)LW1=1T0+T1fT-kVT-λT+dT-r1-TTmax+k¯V+12σ12≤fTT0+T1-λT+dT-r+rTTmax+k¯V+12σ12≤-rTmaxT-T0-λT+dT-r+rTTmax+k¯V+12σ12≤-λT+λT0+k¯V+12σ12. Here, we use the equality dT-r+rT0/Tmax=λ/T0.
We construct aC2-function U1: R+4→R,(20)U1=-ln⁡I+c1+d1TT0+T1-ln⁡T-c2+d2ln⁡V-d3ln⁡L, where c1, c2, d1, d2, and d3 are positive constants which will be confirmed later. Applying Ito^’s formula [33] to U1 and using inequality (19) yield (21)LU1≤-1-ηk¯TVI+c1λT+c2N¯δIV-aLI+d1λT+d2N¯δIV+d3ηk¯VTL+c1+d1λT0+12σ12+k¯V+c2+d2c+12σ42+d3dL+a+12σ22+δ+12σ32≤-3c1c21-ηλk¯N¯δ3+c1λT0+12σ12+c2c+12σ42-4d1d2d3aληk¯N¯δ4+d1λT0+12σ12+d2c+12σ42+d3dL+a+12σ22+δ+12σ32+c1+d1k¯V. Let(22)c1λT0+12σ12=c2c+12σ42=1-ηλk¯N¯δλ/T0+1/2σ12c+1/2σ42,d1λT0+12σ12=d2c+12σ42=d3dL+a+12σ22=aληk¯N¯δλ/T0+1/2σ12dL+a+1/2σ22c+1/2σ42.We calculate that (23)c1=1-ηλk¯N¯δλ/T0+1/2σ122c+1/2σ42,c2=1-ηλk¯N¯δλ/T0+1/2σ12c+1/2σ422,d1=aληk¯N¯δλ/T0+1/2σ122dL+a+1/2σ22c+1/2σ42,d2=aληk¯N¯δλ/T0+1/2σ12dL+a+1/2σ22c+1/2σ422,d3=aληk¯N¯δλ/T0+1/2σ12dL+a+1/2σ222c+1/2σ42. Hence,(24)LU1≤-1-ηλk¯N¯δλ/T0+1/2σ12c+1/2σ42-aληk¯N¯δλ/T0+1/2σ12dL+a+1/2σ22c+1/2σ42+δ+12σ32+c1+d1k¯V=-λk¯N¯δa+1-ηdL+1/2σ22λ/T0+1/2σ12dL+a+1/2σ22c+1/2σ42+δ+12σ32+c1+d1k¯V=-δ+12σ32R0s-1+c1+d1k¯V≔-H+BV, where R0s is defined in (12) and(25)H=δ+12σ32R0s-1,B=c1+d1k¯.
Define aC2-function U: R+4→R, in the following form: (26)UT,L,I,V=MU1-ln⁡T-ln⁡L-ln⁡I+11+θT+L+I+V2N¯1+θ≔MU1+U2+U3+U4+U5, where 0<θ<min1,2dL∧δ∧2c/2σ2. Select an appropriate constant M>0 to meet the following requirement: (27)-MH+C≤-2, where (28)C=supT,L,I,V∈R+4-r2TmaxT2+θ-dL2L1+θ-δ4I1+θ-c2V2N¯1+θ+rTTmax+E+dT+dL+a+δ+12σ12+12σ22+12σ32<∞. It is obvious to get that(29)liminfn→∞,T,L,I,V∈R+4∖UnUT,L,I,V=+∞, where Un=(1/n,n)×(1/n,n)×(1/n,n)×(1/n,n) and U(T,L,I,V) is a continuous function. So, U(T,L,I,V) should have a minimum point (T¯0,L¯0,I¯0,V¯0) in the interior of R+4. Therefore, we formulate the following nonnegative C2-function U¯:R+4→R+:(30)U¯T,L,I,V=UT,L,I,V-UT¯0,L¯0,I¯0,V¯0. According to Ito^ formula,(31)LU2=-λT+dT-r1-TTmax+k¯V+12σ12,LU3=-ηk¯TVL+δL+a+12σ22,LU4=-1-ηk¯TVI+δ-aLI+12σ32,LU5=T+L+I+V2N¯θλ-dTT+rT1-TTmax-dLL-δI2-cV2N¯+θ2T+L+I+V2N¯θ-1σ12T2+σ22L2+σ32I2+σ42V2N¯2≤T+L+I+V2N¯θλ+rT-rT2Tmax-dLL-δI2-cV2N¯+θ2σ12∨σ22∨σ32∨σ42T+L+I+V2N¯θ-1T2+L2+I2+V2N¯2≤λ+rTT+L+I+V2N¯θ-rTmaxT2+θ-dLL1+θ-δ2I1+θ-cV2N¯1+θ+θ2σ2T1+θ+L1+θ+I1+θ+V2N¯1+θ≤-r2TmaxT2+θ-dL2L1+θ-δ4I1+θ-c2V2N¯1+θ+E,where (32)E=supT,L,I,V∈R+4-r2TmaxT2+θ-dL2L1+θ-δ4I1+θ-c2V2N¯1+θ+θ2σ2T1+θ+L1+θ+I1+θ+V2N¯1+θ+λ+rTT+L+I+V2N¯θ<∞. According to inequalities (24) and (31), we obtain that (33)LU¯≤-MH+MB+k¯V-λT+rTTmax-ηk¯TVL-1-ηk¯VTI-r2TmaxT2+θ-dL2L1+θ-δ4I1+θ-c2V2N¯1+θ+E+dT+dL+a+δ+12σ12+12σ22+12σ32.
Next, we define a compact subsetDε which satisfies condition (A2) in Lemma A.1. Define the following bounded closed set:(34)Dε=ε≤T≤1ε,ε3≤L≤1ε3,ε3≤I≤1ε3,ε≤V≤1ε, where ε is a small enough constant. In set R+4∖Dε, we further select ε sufficiently small which satisfies the following conditions:(35)-MH+MB+k¯ε+C≤-1,(36)-λε+F≤-1,(37)-ηk¯ε+F≤-1,(38)-1-ηk¯ε+F≤-1,(39)-c42Nε1+θ+F≤-1,(40)-r4ε2+θTmax+F≤-1,(41)-dL4ε31+θ+F≤-1,(42)-δ8ε31+θ+F≤-1,where (43)F=supT,L,I,V∈R+4-r4TmaxT2+θ-dL4L1+θ-δ8I1+θ-c4V2N¯1+θ+MB+k¯V+rTTmax+E+dT+dL+a+δ+12σ12+12σ22+12σ32<∞. Set R+4∖Dε is divided into the following eight domains: (44)D1=T,L,I,V∈R+4,0<V<ε,D2=T,L,I,V∈R+4,0<T<ε,D3=T,L,I,V∈R+4,0<L<ε3,T≥ε,V≥ε,D4=T,L,I,V∈R+4,0<I<ε3,T≥ε,V≥ε,D5=T,L,I,V∈R+4,V>1ε,D6=T,L,I,V∈R+4,T>1ε,D7=T,L,I,V∈R+4,L>1ε3,D8=T,L,I,V∈R+4,I>1ε3. Obviously, Dεc=D1∪D2∪D3∪D4∪D5∪D6∪D7∪D8.
Case  1. If (T,L,I,V)∈D1,(45)LU¯≤-MH+MB+k¯V+rTTmax-r2TmaxT2+θ-dL2L1+θ-δ4I1+θ-c2V2N¯1+θ+E+dT+dL+a+δ+12σ12+12σ22+12σ32≤-MH+MB+k¯ε+C. Based on condition (35), it indicates that LU¯≤-1 for any (T,L,I,V)∈D1.
Case  2. If (T,L,I,V)∈D2,(46)LU¯≤-λT+MB+k¯V+rTTmax-r4TmaxT2+θ-dL4L1+θ-δ8I1+θ-c4V2N¯1+θ+E+dT+dL+a+δ+12σ12+12σ22+12σ32≤-λT+F≤-λε+F. In view of (36), we get LU¯≤-1 for any (T,L,I,V)∈D2.
Case  3. If (T,L,I,V)∈D3,(47)LU¯≤-ηk¯TVL+MB+k¯V+rTTmax-r4TmaxT2+θ-dL4L1+θ-δ8I1+θ-c4V2N¯1+θ+E+dT+dL+a+δ+12σ12+12σ22+12σ32≤-ηk¯TVL+F≤-ηk¯ε+F. According to (37), we obtain that LU¯≤-1 for any (T,L,I,V)∈D3.
Case  4. If (T,L,I,V)∈D4,(48)LU¯≤-1-ηk¯VTI+MB+k¯V+rTTmax-r4TmaxT2+θ-dL4L1+θ-δ8I1+θ-c4V2N¯1+θ+E+dT+dL+a+δ+12σ12+12σ22+12σ32≤-1-ηk¯VTI+F≤-1-ηk¯ε+F. From condition (38), we deduce that LU¯≤-1 for any (T,L,I,V)∈D4.
Case  5. If (T,L,I,V)∈D5,(49)LU¯≤-c2V2N¯1+θ+MB+k¯V+rTTmax-r4TmaxT2+θ-dL4L1+θ-δ8I1+θ+E+dT+dL+a+δ+12σ12+12σ22+12σ32≤-c4V2N1+θ+F≤-c42Nε1+θ+F. In view of (39), we obtain that LU¯≤-1 for any (T,L,I,V)∈D5.
Case  6. If (T,L,I,V)∈D6,(50)LU¯≤-r2TmaxT2+θ+MB+k¯V+rTTmax-dL4L1+θ-δ8I1+θ-c4V2N¯1+θ+E+dT+dL+a+δ+12σ12+12σ22+12σ32≤-r4TmaxT2+θ+F≤-r4ε2+θTmax+F. By condition (40), we conclude that LU¯≤-1 for all (T,L,I,V)∈D6.
Case  7. If (T,L,I,V)∈D7,(51)LU¯≤-dL2L1+θ+MB+k¯V+rTTmax-r4TmaxT2+θ-δ8I1+θ-c4V2N¯1+θ+E+dT+dL+a+δ+12σ12+12σ22+12σ32≤-dL4L1+θ+F≤-dL4ε31+θ+F. It derives that LU¯≤-1 for any (T,L,I,V)∈D7 if condition (41) is satisfied.
Case  8. If (T,L,I,V)∈D8,(52)LU¯≤-δ4I1+θ+MB+k¯V+rTTmax-r4TmaxT2+θ-dL4L1+θ-c4V2N¯1+θ+E+dT+dL+a+δ+12σ12+12σ22+12σ32≤-δ8I1+θ+F≤-δ8ε31+θ+F. By condition (42), it is obvious to get that LU¯≤-1 for any (T,L,I,V)∈D8. Clearly, from (45)–(52), there exists a sufficiently small ε, such that(53)LU¯≤-1∀T,L,I,V∈Dεc. By Lemma A.1 in the Appendix, we obtain that system (5) is ergodic and has a unique stationary distribution. This completes the proof.Remark 3.
For system (5), if there is no white noise, that is, σ1=σ2=σ3=σ4=0, then the critical condition R0s is the basic reproduction number R0 of its corresponding deterministic model (1). It indicates that the existence of the stationary distribution of our stochastic model is a generalization of its corresponding deterministic model to the stability of the infected equilibrium.


## 4. Extinction

Theorem 4.
For any initial value(T(0),L(0),I(0),V(0))∈R+4, one assumes Tt,Lt,It,Vt be the solution of system (5). If r-dT>σ12/2, then for almost ω∈Ω, the solution Tt,Lt,It,Vt of system (5) is(54)limsupt→∞1tln⁡R0dL+aLt+R0aIt+R0aN¯Vt≤ϕ,a.s., where ϕ=min{dL+a,δ,c}(R0-1)IR0≤1+max{dL+a,δ,c}(R0-1)IR0>1+cR0/T0∫0∞x-T0μ(x)dx-1/2(σ2-2+σ3-2+σ4-2)-1, and the distribution of T(t) converges weakly to a measure which has the invariant density(55)μx=Qx-2+2r-dT/σ12exp⁡-2σ12λx+rxTmax,x∈0,∞, where Q is a constant such that ∫0∞μ(x)dx=1. If ϕ<0, then L, I and V will be eliminated with probability one; namely,(56)limt→∞Lt=0,limt→∞It=0,limt→∞Vt=0,a.s.Proof.
Firstly, we can obtain that the solution of system (5) is positive for any given positive initial value from Theorem 1. Consider the following auxiliary stochastic differential equation:(57)dxt=λ-dTxt+rxt1-xtTmaxdt+σ1xtdB1t.For the initial value x(0)=T(0)>0 of system (57), denote x(t) as its solution. Following the proof process in Theorem 3.1 of article [31], by the law of large numbers theorem [36], it is easy to get that x(t) has the ergodic property with the invariant density(58)μx=Qx-2+2r-dT/σ12exp⁡-2σ12λx+rxTmax,x∈0,∞, where Q is a constant such that ∫0∞μ(x)dx=1. It then follows that (59)limt→∞1t∫0txsds=∫0∞xμxdxa.s. By the comparison theorem of one-dimensional stochastic differential equation [37], we obtain that(60)Tt≤xt∀t≥0a.s.
Next, we denote a matrix(61)F=00ηk¯T0dL+aaδ01-ηk¯T0δ0N¯δc0. It is easy to obtain that there must be a left eigenvector v1,v2,v3 of matrix F corresponding to the value R0 [38], such that(62)v1,v2,v3F=R0v1,v2,v3, where (v1,v2,v3)=R0,δR0/a,cR0/aN¯.
Define aC2-function G: R+4→R,(63)GL,I,V=α1L+α2I+α3V, where α1=v1/(dL+a),α2=v2/δ,α2=v3/c. By applying Ito^’s formula, it yields(64)dln⁡G=Lln⁡Gdt+1Gα1σ2LdB2t+α2σ3IdB3t+α3σ4VdB4t,where (65)Lln⁡G=α1Gηk¯TV-dLL-a+α2G1-ηk¯TV-δI+aL+α3GN¯δI-cV-α12σ22L22G2-α22σ32I22G2-α32σ42V22G2≔J-α12σ22L22G2-α22σ32I22G2-α32σ42V22G2. We have(66)G2=α1σ2L1σ2+α2σ3I1σ3+α3σ4V1σ42≤α12σ22L2+α22σ32I2+α32σ42V21σ22+1σ32+1σ42. Moreover, by using (60), we get(67)J=α1Gηk¯TV-dLL-a+α2G1-ηk¯TV-δI+aL+α3GN¯δI-cV≤α1Gηk¯x-T0V+α1Gηk¯T0V-dLL-a+α2G1-ηk¯x-T0V+α2G1-ηk¯T0V-δI+aL+α3GN¯δI-cV≤k¯α3α1η+α21-ηx-T0+1Gv1dL+aηk¯T0V-dLL-a+v2δ1-ηk¯xV-δI+aL+v3cN¯δI-cV=k¯α3α1η+α21-ηx-T0+1Gv1,v2,v3FL,I,VT-L,I,VT=k¯α3α1η+α21-ηx-T0+1GR0-1v1L+v2I+v3V=k¯α3α1η+α21-ηx-T0+1GR0-1α1dL+aL+δα2I+cα3V≤k¯α3α1η+α21-ηx-T0+min⁡dL+a,δ,cR0-1IR0≤1+max⁡dL+a,δ,cR0-1IR0>1. By (66) and (67), we obtain that (68)Lln⁡G≤mindL+a,δ,cR0-1IR0≤1+max⁡dL+a,δ,cR0-1IR0>1+k¯α3α1η+α21-ηx-T0-12σ2-2+σ3-2+σ4-2-1. From (64), we further derive that(69)dln⁡G≤min⁡dL+a,δ,cR0-1IR0≤1+max⁡dL+a,δ,cR0-1IR0>1+k¯α3α1η+α21-ηx-T0-12σ2-2+σ3-2+σ4-2-1dt+1Gα1σ2LdB2t+α2σ3IdB3t+α3σ4VdB4t. Integrate (69) from 0 to t, and divide by t on both sides, then(70)ln⁡Gtt≤ln⁡G0t+min⁡dL+a,δ,cR0-1IR0≤1+max⁡dL+a,δ,c×R0-1IR0>1+k¯α1η+α21-ηα3t∫0txs-T0ds-12σ2-2+σ3-2+σ4-2-1+1t∫0tα1σ2LsGsdB2s+1t∫0tα2σ3IsGsdB3s+1t∫0tα3σ4VsGsdB4s=ln⁡G0t+min⁡dL+a,δ,cR0-1IR0≤1+max⁡dL+a,δ,c×R0-1IR0>1+k¯α1η+α21-ηα3t∫0txs-T0ds-12σ2-2+σ3-2+σ4-2-1+M1tt+M2tt+M3tt, where (71)M1t=∫0tα1σ2LsGsdB2s,M2t=∫0tα2σ3IsGsdB3s,M3t=∫0tα3σ4VsGsdB4s.M1(t),  M2(t), and M3(t) are local martingales whose quadratic variations are 〈M1,M1〉t=σ22∫0tα1L(s)/G(s)ds≤σ22t,  〈M2,M2〉t=σ32∫0tα2I(s)/G(s)ds≤σ32t, and 〈M3,M3〉t=σ42∫0tα3V(s)/G(s)ds≤σ42t. According to the strong law of large numbers for local martingale [33], we can obtain that(72)limt→∞⁡Mitt=0a.s.,i=1,2,3.Because x(t) is ergodic and ∫0∞xμ(x)dx<∞,(73)limt→∞1t∫0txs-T0ds=∫0∞x-T0μxdx<∞. Taking the superior limit on both sides of (70) and combining with (72) and (73) yield (74)limsupt→∞ln⁡Gtt≤min⁡dL+a,δ,cR0-1IR0≤1+max⁡dL+a,δ,c×R0-1IR0>1+k¯α1η+α21-ηα3limt→∞1t∫0txs-T0ds-12σ2-2+σ3-2+σ4-2-1=min⁡dL+a,δ,cR0-1IR0≤1+maxdL+a,δ,c×R0-1IR0>1+cR0T0∫0∞x-T0μxdx-12σ2-2+σ3-2+σ4-2-1≔ϕa.s., which is the required assertion. In particular, if ϕ<0, then it derives that (75)limt→∞ln⁡Ltt=0,limt→∞ln⁡Itt=0,limt→∞ln⁡Vtt=0a.s.; therefore, limt→∞⁡L(t)=0,  limt→∞I(t)=0, and limt→∞V(t)=0 a.s. In other words, the populations L, I, and V will tend to zero exponentially with probability one. This completes the proof.


## 5. Numerical Simulations

The numerical simulation is performed to investigate the effect of the white noises on model dynamical behavior. Applying Milstein’s higher order method [39], system (5) can be rewritten as the following discrete form: (76)Tm+1=Tm+λ-dTm-rTm1-TTmax-k¯VmTmΔt+σ1TmΔtφ1m+0.5σ12TmΔtφ1m2-Δt,Lm+1=Lm+ηk¯VmTm-dLLm+aLmΔt+σ2LmΔtφ2m+0.5σ22LmΔtφ2m2-Δt,Im+1=Im+1-ηk¯VmTm-δIm+aLmΔt+σ3ImΔtφ3m+0.5σ32ImΔtφ3m2-Δt,Vm+1=Vm+N¯δTm∗-cVmΔt+σ4VmΔtφ4m+0.5σ42VmΔtφ4m2-Δt,where φ1m, φ2m, φ3m, and φ4m,  m=1,2,…,n, are the mth realization of four independent Gaussian random variables with distribution N(0,1), and the time step Δt=0.01.In order to examine the effect of the white noises on model behavior numerically, we choose the same parameter values as the deterministic model (1). For the source of the parameter values, see literature [14] for detail. (77)λ=10μl-1day-1,dT=0.03day-1,r=0.1day-1Tmax=1500μl-1,k=0.0001μl-1day-1,η=0.02,dL=0.001day-1,a=0.1day-1,δ=1day-1,N=1000virions/cell,c=20day-1,nrt=0.4,np=0.5,and the white noises σ1=σ2=σ3=0.08,  σ4=0.4. For this set of parameters, we have the critical condition R0s=1.2726>1, and system (5) has a unique stationary distribution (Theorem 2 is satisfied). For its corresponding deterministic model (1), the basic reproduction number R0=1.7657>1, and the unique infected equilibrium is locally asymptotically stable (result (II) in Section 2 is satisfied). When R0s>1, Figure 1 shows that the unique stationary distribution occurs (see Figures 1(a) and 1(c)), and the corresponding distributions of total T-cells population (the sum of the healthy T-cell T, the latently infected T-cells L, and the actively infected T-cells I) and virus population are also given (see Figures 1(b) and 1(d)).Figure 1
The solution of stochastic model (5) and its histogram with white noises σ1=σ2=σ3=0.08 and σ4=0.4. The black lines are the solutions of total T-cells and virus populations of stochastic model (5), and the red lines are the solutions of their deterministic model, respectively.
(a) (b) (c) (d)We further choose the white noisesσ1=σ2=σ3=0.3,  σ4=1.5, and all the other parameters are the same as the above. We calculate that R0s=0.2524<1, and Theorem 2 does not hold. In Figure 2, we can see that the virus can be eradicated for the stochastic model (5), while the virus is prevalent for the deterministic model (1). From Figures 1 and 2, we examine that the existence of stationary distribution can be maintained if the white noises are small, and the virus can be eradicated if the white noises are big enough.Figure 2
The virus can be eradicated for the stochastic model (5) when the white noises σ1=σ2=σ3=0.3 and σ4=1.5, while the virus is prevalent for its corresponding deterministic model.
(a) (b)


## 6. Conclusions

We have studied a stochastic HIV infection model with latent infection, target T-cells logistic growth, and antiretroviral therapy. The critical conditionR0s of the stochastic model which is related to its corresponding deterministic model is derived. We obtain the existence of the unique stationary distribution when the critical condition is bigger than one. Moreover, the sufficient condition which is dependent on the basic reproduction number is obtained for the extinction of the virus. Numerically, we investigate that small fluctuations can be helpful to the persistence of the virus, while large fluctuations can be beneficial to the extinction of the virus.For the special case, if the target T-cells logistic growth is not considered in our stochastic model, then system (5) can be rewritten as the following:(78)dTt=λ-dTT-1-nrtkVTdt+σ1TdB1t,dLt=η1-nrtkVT-dLL-aLdt+σ2LdB2t,dIt=1-η1-nrtkVT-δI+aLdt+σ3IdB3t,dVt=1-npNδI-cVdt+σ4VdB4t.Model (78) is the corresponding stochastic model of the deterministic HIV basic latent infection model in literature [10–12], and our theoretical analysis results can also be established for this model. When there are no drug therapies in system (5), that is, nrt=np=0, then our theoretical results are also satisfied.In the process of virus infection, cytotoxic T lymphocyte (CTL) is the main immune cell in-host, which can limit the replication of virus in the body and determine the viral load [2, 40]. Wang et al. in article [31] have considered the CTL immune response in HIV stochastic model; therefore, a stochastic HIV model with latent infection and CTL immune response will be our future work. In addition, under the influence of environment, the stochastic perturbation is not always treated as a constant but a periodic behavior [41]; thus HIV latent infection model regarding periodic perturbation is also our research work in the future.



---

*Source: 5175383-2018-06-25.xml*