MATHEMATICAL MODELLING OF INTERNAL HIV DYNAMICS

. We study a mathematical model for the viral dynamics of HIV in an infected individual in the presence of HAART. The paper starts with a literature review and then formulates the basic mathematical model. An expression for R 0 , the basic reproduction number of the virus under steady state application of HAART, is derived followed by an equilibrium and stability analysis. There is always a disease-free equilibrium (DFE) which is globally asymptotically stable for R 0 < 1. Deterministic simulations with realistic parameter values give additional insight into the model behaviour. We then look at a stochastic version of this model and calculate the proba- bility of extinction of the virus near the DFE if initially there are only a small number of infected cells and infective virus particles. If R 0 ≤ 1 then the system will always approach the DFE, whereas if R 0 > 1 then some simulations will die out whereas others will not. Stochastic simulations suggest that if R 0 > 1 those which do not die out approach a stochastic quasi-equilibrium consisting of random ﬂuctuations about the non-trivial deterministic equilibrium levels, but the amplitude of these ﬂuctuations is so small that practically the system is at the non-trivial equilibrium. A brief discussion concludes the paper.


MATHEMATICAL MODELLING OF INTERNAL HIV DYNAMICS
Nirav dalal, David greenhalgh and Xuerong mao Department of Statistics and Modelling Science Livingstone Tower, 26 Richmond Street, Glasgow G1 1XH, U.K.
Abstract. We study a mathematical model for the viral dynamics of HIV in an infected individual in the presence of HAART. The paper starts with a literature review and then formulates the basic mathematical model. An expression for R 0 , the basic reproduction number of the virus under steady state application of HAART, is derived followed by an equilibrium and stability analysis. There is always a disease-free equilibrium (DFE) which is globally asymptotically stable for R 0 < 1. Deterministic simulations with realistic parameter values give additional insight into the model behaviour.
We then look at a stochastic version of this model and calculate the probability of extinction of the virus near the DFE if initially there are only a small number of infected cells and infective virus particles. If R 0 ≤ 1 then the system will always approach the DFE, whereas if R 0 > 1 then some simulations will die out whereas others will not. Stochastic simulations suggest that if R 0 > 1 those which do not die out approach a stochastic quasi-equilibrium consisting of random fluctuations about the non-trivial deterministic equilibrium levels, but the amplitude of these fluctuations is so small that practically the system is at the non-trivial equilibrium. A brief discussion concludes the paper.
benefits of using these models is that it helps us in identifying important parameters and factors which have dominant effect on the development, transmission and spread of the disease. This gives a better understanding for development of treatment strategies and we can draw biologically relevant interpretations.
In this paper we shall study a mathematical model for HIV internal viral dynamics. HIV viral dynamics is defined as the interaction of the human immune system with HIV. Infection with HIV causes depletion of CD4 T cells, also known as the helper T cells. CD4 T cells play a central role in our immune system. They are involved in protecting against viral, fungal and protozoal infections. These cells normally orchestrate the immune response, signalling other cells in the immune system to perform their special functions [18]. The depletion of CD4 T cells plays a pivotal role in the degeneration of the immune system in HIV infected individuals.
Two types of HIV infect humans: HIV-1 and HIV-2. HIV-1 is the more virulent and easily transmitted, and is the source of the majority of HIV infections throughout the world. Various methods have been suggested to effectively tackle HIV-1. Early detection and treatment of the disease is one of them [5]. Treatment with a single drug is usually not effective in containing the virus. However individuals infected with HIV can be treated with a combination of antiretroviral drugs. This is known as highly active antiretroviral therapy (HAART). HAART works by reducing the HIV viral load in the plasma (usually measured in the blood) to an undetectable level. This delays the onset of the debilitating and ultimately fatal symptoms of AIDS for a very long time, perhaps indefinitely. Nowadays HAART is usually effective but it is expensive, can involve a complicated medication regime including many drugs and has unpleasant side effects. Moreover as HIV mutates rapidly resistance is an increasing problem. Treatment with HAART greatly affects the HIV internal viral dynamics. Viral dynamics under the influence of HAART has been a topic of some major research work recently [2], [4], [6] and [29].
In the past mathematical models have been an important tool in modelling viral dynamics. Several models have already been developed to describe the internal viral dynamics of HIV-1 infection within an individual. Bonhoeffer et al. [5] started with a basic model and went on to include the development of drug resistant virus. Work has also been done on characterising HIV-1 dynamics incorporating virus resistance [21]. Kamina et al. [14] have worked on a stochastic model for early HIV-1 population dynamics based on a multi-dimensional diffusion process. Tuckwell and Le Corfec [26] use a system of stochastic differential equations to model virus dynamics. The model was formed to model the initial stages of infection. The aim was to model the intrinsic variability in HIV-1 growth and to explore effects of perturbation in the parameter values. Wick and Self [31] used a branching process to model the early events in HIV-1 infection and study the influence that the time of appearance of virus specific antibodies or the administration of antiretroviral drugs has on probability of progression to chronic infection. One set of stochastic models have looked at the effects of increasing variability among viral strains, as a means of escaping control by the immune system in progression to AIDS [19,20]. Vergu et al. [17] showed the impact of viral diversity on the immune response and disease dynamics. Nelson et al. [17] included less than perfect drug effects and a delay in the initiation of virus production. A system of delay differential equations was used to model the system. It is widely accepted that HAART leads to maintenance of low levels of viral load. It is important to identify the factors which lead to depletion of HIV-1 virus under HAART. Wu et al. [32] described the reasons for the decay of HIV-1 virus in terms of decay rates. They observed that in persons with moderately advanced HIV-1 infection treated with HAART, there is significant inter-subject variability in the decay rates of HIV-1 virus. These differences were not related to reported or measurable differences in adherence or drug levels. Stafford et al. [25] suggested that cytotoxic T lymphocyte (CTL) destruction of productively infected cells and suppression by CD8+ T cell antiviral factor (CAF) are responsible for lowering the viral load.
The main aim of this paper is to improve existing mathematical models to describe internal HIV-1 viral dynamics in the presence of HAART treatment. We start by formulating the model which includes drug treatment effects into the basic model developed by Bonhoeffer et al. [5] and identify the possible equilibria. The basic reproductive number R 0 is a key parameter of the model. We perform a global stability analysis and supplement our theoretical results by numerical simulations of the model. We look at stochastic effects and the probability of extinction in the model. Again our theoretical results are supplemented by numerical simulations. A brief discussion concludes the paper.
2. Model formation. The basic model of Bonhoeffer et al. [5] which describes HIV-1 viral dynamics is given by the following set of differential equations: Here x represents the number of uninfected cells, y the number of infected cells and v the number of virus particles at time t. λ is the total rate of production of new uninfected cells per unit time, d, a and u the per capita death rates of uninfected cells, infected cells and virus particles respectively. Each infected cell produces free virus particles at rate k per unit time. Infected cells are produced from uninfected cells and free virus particles at total rate βxv per unit time. Bonhoeffer et al. [5] defined the basic reproductive number R 0 = βλk adu to be the expected number of secondary virus particles produced by a single new virus particle entering a population of uninfected cells at equilibrium with neither infected cells nor virus present. For the above model two equilibria existed. The disease free equilibrium (DFE) was x =x = λ d , y =ŷ = 0, v =v = 0 and the endemic equilibrium was given by They went on to include the dynamics of resistant strains into the basic model. The authors were of the opinion that treatment should start as early as possible and with as many drugs as clinically possible. They also estimated the length of the treatment required to eliminate or minimise the viral load. Ways were suggested to reduce the chance of emergence of multiple drug resistant virus.
Nowak et al. [21] also included strains of resistant virus in their model. Analytical approximation for the rate of emergence of resistant virus was studied in this paper. The decline of wild-type virus and the rise of resistant mutant virus in different compartments of the virus populations such as free plasma virus, cells infected with actively replicating virus, long lived infected cells and cells carrying defective provirus was studied. The aim of the paper was to determine the rate of production of plasma virus by cells, half life of infected cells and half life of free plasma virus.
Culshaw and Ruan [9] obtained a restriction on the number of viral particles released per infectious cell in order for infection to be sustained. They also obtained sufficient conditions on the parameters for the stability of the infected steady state. After introducing a time delay into the model which described the time between infection of a CD4 T cell and the emission of viral particles on a cellular level, they derived that the same restriction on the number of viral particles released per infectious cell is required. Stability conditions for the infected steady state in terms of the parameters and independent of the delay were obtained.
Mittler et al. [16] analysed a model for the interaction of HIV-1 virus with target cells and included a time delay between initial infection and the formation of productively infected cells. They were of the opinion that it was theoretically possible to estimate not only the viral clearance rate and the death rate of productively infected cells, but also the mean and the variance of the infection delay. They obtained an analytical solution for a particular case of delay in which the infection delay conforms to a gamma distribution. They also showed using simulated datasets that the non-linear least-square regression methods used to analyse clinical data can, in principle, extract viral dynamic parameters from models in which distributed delays are present.
Huang et al. [12] developed a viral dynamic model to evaluate how time-varying drug exposure and drug susceptibility affected antiviral response. Plasma concentrations were modelled using a standard pharmacokinetic one-compartment open model with first order absorption and elimination as a function of fixed individual pharmacokinetic parameters and dose times. They observed that poor adherence may result in early viral rebound and a higher set point after treatment failure. Particular patterns of non-adherence affect responses differentially. Also longer sequences of missed doses increase the chance of treatment failure and accelerate the failure.
There does not seem to be much work which provides analytical conditions for decay in viral load levels. We wish to provide some insight about the parameters and conditions which would have major ramifications on the viral load levels. Tuckwell and Wan [27] analyse a very similar model to ours. However, they do not consider HAART in their model. In addition, we provide a global stability proof and a probabilistic approach for this model.
We propose the following three dimensional model to describe the viral dynamics in the presence of HIV-1 infection and HAART: with suitable initial conditions.
The above model also includes the effect of HAART. HAART is generally a combination of reverse transcriptase inhibitor (RTI) drugs and protease inhibitor (PI) drugs. RTI drugs are designed to prevent the conversion of HIV RNA to DNA in early stages of HIV replication. Thus RTI drugs block conversion of uninfected cells to infected cells. PI drugs are designed to intervene in the last stage of the virus replication cycle to prevent HIV from being properly assembled, and thus cause the newly produced virus to be noninfectious [10]. The treatment effects are accordingly incorporated in the above model.
The variables and parameters in the model are described as follows: x(t) is the concentration of uninfected cells; y(t) is the concentration of infected cells; v(t) is the concentration of virus particles; (1-γ) is the reverse transcriptase inhibitor drug effect; (1-η) is the protease inhibitor drug effect; λ is the total rate of production of healthy cells per unit time; δ is the per capita death rate of healthy cells; β is the transmission coefficient between uninfected cells and infective virus particles; a is the per capita death rate of infected cells; N is the average number of infective virus particles produced by an infected cell in the absence of HAART during its entire infectious lifetime; and u is the per capita death rate of infective virus particles.
Note that in the absence of HAART, assuming that the lifetime of an infected cell has an exponential distribution each infected cell is infectious for time (1/a). During its total infectious lifetime it produces (N a)/a = N virus particles.
Also observe that when a single infective virus particle infects a single uninfected cell the virus particle is absorbed into the uninfected cell and effectively dies. Hence the term (1−γ)βxv appears in equation (3) as well as equations (1) and (2). Perelson and Nelson [23] ignored the term (1 − γ)βxv in equation (3). The arguments given were first that in examining data from patients with different T-cell counts there was no statistically significant correlation between the T-cell count and the clearance rate. However this may not always be true for all patients in all situations. The second argument given is that if x is approximately constant then the clearance rate can be redefined. As we wish to examine a dynamic model where x may vary it is better to have this term in the model. Also given that there is a biological justification for this term there is no harm done by including it as the model will be more accurate with this term in than not. One of the models studied by Tuckwell and Wan [27] also includes this term as does a stochastic model by Tuckwell and Le Corfec [26]. However Tuckwell and Le Corfec's work is relevant only to the early stages of HIV infection.
3. Equilibrium and stability results. We begin our analysis by determining the equilibria that exist and the conditions for their existence. We then move on to investigating the conditions under which those equilibria are possible. This would give us an insight into the long-term behaviour of our model.

For our model the equilibrium with no infected cells and no infected virus particles is
Consider a single infected virus particle entering this equilibrium. Such a particle is infectious for time (u + (1 − γ)β λ δ ) −1 and during this time infects δu+βλ(1−γ) uninfected cells which each survive for average time 1 a and during this time produce infective virus particles at rate (1−η)N a. Hence the basic reproductive number, defined as the number of secondary infective virus particles caused by a single infective virus particle entering this equilibrium is Note also that a similar argument shows that R 0 can be interpreted as the expected number of secondary cells infected by a single infected cell entering the DFE.
Theorem 3.1. The model (1)-(3) describing the viral dynamics has a unique DFE which is globally asymptotically stable for R 0 < 1.
Proof. Clearly our model has a unique DFE x, We divide the proof into three parts. We first show that limsup x(t) ≤ λ/δ, then we prove that y(t), v(t) both tend to zero, then that x(t) → λ/δ.
Hence we can write Next we show that y and v tend to zero as t → ∞.
Now choose ǫ 1 small enough so that, We can do this because the right hand side of (4) tends to R 0 < 1 as ǫ 1 → 0. Note that with this value of ǫ 1 There exists t 0 such that t ≥ t 0 implies that If π(0) > 0 then it is straightforward to show that y(t) and v(t) are strictly positive for t > 0. If π(0) > 0 then we can write where ξ 3 = min ξ1 ψ , ξ 2 > 0. This implies that On the other hand if π(0) = 0 then it is straightforward that y(t) and v(t) → 0 as t → 0.
Tuckwell and Wan [27] analysed a functionally similar model, but without HAART, corresponding to γ = η = 0 in our model. If in our model we define new parameters k = (1 − γ)β, c = (1 − η)N a, s = λ, µ = δ andγ = u and substitute the new parameters into our model we get exactly Tuckwell and Wan's model. Hence by using this reparameterisation we can apply their analytical results directly to our model. Doing this we deduce that the threshold condition ks(c − a)/µaγ < 1 in Tuckwell and Wan's parameters becomes (1 − γ)βλ((1 − η)N − 1)/δu < 1 in our parameters (i.e. R 0 < 1). Applying Tuckwell and Wan's results to our model we deduce that for R 0 < 1 there is only the DFE which is locally asymptotically stable (LAS). For δu < (1 − γ)βλ((1 − η)N − 1) (R 0 > 1) the DFE is unstable and there is a unique non-trivial equilibrium with infected cells and virus particles present which is LAS. In brief summary for our model for R 0 ≤ 1 there is only the DFE which is LAS for R 0 < 1. There is a unique equilibrium with infected cells and virus particles present if and only if R 0 > 1 and this equilibrium is LAS when it exists. Moreover the DFE is unstable for R 0 > 1.
Theorem 3.1 extends the results of [27] in that it shows global not local stability of the DFE for R 0 < 1 using completely different mathematical methods than [27]. It also incorporates the effect of HAART and identifies R 0 as the basic reproduction number of the virus (alternatively the invasion of infected cells) into the DFE. 4. Deterministic simulations. As we have seen our analytical results suggest the existence of two equilibria for the system (1) -(3). An equilibrium with only uninfected cells present and no infected cells or virus particles always exists and is globally asymptotically stable for R 0 ≤ 1 and is unstable for R 0 > 1. A nontrivial equilibrium with uninfected cells, infected cells and virus particles exists for R 0 > 1 and is locally stable when it exists. We now use simulations with realistic parameter values to explore the likely behaviour of the model further. The numerical integration package SOLVER was used to integrate the differential equations. We required two sets of parameter values for the simulations. One set was to be used to get a value of R 0 < 1 and the other one for R 0 > 1. Hence the parameter values were sourced from different publications so that they would be realistic but need not be significantly altered from the published values.
Our parameter values have all been taken from published literature. λ, β, δ, u, a have all been taken from [5] and N from [6]. The parameter values for As can be seen from Figure 1, for R 0 ≤ 1 the system converges to the equilibrium with only uninfected cells present. The infected cells and the infective virus particles die out and the number of uninfected cells converge to λ/δ. For realistic parameter values even with large initial numbers of cells and virus particles the number of infected cells dropped to 10% of its initial value in around six days. The number of virus particles rises sharply in the beginning but quickly peaks and drops to 10% of its value in approximately nine days. The level of both was virtually zero after only ten days. The number of susceptible cells took slightly longer to reach its very large equilibrium value of 1 × 10 8 susceptible cells, but had reached 90% of this value by roughly 25 days and the equilibrium level had virtually been reached after 45 days. The simulations were repeated for a large number of different starting values and in each case the system approached the equilibrium with only uninfected cells present.
For Figure 2 the parameter values λ, a, u have been taken from [17] and β, δ are taken from [21]  Deterministic HIV viral dynamics model in vivo approaches the disease-free equilibrium for R 0 < 1 of infected cells and virus particles which causes a corresponding sharp decrease in the number of susceptible cells. After this the infected cells and virus particles decline and then tend to their equilibrium values. This equilibrium is reached very quickly after about only 60 days. This simulation was repeated with a variety of different initial conditions and in each case x → x * , y → y * and v → v * as t → ∞ suggesting that provided that at least one uninfected cell or one infective virus particle is initially present the system will tend to this equilibrium. In [27] it is shown that for the model without HAART that this equilibrium is LAS when it exists. Our results strengthen these analytical results in two ways: First our simulations suggest that the equilibrium with infected cells and virus particles is actually globally stable when it exists in the sense that if R 0 > 1 and at least one infected cell or virus particle is initially present then the system of cells and virus particles will tend to this equilibrium.
Second the model has been extended to include the effect of HAART. Usually initially in an infection individuals are not treated with HAART but the extension to include HAART is still useful, as some individuals such as healthcare workers exposed to needlestick injuries will be on HAART from the beginning. Moreover many other individuals will be diagnosed later and start HAART then. The model is then applicable to the course of their internal HIV viral dynamics both before HAART is started, and with suitably modified parameters, after HAART is started.  Deterministic HIV viral dynamics model in vivo approaches the endemic equilibrium for R 0 > 1 5. Probabilistic approach. While deterministic models have their uses a different way of looking at the problem of HIV internal viral dynamics would be to use a stochastic model. This can tell us about quantities such as the probability of extinction which deterministic models cannot. A straightforward stochastic analogue of the deterministic model is to assume that the individual stochastic events occur according to a Poisson process with the corresponding rate being the same as in the deterministic model. Further details are given below. Then we can use standard branching process theory to calculate the probability of extinction of the virus. Wick and Self [31] also use a branching process to model the HIV virus dynamics. The difference between their and our work lies in the assumption of 'no competition or cooperation' between the cells and virus particles in Wick and Self. We do not make any such assumption. Also their results are mainly numerical.
When an individual is initially infected with HIV, the usual situation is that they will not be treated with HAART. However there are cases where an exposed individual will be treated with HAART, for example healthcare workers exposed to needlestick injuries [7,13]. Hence we give the theoretical results with HAART included and then state them without HAART. In the numerical examples we suppose that HAART is not given, as this is the more usual situation.
We are interested in estimating the probability that the virus invades the DFE. To calculate this we shall assume that the number of susceptible cells is kept fixed at the deterministic equilibrium level λ/δ. This approximation will be valid whilst the number of virus particles and infected cells remains small. This approximation will be valid near the DFE and is a similar approximation to that used to calculate the probability of a major epidemic outbreak in Whittle's simple stochastic epidemic model [30]. In that paper it is assumed that the number of susceptibles remains constant and the number of infectives is a stochastic birth and death process in order to calculate the probability of a major epidemic outbreak. Once the stochastic process moves away from the DFE this approximation will no longer be valid but it can be used to estimate the probability that the virus and infected cells invade the DFE.
The stochastic model assumes that the probability that each single infected cell dies in [t, t + δt] is aδt + o(δt). The probability that each single virus particle dies in [t, t+δt] is uδt+o(δt). Each virus particle independently infects each susceptible cell in [t, t+ δt] with probability (1− γ)βδt+ o(δt). When this happens the virus particle is absorbed into the infected cell. Also each infected cell produces a virus particle in [t, t + δt] with probability (1 − η)N aδt + o(δt). All infected cells, virus particles, death processes, infections and production of new virus particles are independent.
The translation of the deterministic model into the stochastic model is not a unique process. As virus particles are produced continuously by infected cells we have chosen to assume that each infected cell produces a virus particle in [t, t + δt] with probability (1 − η)N aδt + o(δt), and otherwise produces no virus particles in [t, t + δt]. One possible alternative assumption might be to assume that each infected cell simultaneously produces N virus particles in [t, t + δt] with probability (1 − η)aδt + o(δt) and otherwise produces no virus particles in [t, t + δt]. This alternative would lead to a qualitatively different solution.
We consider the stochastic model at the equilibrium consisting only of uninfected cells x = λ/δ, y = 0 and v = 0 and consider an initial small infective population ofã infective virus particles andb infected cells entering this equilibrium, whereã andb are small. Define q 1 = (1−γ)βλ (1−γ)βλ+uδ and q 2 = (1−η)N (1−η)N +1 . q 1 is the probability that a single infective virus particle at this equilibrium infects a cell before the virus particle ceases to be infective. q 2 is the probability that a single infected cell at the above equilibrium creates one infective virus particle before it dies. Note that Theorem 5.1. Consider a small initial infective population ofã infective virus particles andb infective cells entering the equilibrium consisting only of uninfected cells. The probability that the population of infective virus particles and infective cells will go extinct without creating a major outbreak is 1 for R 0 ≤ 1 and is for R 0 > 1.
Proof. Suppose first that one infective virus particle enters the equilibrium consisting only of uninfected cells. This particle either dies out at a rate u or ceases to be infective by infecting an uninfected cell at rate (1− γ)β λ δ . Let X denote the number of initially uninfected cells infected by the virus particle before it either loses its infectivity or dies. Then Recall that when a single infective virus particle infects a single uninfected cell the virus particle is absorbed into the uninfected cell and effectively dies. Hence there are only two possible cases that can happen. The probability generating function (p.g.f.) of X is given by Now suppose that a single infected cell enters the equilibrium consisting only of uninfected cells. Let Y be the number of infected virus particles created before it dies. Then for k = 0, 1, 2, 3, . . .
Now let Z 1 denote the number of secondary infective virus particles infected directly by one initial infective virus particle entering the equilibrium consisting only of uninfected cells. The p.g.f. of Z 1 is given by where X 1 is the number of cells infected by this infective virus particle. The probability of extinction is given by the smallest root of s = G Z1 (s) in [0,1], [8]. Hence we get (s − 1)((1 − q 1 q 2 ) − q 2 s) = 0. We obtain the two roots of the above equation as 1 and 1−q1q2 q2 . If R 0 ≤ 1, 1−q1q2 q2 ≥ 1 and the probability of extinction is 1, whereas if R 0 > 1, 1−q1q2 q2 < 1 and the probability of extinction is 1−q1q2 q2 . Similarly let Z 2 be the number of secondary infective cells produced by a single infective cell entering the DFE.
Here Y 1 is the number of infective virus particles produced by this infective cell. 13 Again solving s = G Z2 (s) we get Hence if R 0 ≤ 1, 1−q2 q1q2 ≥ 1 and the probability of extinction is 1, whereas if R 0 > 1, 1−q2 q1q2 < 1 and the probability of extinction is 1−q2 q1q2 .
So the probability of ultimate extinction of a small population ofã infective virus particles andb infected cells, entering the DFE is 1 if R 0 ≤ 1 and The proof is thus complete.
In the more practically useful case where HAART is not present the same formula holds where q 1 = βλ βλ+δu and q 2 = N N +1 .
6. Stochastic simulations. Usually newly infected individuals will not be treated with HAART, but occasionally they will be. Hence we keep HAART in theoretical calculations but omit them in the following numerical examples and simulations.
We have looked at the simulation of the stochastic system with the same set of parameter values as the deterministic system. The results look very similar.
According to the analytical results in the previous section the probability of extinction ofã initial virus particles andb initial infected cells is 1 for R 0 ≤ 1 and is for R 0 > 1. To confirm these observations we carried out stochastic simulations. The simulations assume that the parameter values are the same fixed values as in the deterministic simulations but that events occur at random times given by an exponential distribution with the appropriate rate as described in the detailed description of the stochastic model given in Section 5. Our simulations program was written in FORTRAN and comprehensively verified using detailed output from a large number of runs.
The results obtained were as expected. For those simulations where the number of infected cells and virus particles died out (including those with R 0 ≤ 1) the stochastic effects were clearly visible on the trajectories of the number of virus particles and the infected cells, but although stochastic effects were present in the trajectory of uninfected cells, the scaling used meant that this trajectory looked much smoother and more like the deterministic trajectory. Figure 3 shows such a simulation with parameters and initial values as in Figure 1 except γ = η = 0. For R 0 > 1 we performed simulations with parameters and initial values as in Figure  2 except γ = η = 0. Some of these simulations died out and some did not. Those simulations which did not die out tended to the deterministic endemic equilibrium values. Although stochastic effects could be seen if these trajectories were plotted on a very fine scale those trajectories which took off looked practically identical to the deterministic ones in Figure 2, if the vertical axis scaling for Figure 2 was used.
For two sets of parameter values we also calculated the probabilities that the simulations died out starting with one infective virus particle entering the stochastic   Table 1. Calculated extinction probabilities, 95% confidence intervals and observed number of 1000 simulations which go extinct for two parameter sets DFE, that is y = 0, v = 1, and one infected cell entering the stochastic DFE, that is y = 1, v = 0. One set was the values used for Figure 2. The only change was the change in treatment effects, that is γ = η = 0. For the second set of values λ, u were taken from [17], β from [21], δ from [12] and a from [6]. N remained the same. The values were β = 1 × 10 −8 day −1 dm 3 , λ = 10 6 day −1 dm −3 , γ = η = 0, a = 0.7 day −1 , δ = 0.072 day −1 , u = 3 day −1 and N = 100. The theoretically calculated probabilities of extinction were checked numerically by running one thousand independent simulations for each starting value. A 95% confidence interval for the number of the thousand simulations to go extinct was calculated using the Normal approximation to the Binomial. In both cases the actual number of simulations which went to extinction lay in the confidence interval. The numbers are given in Table 1. 7. Discussion of results. We have discussed an ordinary differential equation model which describes the internal HIV viral dynamics in the presence of treatment such as antiretroviral drugs. We derived an expression for the basic reproductive number R 0 and found that if R 0 was less than or equal to one then the equilibrium with only uninfected cells was the unique equilibrium and was globally asymptotically stable. This paper gives analytical results for the model (1)- (3). Most other work in this area concentrates on obtaining numerical solutions to the problem under consideration, which depend on parameter values used [27,31]. Instead we have established some concrete mathematical results for our model as these will be true for all possible parameter values.
Some authors have previously neglected the term (1 − γ)βxv in the model in equation (3) as numerically insignificant. This was motivated by statistical studies on viral clearance rates in patients. However this term is biologically realistic as each time that the virus infects a cell it must enter that cell, which means that it cannot infect any more cells. Moreover other authors have also included this term [27,31]. It is of interest to include this term in the model as it is biologically sensible and it may have an effect on the dynamics in some situations and for certain parameter values.
We have also taken a stochastic approach in this paper. The processes of virus replication and the increase of infected cells are highly interdependent processes. Separation of one from the other is not possible. Our branching process approximation reflects this interdependence. These results give us the condition for the eradication of the virus and the infected cells in terms of the parameters of the model.
Some recent research has suggested the presence of latent reservoirs in the body in which HIV-1 persists even after the administration of HAART [15,24,22]. These reservoirs are not biologically significant when the level of virus in the body is high but become significant at low virus loads. They are not important in the deterministic model where viral loads are high. These might possibly be significant when considering the probability of extinction of the virus, however the timescale on which these reservoirs are established is not yet clear.
In summary our paper has extended the results in [27] to include the effect of HAART. We have additionally shown that the threshold value derived in [27] corresponds to the basic reproduction number, R 0 , even when the model includes the effect of HAART. For the deterministic model both with and without HAART for R 0 ≤ 1 there is only the DFE which is LAS if R 0 < 1. For R 0 > 1 the DFE is unstable and there is a unique equilibrium with infected cells and virus particles present which is LAS. Moreover we have shown global stability of the DFE when R 0 < 1 which was not shown in [27]. Then we examined the probability of extinction using a branching process model. We have also obtained numerical results for both the deterministic and the stochastic versions of this model. The deterministic simulations suggest that for the deterministic model (either with or without HAART) if R 0 > 1 and at least one infected cell or virus particle is initially present then the system will tend to the equilibrium with both infected cells and virus particles present.