A Mathematical Model for Treatment-Resistant Mutations of HIV

In this paper, we propose and analyze a mathematical model, in the form of a system of ordinary differential equations, governing mutated strains of human immunodeficiency virus (HIV) and their interactions with the immune system and treatments. Our model incorporates two types of resistant mutations: strains that are not responsive to protease inhibitors, and strains that are not responsive to reverse transcriptase inhibitors. It also includes strains that do not have either of these two types of resistance (wildtype virus) and strains that have both types. We perform our analysis by changing the system of ordinary differential equations (ODEs) to a simple single-variable ODE, then identifying equilibria and determining stability. We carry out numerical calculations that illustrate the behavior of the system. We also examine the effects of various treatment regimens on the development of treatment-resistant mutations of HIV in this model.

1. Introduction.Human immunodeficiency virus (HIV) is a pathogen that infected approximately five million people worldwide in 2003 [34].Progression from HIV infection to acquired immunodeficiency syndrome (AIDS) typically occurs over a decade or two.In 2003, almost three million people died from AIDS; since the first cases were identified in 1981, more than 20 million have died from AIDS.The immune system and, in particular, the T cells, play a central role in HIV population dynamics, including the progression to AIDS [37,8].There are various treatments in use, but no available treatment protocol results in clearance of the virus from patients.Current drug treatments can extend the healthy life span of infected patients by years.
One of the main problems is the emergence of drug-resistant strains of HIV in patients undergoing treatment [31].These mutant strains result in the resurgence of viral loads after long-term suppression of the loads from treatment.Viral-load resurgence correlates with T cell depletion and with the progression of patients to AIDS [37].In this paper, we apply mathematical modeling and analysis techniques to examine the properties of a model that includes drug-resistant strains.The understanding gained from this work will be applied in future work to optimize treatment regimes in this setting, with the goal of maintaining patient health for as long as possible.
There is an extensive body of work that develops models of this type for the interaction of T cells with HIV.Alan Perelson and Patrick Nelson [29] and Martin Nowak and Robert May [25] have written detailed surveys of the main ideas developed through such models.Several other diseases with similar pathogen-immune system interaction dynamics have been modeled in this way as well.These include hepatitis B (cf. [23]), hepatitis C (cf. [22]), tuberculosis (cf.[38]), and chronic myelogenous leukemia (CML) (cf.[21]).
Among the previous models that have considered mutant strains of HIV, most consider only a single drug-resistant strain, and very few allow for mutations between strains.The contribution of the model presented in this paper is the analysis of a model that includes four categories of differing drug resistance, with mutations allowed between all of the categories.The virus categories we consider are wildtype, resistant to treatment in the protease inhibitor class, resistant to treatment in the reverse transcriptase inhibitor class, and resistant to treatment in both classes.Our analysis is intended to add to a broader understanding of the dynamics of multiple drug-resistant strains, which could help guide future combination treatments in this complex situation.
In addition to the analytic results we present, we also obtain numerical results, first without any treatment, and then with various treatments.We consider constant dosing protocols of various strengths, as well as structured treatment interruptions, with various interruption lengths.Our intent here is to demonstrate some of the possible outcomes over fixed time periods.In agreement with models that do not consider the two drug treatment classes separately (for example, [15]), we find that the length of the cycle in the treatment interruption schedule does not qualitatively alter the resulting levels of HIV or T cells.We also find that treatment appears to exert selective pressure in favor of drug-resistant strains.
2. Background and assumptions for the model.We use differential equations to model the dynamics in the peripheral blood of a hypothetical patient.The equations give rates of change for various T cell populations, with parameter values obtained from available experimental data and from estimates.
We consider healthy T cells (which are not infected with HIV) and T cells that have been infected with several different strains of HIV.The HIV strains we consider are wild-type virus (which is sensitive to the two drugs considered here), virus that is resistant to one type of treatment, virus that is resistant to a second treatment, and virus that is resistant to both drugs.We assume that resistance is acquired through mutations, mainly from errors that occur during the transcription process.Later, the drugs we consider will be one from the class of protease inhibitors (PIs) and one from the class of reverse transcriptase inhibitors (RTIs).
We do not include latently-infected T cells (those that will begin production of virions after a significant delay) or nonproductively infected T cells (those that will die before producing very many virions).These populations have been considered in other models, but do not appear to significantly alter the dynamics of the system.(Cf.[27] and [28] for the consideration of latently-infected T cells; cf.[30] in the case of nonproductively infected T cells.)We also do not include T cells that are infected by more than one type of HIV represented here, with the assumption that this number is much smaller than the other populations.
Throughout this paper, we refer to the T cell populations infected with different HIV strains using binary notation: quantities related to wild-type, drug-sensitive HIV have the subscript " 00 "; those related to the first class of strains have the subscript " 01 "; those related to the second class of strains have the subscript " 10 "; and those related to the HIV population resistant to both drugs have the subscript " 11 ".Later, we will have the first class of strains (01) represent those resistant to a protease inhibitor and the second class of strains (10) represent those resistant to a reverse transcriptase inhibitor.
Our model is based in the circulatory, or peripheral, blood system, since we wish to tie it as closely as possible to available data from blood samples.We include a source term for new T cells entering the circulatory blood from other compartments (such as the bone marrow, lymph nodes, and thymus).However, we assume that preexisting T cells entering and exiting the blood due to diffusion give net changes that are approximately zero.We also assume that the concentrations of T cells infected with the various strains of virus are approximately proportional to the concentration of free virus in the circulatory blood.This is a good approximation for a system at or near equilibrium, and experimental support for this approximation in the case of HIV appears in [13] and [27].This assumption of an equilibrium state will allow us to consider infected T cell populations, without including the corresponding free virus populations.To make this assumption, we restrict our model to the middle stages of HIV (after the acute phase) and before the initial stages of AIDS.
We assume the T cells come into contact with virus in the blood in a random fashion.Because the encounters take place in the blood, we use the "law of mass action", which says that the total number of encounters between members of the two populations is proportional to the product of the sizes of the two populations.See [9] for more discussion of the law of mass action.
3. Details and explanation of the model.We consider five populations of T cells in the circulating blood system.The first is the population of T cells that is uninfected with any HIV.The other four are T cells that are infected with various populations of HIV, as described below.All cell populations are measured in concentrations of cells per µl and are functions of time, t, which is measured in days: T = uninfected T cells, I 00 = T cells with drug-sensitive wild-type virus, I 01 = T cells with virus with mutations of type 1 and not of type 2, I 10 = T cells with virus with mutations of type 2 and not of type 1, I 11 = T cells with virus with mutations of types 1 and 2.
The system of differential equations is given below, followed by explanations of the terms: Each equation represents the rate of change, with respect to time, of one of the populations.The lowercase coefficients, or parameters, (e.g., s and µ 01 ) are all taken to be constants, as is T max .Later in this paper, we vary the parameter values and examine the resulting changes in the system.
Figure 1 shows the cell population diagram for the system (1)-(5).We assume the changes in populations due to diffusion are approximately zero.That is, we assume that the numbers of the various populations that move into or out of other compartments (other than the peripheral blood) are approximately equal.Thus, there are no terms in the equations for diffusion of T cells into or out of the blood.
The first term on the right-hand side of equation ( 1) is a source term for new T cells entering the blood system.We approximate this as being a constant, s, during most of the intermediate stages of infection with HIV.The second term is a logistic term with logistic growth rate p and limiting value T max , as originally proposed in [14].We use all of the T cell populations considered here (uninfected and infected with HIV) in the logistic population count.The third term is the loss due to the natural attrition in the T population.The factor δ T is the death rate constant of the T population, and is equal to the reciprocal of the average life span of cells in the T population.In the absence of infection, these three terms maintain homeostasis of the T cells in the peripheral blood [32].
The last four terms in equation ( 1) represent losses of T cells from the peripheral blood due to interaction and subsequent infection with various types of HIV.Each of these terms is in mass-action form, as they represent contributions due to encounters between T cells and HIV occuring in the peripheral blood, which we assume is well mixed.Recall that we assume the concentration of free virus of each type is proportional to the concentration of T cells infected with each type of HIV.Hence, β ij I ij includes a constant of proportionality that gives the concentration of free virus of type ij , as well as the law of mass action constant, which incorporates the expected frequency of encounters between T cells and HIV, and the expected fraction of these encounters that subsequently lead to infection by HIV.
We now examine the four equations ( 2)-( 5), which describe the rates of change of T cells infected with the four types of strains of HIV we consider here.In each of these equations, the first term is a contribution that equals a loss term for T cells in equation (1).For example, in equation ( 2), the first term is equal to the rate of loss of T cells that become infected by wild-type HIV.The second term in each of the four equations is a loss term δ ij I ij due to the natural life span of T cells, with δ ij the death-rate constant for T cells infected with HIV of type ij.We assume that δ ij is at least as large as the death rate of healthy T cells, δ T .
The constants µ ij , for ij = 01 or ij = 10, give the average mutation rates of the mutations that lead to HIV with mutations of type 1 or type 2, respectively.These mutations are assumed to be independent, and hence the mutation rate for HIV with mutations of both type 1 and type 2 is given by µ 01 µ 10 .The factors (1 − µ 01 ) and (1 − µ 10 ) give the fractions of virus that are expected not to mutate to strains of type 1 or type 2, respectively, at any given time.The constants µ ij , for ij = 01 or ij = 10, represent the backward mutation rates.For example, µ 01 represents the mutation rate of strains of type 1 back to wild-type virus.
The constants κ ij are the net-gain scaling factors, or effective burst-rate constants.These take into account the rate at which bursts release large numbers of free virus, as well as the rate of any releases of smaller scales, which may be considered to take place continuously.As an example, the term κ 00 (1 − µ 01 )(1 − µ 10 )I 00 in equation (3) gives the contribution to the I 00 population due to the "bursting" of the infected T cells I 00 that do not mutate to one of the other types of strains considered here.Similarly, the term κ 01 µ 01 I 01 in equation (3) gives the contribution to the I 01 population due to the "bursting" of I 01 T cells that mutate from strains of type 1 to wild-type virus.
It is the assumption that the distribution of virions is at or near equilibrium that allows us to avoid equations for the free virus.This assumption is used in the last four terms of equation ( 1), when we use I ij in place of a constant times the corresponding virus population.We also use this assumption in equations ( 2)-( 5), when we write the contribution terms from production of virions as if the "bursting" of the infected T cells contributes directly to the infected T cell populations.In this instance, the constants κ ij include factors that give the appropriate rate at which virions could move from being produced by I ij to producing new infected cells I ij .
4. Discussion of parameter estimates and initial population values.The values of the constants used in the system of differential equations above appear in Table 1.This table gives a brief description of the constants, the values used as initial estimates for the constants, the ranges that are used later in the paper, and the units and references for the constants.We use the following initial values for the populations, where t = 0 is the starting time for the model: The value we use for s, the source term for the uninfected T population, is calculated using data from [19].In the data collected from subjects with HIV, the average contributions from a source were found to be 0.025 CD4 + cells per µl per day, and 0.023 CD8 + cells per µl per day.This gives a net contribution of 0.048 cells per µl per day, which we use as the source term for the population,T .
We base our logistic growth rate p on the same calculation used in [35].Namely, we choose p so that it agrees with the increase by 95 cells per µl in six weeks in [27], assuming that treatment is highly effective.Here, we are taking into account that s = 0.048.
For T max , the expected maximum number of uninfected T cells per µl, we note that Hiroshi Mohri et al. [19] found average numbers of CD4 + plus CD8 + cells to be close to 1,700 cells per µl in healthy patients.The highest count they found was 2,017 cells per µl, so we use 2,000 as our best estimate, with a range that allows values up to 2,500 cells per µl.
The value we use for the death rate constant, δ T , of uninfected T cells is the same as that used in [35], which assumes that T cells live approximately two years on average.For the death-rate constants δ ij for the infected T cell populations, we use values close to the mean reported in [30] for patients with HIV.We assume that the T cells infected with wild-type virus have the lowest death rate of the four.
The values for β ij are based on several estimates.We start with the estimate in [35] that the ratio of virions to infected T cells in the peripheral blood is approximately 112.We make very rough estimates of mass-action encounter rates, as in other T cell-pathogen models, such as [10].Finally, we incorporate the estimated infectivity rate of 3.43 × 10 −8 in [35].
For the mutation rates, Louis Mansky and Howard Temin [16] showed that the forward mutation rate for strains of HIV resistant to RTIs is approximately 3×10 −5 per base pair per cycle.As in other models of drug-resistant mutations (e.g., [3] and [35]), we assume the forward and backward mutation rates are the same.With this assumption, the pressure for forward mutations occurs because of the increased fitness of the drug-resistant strains under treatment conditions.We also assume that µ 01 < µ 10 .There is some evidence that mutations resistant to RTIs are more prevalent than mutations resistant to PIs. (See, for example, [1] and [2].)Although this could alternatively be achieved by assuming that RTI-resistant virions have higher fitness than PI-resistant virions (i.e., that η RTI < η PI ), our simulations favor the former assumption, namely that µ 01 < µ 10 .
In [19], the average baseline concentration of CD4 + plus CD8 + T cells in HIV patients was approximately 1,200 cells per µl, which we use for T (0).We arbitrarily choose the initial values of the infected T cells, with the assumption that there are more infected with wild-type than with resistant strains.

Analytic results.
To study the behavior of the system of differential equations, we first show how to use the system to get a (nonlinear) second-order ordinary differential equation (ODE) in T , with two initial conditions.This ODE is used later to obtain the fixed points of the complete system.Once we find any fixed points, we can examine the stability at each of them.For ease in later computations, we write this equation in vector form as where we use bold type to indicate the vectors α = (α 00 , α 01 , α 10 , α 11 ) and I = (I 00 , I 01 , I 10 , I 11 ).Next, we rewrite the four equations ( 2)-( 5) in matrix form as shown: Here, T = T (t) is the healthy T cell population as before, I = I(t) is the vector defined above, and A and B are the following constant matrices: where Let M(t) = A + T (t) B. Since A is in M 4×4 (R), there exists an invertible matrix P ∈ M 4×4 (R) such that P −1 A P = J , where J ∈ M 4×4 (R) is the Jordan canonical form of A. So A = P J P −1 , and we can rewrite the matrix M(t) as follows: M(t) = P J P −1 + T (t) B. The only assumption we make to proceed here is that β 00 ≈ β 01 ≈ β 10 ≈ β 11 ≈ β.This gives the following: where I is the 4 × 4 identity matrix.This allows us to rewrite equation ( 8) as follows: Since A is a constant matrix, P does not depend on t, so Letting X(t) = P −1 I(t), we get In general, J is diagonal, since the the set of diagonalizable matrices is dense in the set of real n × n matrices.We will assume that this is the case so that the diagonal entries of J are the eigenvalues of A. We will denote these eigenvalues by λ 1 , λ 2 , λ 3 , and λ 4 .Note that the method below still works if J is not diagonal, since we can still find the eigenvalues of A, but the calculations are more cumbersome.Now, we will focus on solving (15).Once we solve (15), we will know X(t) in terms of T (t), and hence we will know I(t) in terms of T (t).We can then substitute the expression for I(t) in terms of T (t) back into equation (8), to obtain a system of ordinary differential equations in one variable (T ), which we can then solve.
Again, using the assumption that Using commutative properties of diagonal matrices and the fact that the exponential of a diagonal matrix is a diagonal matrix, it can be checked that is a solution to equation ( 16), where τ is a dummy variable and C is a constant vector.To find C, we plug in t = 0 to (17) to get X(0) = C and recall that X(t) = P −1 I(t).So C = P −1 I(0), and From this expression for X(t) and the fact that I(t) = P X(t) we get Let F (t) = t 0 T (τ )dτ .Then I(t) = P e t J e (F (t) B) P −1 I(0) Let Then Plugging this expression for I = I(t) into equation ( 7) and letting α ij = α 2i+j+1 gives We use a Taylor series to expand the exponential factors E k as and plug these back into (22): For T = 0, we can divide both sides by T : Recall that we defined F (t) = t 0 T (τ )dτ .By the fundamental theorem of calculus, we have F (t) = T (t).We will use this as we take the derivative of both sides of (25): Multiplying both sides by −T 2 and rearranging, we get where We let T (0) = T 0 , which gives one initial condition.To get a second initial condition, we use equation ( 22) and the fact that F (0) = 0: Equation ( 27) and the two initial conditions T (0) = T 0 and (29) give a well-posed ODE initial value problem.Once we obtain T values (by solving analytically or numerically), we then can use (20) or (21) to find I(t).

Equilibrium analysis.
We now set each of the five derivatives in (1)-( 5) equal to zero and solve for T and I ij , for i, j = 0, 1.This gives the fixed points, or equilibrium solutions; that is, it gives values of T and I ij for which the system will no longer change (since all of the derivatives, or rates of change, will be zero).
We first consider the equilibrium solution(s) for which I ij = 0, for i, j = 0, 1.In this case, equation (1) implies that This is the unique healthy equilibrium, which we denote by W 0 = (T h , 0, 0, 0, 0).For other equilibria, if I = 0, we use (8) to get Since we assume that the β ij 's are nonzero, B is invertible, and we have the following eigenvalue problem: where C = B −1 A, I = 0. We will assume that T = 0, since T = 0 is outside of the biologically relevant region we wish to consider.The matrix C is a 4 × 4 matrix with real-valued entries, and so it has eigenvalues that we will denote by σ 1 , σ 2 , σ 3 , and σ 4 .From (32), we know that the values taken by −T must be σ k , for k = 1, 2, 3, 4. (Recall that we have set T (t) = 0, so T is now a constant.)Let be an eigenvector associated with eigenvalue σ k .Since any nonzero multiple of V k will also be an eigenvector with eigenvalue σ k , we have     for 0 = r k ∈ R.
Since T (t) = 0, equation (7) implies that which implies that With this value of r k , we now have enough information to find the specific eigenvector r k V k for the eigenvalue σ k for the eigenvalue problem above.From this, we get a total of four possible fixed points, given by for k = 1, 2, 3, 4.

Stability analysis.
To determine the behavior of the cell populations near the equilibrium solutions, we need to compute the linearization of the system, which is obtained from the Jacobian matrix of the system.See, for example, [6, p. 446] for more information about such calculations.If we let A = a kl , then for the system of equations ( 1)-( 5) the Jacobian Df is the following matrix: We first consider stability of the fixed point W 0 = (T h , 0, 0, 0, 0).In this case, the Jacobian Df (W 0 ) is the following matrix: For this fixed point, using the parameter values in Table 1 and the expression for T h given in equation (30), one of the eigenvalues of this matrix is positive.This is still true even when various combinations of the parameter values are varied over an order of magnitude or so.This indicates that the healthy equilibrium is unstable for this model, as is the case in patients: a small viral load is expected to eventually lead to full-blown AIDS.We numerically compute the other possible fixed points, W 1 , W 2 , W 3 , and W 4 , using the parameter values in Table 1.To do this, we first compute the eigenvalues and eigenvectors σ 1 , σ 2 , σ 3 , and σ 4 and V 1 , V 2 , V 3 , and V 4 , respectively, of the matrix C = B −1 A, where A and B are defined as in ( 9) and (10).We plug these into equation (36) and then compute the equilibrium solutions W k as in (37).We find there is exactly one equilibrium value, which we will denote by W + , that has a nonnegative value for each of the cell populations.In this case, the I 00 population dominates all others.The eigenvalues of the Jacobian for the solution W + all have negative real parts for a wide range of parameter values tested, so we conclude that this equilibrium is asymptotically stable.This matches the behavior shown by the numerical solutions in Figure 2.
6. Model with treatment.Here, we consider the same model as above but with treatment.We let η RTI denote the efficacy of treatment with a reverse transcriptase inhibitor and let η PI denote the efficacy of treatment with a protease inhibitor.Since 7.1.Without treatment.First, we show graphs of the solution curves to equations (1)-( 5), with no treatment.All parameter values are those that appear in Table 1, unless otherwise noted.All T cell counts are per µl along the vertical axes, and the horizontal axes show time in days.In agreement with our analytic results in section 5.3, we see that T cells infected with the wild-type virus dominate all others and approach a stable equilibrium value.7.3.Discussion of results of numerical analysis.The results of the numerical analysis in the previous section are representative of the variety of behaviors under different treatment regimes, and with various parameter values.The behavior of the graphs could vary greatly if parameter values were varied, even if by only a few percentage points, regardless of the treatment scheme.Surprisingly, we found that structured treatment interruptions (STIs) did not differ qualitatively from constant-dose treatments.The numerical solutions obtained using STIs are very similar to those with constant dosing schemes, as shown in Figures 3-5.We found this to be true for STIs of varying lengths, from cycles of seven days on and seven days off, to cycles of 120 days on and 120 days off.
For comparison, we cite two studies that confirm such results in the clinic.The retrospective clinical study by José Moltó et al. [20] found no statistically significant difference after forty-eight weeks in outcome between a group of patients that had undergone STIs (six cycles of two weeks off and four weeks on) and a group that had been on constant-dose therapy without STIs.A study by Annette Oxenius and Bernard Hirschel [26] compared cohorts on STIs with those not on STIs.Their findings also indicate that STIs does not lead to a statistically significant difference in outcome for patients who started STIs during the chronic phase of HIV infection.However, they did find statistically significant differences in outcome between cohorts in the acute phase of infection, with STI patients faring better.This indicates that our model would need to be changed to model the acute phase of HIV infection accurately.8. Conclusions.We have proposed a new model for HIV with two different types of resistance to treatments.We have performed the preliminary mathematical analysis on the model, and have tested the model qualitatively using numerical simulations.Data collection to determine more of the unknown parameter ranges would enhance the results.This model could also be used and refined to include other aspects of HIV dynamics interaction.One example is the incorporation of the effects of immune system-boosting treatments in patients with HIV.
Having this model now allows us to set up the optimal drug-dosing problem for PI and RTI treatments, which is the subject of our next paper.Given estimates of the mutation rates of HIV, the most effective dosing levels over time can be predicted for treatments.In the case of unlimited resources, optimal dosing could substantially increase the period in which a patient remains healthy, which is especially important for treatments to which diseases eventually develop resistance.Since drugs for HIV can cost thousands of dollars per month, a constrained optimal dosing solution, for limited resource scenarios, would be valuable as well.Further work predicting optimal combinations of PIs and RTIs would also be of great benefit in the treatment of HIV.

Figure 1 .
Figure 1.Cell population diagram.In this figure, solid curves indicate changes solely dependent on the originating population; dashed curves indicate changes that depend not only on the originating population but also on the destination population.The movement of uninfected T cells into the infected populations takes place due to infection of the T cells by virus with the respective mutations.Movement between the populations of T cells that are infected with various strains of HIV occurs due to mutation in the viral genome.

Figure 2 .
Figure 2. Qualitative T cell population behavior in the absence of treatment.Uninfected T cells are shown in black, T cells infected with wild-type virus are shown in purple, T cells infected with strains resistant to type 1 treatment are blue, those resistant to type 2 treatment are green, and those resistant to both types of treatments are red.Cells per µl are on the vertical axis, and days are on the horizontal axis for this and all following graphs.

7. 2 .Figure 3 .
Figure 3. T cell populations with constant-dose treatment.Here, we set η PI = 0.1 and η RTI = 0.25.Uninfected T cells are black, those infected with wild-type virus are purple, those resistant to type 1 are blue, those resistant to type 2 are green, and those resistant to both are red.

Figure 4 .Figure 5 .
Figure 4. T cell populations with stronger constant-dose treatment.Here, κ 00 = 0.57, κ01 = 0.44, and κ10 = 0.44, and we set η PI = 0.2 and η RTI = 0.25.Uninfected T cells are black, those infected with wild-type virus are purple, those resistant to type 1 are blue, those resistant to type 2 are green, and those resistant to both are red.