A stochastic model of the dynamics of HIV under a combination therapeutic intervention

Drug resistance to single therapeutic treatment in HIV infected individuals has promoted research into combined treatments. In this paper we propose a stochastic model under combined therapeutic treatment by extending the model of HIV pathogenesis under treatment by anti-viral drugs given in [Perelson AS, Neumann AU, Markowits M, Leonard JM & Ho DD, 1996, HIV-1 dynamics in vivo virion clearance rate, infected cell life span, and viral generation time, Science New Series, 271, pp. 1582–1586]. Variance and co-variance structures of variables are obtainable via this approach in addition to the mean numbers of free HIV, infectious free HIV and non-infectious free HIV that were obtained by Perelson et al. Comparing simulated data for before and after treatment indicates the importance of combined treatment and its overall effect(s).


Introduction
In HIV infected individuals, the infection exhibits a long asymptomatic phase (after the initial high infectious phase) of approximately 10 years on average before the onset of AIDS.During this incubation period which some call the clinical latency period, the individuals appear to be well and may contribute significantly to the spread of the epidemic in a community.Some clinical markers such as the CD4 cell count and the RNA viral load (viraemia) provide information about the progression of the disease in infected individuals.Also, the clinical latency period of the disease may provide a sufficiently long period during which to attempt an effective suppressive therapeutic intervention in HIV infections.
The knowledge of principal mechanisms of viral pathogenesis, namely the binding of the retrovirus to the gp120 protein on the CD4 cell, the entry of the viral RNA into the target cell, the reverse transaction of viral RNA to viral DNA, the integration of the viral DNA with that of the host, the viral regulatory processes mediated through regulatory proteins such as tat and rev, and the action of viral protease in cleaving viral proteins into mature products have led to the design of drugs (chemotherapeutic agents) to control the production of HIV.Principal directions along which drugs (such as AZT and Ritonavir (see Shafer, et al. (2001)) are currently attempted include inhibition of the reverse transcriptase of HIV and inhibition of the protease of HIV.
With the spread of the HIV-AIDS pandemic and in the absence of an "effective" vaccine or cure, therapeutic interventions are still heavily relied on.Several research studies have been carried out in the recent past, both theoretically and experimentally, to analyse the impact of therapy on the viral load in HIV infected persons in order to ascertain the effectiveness of the treatment (see, for example, Nelson & Perelson (1995), Wei, et al. (1995), Perelson, et al. (1996), Mellors, et al. (1997), Nijhuis, et al. (1998), Tan and Xiang (1999) and Bangsberg, et al. (2004)).Nelson and Perelson (1995) proposed a mathematical model of therapeutic intervention to delay the onset of AIDS by the stimulated production of genetically engineered defective interfering virus (DIVs) that interferes with the HIV replication process.A DIV is a deletion mutant and it is incapable of replicating by itself in a host cell (CD4 cell), but may replicate if the host cell is co-infected with HIV.Assuming that DIV depends on HIV to multiply, Nelson and Perelson (1995) constructed a mathematical model describing the interaction among HIV, DIV and uninfected CD4 cells, and they analysed the co-evolution of DIV and HIV in a single compartment.Their model is essentially a system of ordinary differential equations involving eight variables and several parameters representing the activities of DIV and HIV.By considering a higher level of DIV activity in the production of co-infected CD4 cells, they investigated the possibility of blocking the production of HIV so that the burden of HIV on the population of CD4 cells is reduced.
In the paper of Wei, et al. (1995), based upon the results of several experimental studies of the dynamics of HIV replication in the presence of antiretroviral agents, it was reported that HIV had enormous potential in showing resistance to drugs by undergoing several mutations and a rapid and virtually complete replacement of wild-type HIV by a drug resistant virus occurred when anti-viral drugs were administered.Nijhuis, et al. (1998) noticed high drug resistance and unique combination of mutation in individuals when they proposed a stochastic model to test the resistance to protease inhibitors, although there was a reduced effective free HIV population (500-15 000).Perelson, et al. (1996) presented a mathematical model for analysing the kinetic picture of HIV pathogenesis subject to the administration of a drug, called Ritonavir, to inhibit potently the protease of HIV.In their paper, they represented the dynamics of cell infection and viral production after treatment with Ritonavir by means of a set of ordinary differential equations, adopting a deterministic approach and assuming that the viral inhibition of Ritonavir was 100% so that all newly produced virions after the treatment with Ritonavir were non-infectious.By using the mathematical model and non-linear least squares fitting of the viral load data of five HIV-1 infected patients, they were able to obtain estimates of the rate of viral clearance, the infected cell life-span and the average viral generation time.
Tan and Xiang (1999) developted a state-space model of HIV pathogenesis in HIV infected individuals undergoing combination treatment (i.e. a treatment with a combination of antiviral drugs such as AZT and Ritonavir which can inhibit either the reverse transcriptase or the protease of HIV).Their model included a mathematical description of the production of infectious free HIV and non-infectious free HIV, by extending the model of Perelson, et al. (1996) and developing procedures for estimating and predicting the number of uninfected CD4 cells, infectious free HIV, non-infectious free HIV and HIV infected CD4 cells.They not only extended the model by Perelson, et al. (1996) to a stochastic model, but also applied their model to data of patients considered by Perelson, et al. (1996).Their model is discrete in time and comprises a system of stochastic difference equations which were derived from the biological specifications of the HIV-replication cycle.
Since Perelson, et al. (1996) considered a deterministic model, and Tan and Xiang (1999) developed a state-space model, the present authors consider a stochastic model of the growth of the HIV population which carries over the principle of the virology of HIV, the life-cycle of HIV and allows the production of non-infectious (defective) free HIV to reduce the severity of HIV in a HIV-infected individual undergoing a combinationtherapeutic treatment.Our aim in this paper is to use a stochastic model obtained by extending the model of Perelson, et al. (1996) to determine the number of uninfected T4 cells, infected T4 cells and free HIV in an infected individual by examining the combined antiviral treatment of HIV.This is important because it helps in determining the efficacy of methods used in the research areas of pathogenesis, progression and combined treatment of HIV.Obtaining the variance and co-variance structures of variables representing the number of virus producing cells, and the levels of infectious free and non-infectious free HIV is one of the main contributions of this paper.These variance and co-variance structures may sometimes be difficult to obtain for stochastic models (unlike the case with the mean structure), but being able to obtain these structures is expected to shed light on the type of relationship between the variables.Based on the model, we obtain the expected numbers of HIV infected cells, infectious free HIV and non-infectious free HIV at any time t, and derive conclusions for the reduction or elimination of HIV in HIV-infected individuals.
The organisation of this paper is as follows: In Section 2, we formulate our stochastic model describing the production and the clearance of virus producing cells, infectious free HIV and non-infectious free HIV in a therapeutic environment.In Section 3, we derive a system of differential difference equations for the probability function associated with the process and also obtain a partial differential equation for the probability generating function of the numbers of HIV-infected CD4 cells, infectious free HIV and non-infectious free HIV at time t.The population measures are derived in Section 4. In Section 5, we provide a numerical illustration to demonstrate the impact of using combination therapy to control the progression of HIV, and we also obtain variance and co-variance structures of the variables.Some concluding remarks follow in Section 6.
Assume that at time t = 0, a combination therapy treatment is initiated in an HIV-infected individual.We assume that the therapeutic intervention inhibits either the enzyme action of reverse transcriptase or that of the protease of HIV in a HIV-infected cell.A HIVinfected cell with the inhibited HIV-transcriptase may be considered a dead cell as it cannot participate in the production of the copies of any type of HIV.On the other hand, an HIV-infected cell in which the reverse transcription has already taken place and the viral DNA is fused with the DNA of the host, but the enzyme activity of HIV-protease is inhibited, undergoes a lysis releasing infectious free HIV and non-infectious free HIV.A non-infectious free HIV cannot successfully infect a CD4 cell.Accordingly, at any time t, the blood of the infected person contains virus-producing HIV-infected cells, infectious free HIV and non-infectious free HIV.
A virus producing cell existing at time t in the therapeutic environment undergoes one of the following transitions during the interval (t, t + ∆): 1. it splits into two HIV-infected cells with probability λ 1 ∆ + o(∆); 2. it undergoes a lysis with probability υ∆ + o(∆), producing a random number K 1 of infectious free HIV and a random number K 2 of non-infectious free HIV; 3. it dies with probability µ∆ + o(∆); 4. it remains as it is with probability 1 -( We assume that K 1 and K 2 have the joint probability generating function h(s 1 , s 2 ) defined by where π lm represents the probability that l infections free HIV and m non-infectious free HIV are released at the lysis occurring at any time.An infectious free HIV existing at time t in the blood of the individual may undergo one of the following transitions during the interval (t, t + ∆): 1. it infects a T4 cell with probability λ 2 ∆ + o(∆), converting the cell into a virus producing cell; 2. it dies with probability c∆ + o(∆); 3. it remains as it is with probability 1 -( The population of non-infectious free HIV does not grow by replication of its members, but grows by admitting bulk immigrations which occur at the lysis of HIV-infected cells.A non-infectious HIV existing at time t dies during the interval (t, t + ∆) with probability c∆ + o(∆).
Let X(t) be the number of virus producing cells (these are cells that produce copies of the virus to infect other cells) at time t.Let V (t) and D(t) be respectively the number of infectious free HIV (these are HIV in the body that infect cells in the body) and the number of non-infectious free HIV (these are HIV in the body that do not infect cells in the body) at time t.For simplicity, we assume that X(0) = N , V (0) = n, D(0) = 0. We proceed to discuss the probability generating function of the vector process (X(t), V (t), D(t)) in the next section.

Probability generating function
The probability generating function of (X(t), To derive an equation for G(u 1 , u 2 , u 3 ; t), we require the probability function which is defined for any time t by p(i, j, k; t) = P r{X(t) = i, V (t) = j, D(t) = k}, where i, j, k = 0, 1, 2. . .We proceed to derive a system of differential-difference equations for the function p(i, j, k; t).To this end we list below the exhaustive and mutually exclusive events that occur during the interval (t, t+∆) given that X(t) = i > 0, V (t) = j > 0 and D(t) = k > 0: 1. one HIV infected cell splits into two HIV-infected cells -the probability of this event occuring is iλ 1 ∆ + o(∆); 2. one HIV-infected cell undergoes a lysis -the probability of this event occuring is iυ∆ + o(∆); 3. one HIV-infected cell dies -the probability of this event occuring is iµ∆ + o(∆); 4. one infectious free HIV virus infects one CD4 cell, making the CD4 cell an HIVinfected cell -the probability of this event occuring is jλ 2 ∆ + o(∆); 5. one infectious free HIV virus dies -the probability of this event occuring is jc∆ + o(∆); 6. one non-infectious free HIV virus dies -the probability of this event occuring is kc∆ + o(∆); 7. none of the above occurs.
Now, we have and so, by using ( 2) and ( 3), we obtain Equation ( 4) is not solvable, even for any simple form of h(u 2 , u 3 ; t).However, it is possible to obtain the moment-structure of (X(t), V (t), D(t)).We do this in the next section.
4 The moment-structure of (X(t), V (t), D(t)) We introduce the variables together with the parameters m 1 , m 2 , m 11 , m 22 and m 12 which are respectively the derivatives ∂h(s 1 ,s 2 ) Differentiating (4) once with respect to u 1 , u 2 and u 3 at A, we obtain respectively.Differentiating (4) twice with respect to u 1 , u 2 and u 3 at A, we obtain

∂M
(2) respectively.Finally, differentiating (4) with respect to u 1 and u 2 , with respect to u 2 and u 3 and with respect to u 1 and u 3 at A, we obtain (2) (2) respectively.Writing these equations in the matrix form, we therefore obtain the matrix differential equation where The characteristic equation of the matrix R is given by Solving (15), we obtain the characteristic values of R, which are real and distinct.The corresponding characteristic vectors are respectively, where Accordingly, the general solution of ( 14) is where C 1 , C 2 and C 3 are constants.In our model, we have assumed that X(0) = N , V (0) = n, D(0) = 0 and so we have the following initial conditions: M X (0) = N , M V (0) = n and M D (0) = 0. Consequently, the constants C 1 , C 2 and C 3 satisfy the system of linear equations Solving this system, we obtain , and We have not obtained explicitly the expressions for M (2) D (t), M XV (t), M V D (t) and M XD (t).However, we are able to solve completely the system (5)-(13) in the special case where no infectious free virus is released at the lysis of every HIV-infected cell treated with combination therapy.It holds for this special case that m 1 = 0, m 11 = 0, m 12 = 0. Consequently, M XD (t) = A 1 e −αt + A 2 e −βt + A 3 e −2αt + A 4 e −2βt + A 5 e −(α+c)t + A 6 e −(β+c)t + A 7 e −(α+β)t , and where Although the above expressions for M XD (t) and M D (t) are somewhat complicated, we have presented them here for the sake of completeness.However, for the purpose of numerical illustration considered in the next section, we prefer the integral expressions

Numerical illustration of the model
For the purpose of numerical illustration, we have extrapolated estimates from Perelson, et al. (1996) and Tan and Xiang (1999), and we consider three cases: Case (i): Both the mean numbers of the infectious free HIV (m 1 ) and non-infectious free HIV (m 2 ) produced by a virus producing cell at the time of its lysis are greater than zero.
Case (ii): m 1 = 0 and m 2 = 0 and we obtain values of the means M X (t), M V (t) and M D (t) for the values of t ranging from 0.0 to 2.5 in steps of 0.5 for all the cases.The results are shown in Figure 1.

Case (iii):
The second moments are evaluated by adopting Simpson's one-third rule for the evaluation of integrals in (28)-(33).The assumed values of the parameters are given in Table 1 and the results are shown in Tables 2-6.
For all simulated results, we take 1 hour as 0.5 time unit.

Case (i)
From Figure 1 and Table 2, it is easily noted that as t increases, the values of M X (t), M V (t) and M D (t) also increase for λ 1 = 5.When λ 2 = 10 (the rate at which an HIV infected cell splits into two), the values of M X (t), M V (t) and M D (t) also increase with an increase of t (see Table 3).This shows that as the value of λ 1 increases, the values of M X (t), M V (t) and M D (t) increase significantly as functions of time, before treatment.

Case (ii)
Assume m 1 = 0, m 2 = 100 and m 22 = 9 900. Figure 1 shows the curves for M X (t), M V (t) and M D (t) fitted before and after treatment.It is observed that there is a remarkable difference between the values obtained before and after treatment, especially after t = 1.5.This shows the effectiveness of the treatment.As such, the expected number of virus producing cells and the expected number of non-infectious free HIV decrease significantly after treatment (effect of reverse transcriptase drugs).Furthermore, the expected numbers of infectious free HIV is reduced to almost zero at t = 2.5, because of the effect of protease inhibitor drugs, as they reduce the generation of infectious free HIV at the death of actively infected T4 cells.

Case (iii)
Assume that m 1 = 0, m 2 = 100 and m 22 = 9 900.The values of the second order moments, namely M (2) V (t), M XD (t), M XV (t) and M V D (t), are provided in Tables 5  and 6.The variances of virus producing cells and non-infectious free HIV are so large in comparison to those of infectious free HIV that their values increase significantly as t increases, unlike those of infectious free HIV which decrease significantly after treatment.The co-variance results show that there is a positive relationship between virus producing cells and infectious free HIV.Table 4: M X (t), M V (t) and M D (t) as a function of t (after treatment) with C = 3, N = 412, n = 98 000, υ = 0.02, µ = 0.49, λ 1 = 5, λ 2 = 1, m 2 = 100 and m 22 = 9 900.

Concluding remarks
In this paper we have shown the usefulness of our stochastic approach towards modelling combined HIV treatment by obtaining the variance and co-variance structure of the number of virus producing cells at time t, the number of infectious free HIV and the number of non-infectious free HIV at time t.In the models by Perelson, et al. (1996) and Tan and Xiang (1999), the variance and co-variance structures were not obtained; only the expected variable values and their estimates were obtained.The numerical results reported   in Section 5 showed the efficacy of our model.We have not included t = 0 (after treatment) in our results; this is the time of pharmacokinetic delay, i.e. the time required for drug absorption, drug distribution and penetration into target cells (it varies from person to person) (Perelson, et al. (1996)).
We have used estimates extrapolated from clinical data of Perelson, et al. (1996) and Tan

Figure 1 :
Figure1: Simulated mean number of free HIV, infectious free HIV and non-infectious free HIV before and after combined therapeutic treatment.(Units: On the horizontal axis 0.5 is 1 hour, while on the vertical axis viral counts are measured in copies/m of blood.)

Table 1 :
Assumed values of parameters used to obtain numerical results.

Table 6 :
M XV (t), M XD (t) and M V D