Optimal Antiviral Switching to Minimize Resistance Risk in HIV Therapy

The development of resistant strains of HIV is the most significant barrier to effective long-term treatment of HIV infection. The most common causes of resistance development are patient noncompliance and pre-existence of resistant strains. In this paper, methods of antiviral regimen switching are developed that minimize the risk of pre-existing resistant virus emerging during therapy switches necessitated by virological failure. Two distinct cases are considered; a single previous virological failure and multiple virological failures. These methods use optimal control approaches on experimentally verified mathematical models of HIV strain competition and statistical models of resistance risk. It is shown that, theoretically, order-of-magnitude reduction in risk can be achieved, and multiple previous virological failures enable greater success of these methods in reducing the risk of subsequent treatment failures.


Introduction
The development of multi-drug regimens for HIV therapy has resulted in HIV infection becoming a chronic, manageable disease in first world countries [1]. The necessity of a three-drug regimen, where each drug in the regimen targets separate viral epitopes, is due to the extremely high replication and mutation rates characteristic of HIV infection [2]. These make the evolution of viral strains resistant to a single drug inevitable. Three drugs, however, present a mutational barrier high enough to make such an evolutionary occurrence unlikely [3,4]. These combinations contain three drugs from at least two separate classes of antivirals, including the nucleoside/nucleotide analog reverse-transcriptase inhibitors (nRTI), non-nucleoside reverse-transcriptase inhibitors (NNRTI), protease inhibitors (PI), and integrate inhibitors (II). While these three-drug regimens, known as highly active antiretroviral therapy, or HAART, are highly effective at suppressing the virus in the long term, some patients nevertheless experience viral load rebound, driven by the emergence of a viral mutant resistant to all three components of their HAART regimen.

Mutation
Mutation events in HIV replication appear to be dominated by point-substitution events, which occur with very high frequency. This, coupled with the high turnover rate of HIV in uncontrolled infection, create a situation in which multi-drug resistant virus develops frequently. When a resistant mutant emerges, it becomes necessary to switch to a new three-drug regimen, whose components exhibit no cross-resistance with the failed three-drug regimen [5]. There are a limited number of independent drug combinations. A patient who has developed viral strains resistant to all such combinations is called Multi-Drug Resistant or MDR, and such patients are left with few viable treatment options. It is critical, therefore, to preserve the remaining pool of independent HAART regimens, especially for patients who have experienced virological failure on more than one previous regimen.
Attempts have been made to re-sensitize the virus to previously failed regimens through the use of treatment interruptions; the theory is that the wild-type virus, which enjoys a competitive advantage in the absence of therapy, would re-establish dominance and potentially drive the resistant virus extinct through competition [6]. Although these studies showed a brief return of susceptibility, the resistant strain quickly returned upon re-introduction of the drug regimen, and overall patient outcome was worse than a non-interrupted control group. More recent approaches have focused on changing the genetic makeup of the viral pool in MDR patients in preparation for 4-9 drug rescue regimens known as Mega-HAART or giga-HAART [7,8], [9], [10], [11], [12]. These approaches showed mixed results, mostly with poor clinical outcomes. All of these previous approaches attempted to use treatment interruptions to manipulate the susceptibility of the virus to regimens consisting of drugs to which resistant virus had already emerged. None of these addressed the possibility of using interruptions to preserve the usefulness of a naive antiviral regimen. Also, the antiviral regimen introduced following the interruption was always novel, implying an attempt to manipulate susceptibility by genetic profile alone, as opposed to manipulating viral load and genetic profile.
Attempts have also been made to use previously failed drugs in novel combinations in order to preserve some usefulness from previously failed treatments in MDR patients. The problem with this, however, is that the existing mutations represent a lowering of the mutational barrier. The only way to overcome this in the longterm seems to be an increase in the number of drug components used, which may succeed at the goal of reducing viral load at the cost of raising the side effects of the HAART drugs to an unacceptable level. A notable exception to this was the recent Tetriz study [13]. In this study, a drug combination using components from previously failed regimens, including two drugs for which the resistance mutations were known to be antagonistic. Despite the use of only four previously failed components, this regimen succeeding in inducing complete viral suppression in a significant portion of the study group, strongly suggesting the usefulness of permuted regimens. The importance of preserving suppressive regimens has driven a number of clinical studies, including the SWATCH study [14], which showed reduced incidence of virological failure in patients undergoing a preemptive switching schedule based on mathematical models of risk similar to those described in the Analysis section.

Competition and Selection
The development of drug resistance in HIV infection is driven by two phenomenon: mutation and selection. Mutation in HIV replication occurs at a well-characterized, relatively constant rate of approximately 3 Ã 10 {5 substitutions per base-pair per replication cycle [15]. Other mutation types, such as deletions, insertions, and rotations, happen with considerably lower frequency, and do not usually contribute to the development of resistance. Despite this relatively high rate of mutation, the population of virus in a treatment-naive patient contains only virus with very few genetic changes from the nominal, or ''wild-type'' HIV sequence. This is due to the influence of selection; the wildtype dominates in the absence of treatment because it is usually the fittest virus in that environment. Many mutations carry a fitness cost when compared to the wild-type sequence; viruses carrying these mutations do not replicate as efficiently as the wild-type virus. The various virus subtypes compete for target cells, so selective pressures tend to drive extinct virus variants that carry mutations.
When antiviral medication is used, the wild-type is no longer the fittest virus; their interference with the virus' ability to replicate is such that the virus population will shrink exponentially. Various mutations exist that, if present, reduce the ability of the antiviral drugs to interfere with HIV replication; if they interfere to the extent that the mutant virus population is able to grow overall, the mutation can provide clinically significant resistance.

Genetic Distance and Fitness Cost
The likelihood of a particular resistance mutation emerging is influenced by two major factors. The first is the relative fitness of the mutation under the current treatment. This may be calculated by considering the relative effectiveness of the mutation at negating the effect of the drug and the relative fitness cost of the mutation. Fitness cost, in this sense, means the decrease in the viruses' ability to effectively replicate in the absence of treatment as a result of the mutation.
The second factor influencing likelihood of emergence is the genetic distance of the resistance mutation from the existing virus pool. This is the number of point mutations necessary to generate the resistance mutation. If the HIV genome is considered to reside in a sequence space, the genetic distance is equivalent to the Hamming distance. Because mutation is a random process, mutations with a high genetic distance are very unlikely to emerge.
Example Strains from the Stanford database. Extensive data on resistance mutations to the antivirals listed above has been collected online at the Stanford HIV database [16]. An example from the database can illustrate the genetic distance calculations referenced above. Consider a patient who developed viral resistance to an initial therapy consisting of the NRTIs abacavir and lamivudine and the NNRTI nevirapine (this is a standard firstline therapy). According to the database, one set of mutations yielding significant resistance to these three drugs is (74V,103N,184V), that is, a substitution of valine for leucine in position 74 of the viral reverse transcriptase protein, a substitution of asparagine for lysine in position 103, and a substitution of valine for methionine in position 184. Together, these mutations require at least three point substitutions from the wild-type genome (the number may be higher, as multiple sequences may code for the same amino acid), giving us a genetic distance of 3.
Since the (103N) mutation gives broad class resistance against all NNRTIs, any follow-up therapy will not use NNRTIs. Neither (74V) nor (184V) exhibit strong cross-resistance patterns with any other NRTIs, so a possible follow-up therapy would be the two NRTIs tenofovir and zidovudine together with the PI nelfinavir. Clinically significant resistance to these three drugs could be conferred by the set of mutations (41L,210W,215Y) on the viral RT protein and the mutation (30N) on the viral protease protein. This set of mutations has a genetic distance of 4 from the wildtype, but a genetic distance of 6 from the mutant virus generated in the first round of treatment. This is because the inclusion of either the (74V) or the (184V) mutation increases the susceptibility of the virus to both tenofovir and zidovudine. Consequently, any resistant virus arising to the second treatment will probably arise from the wild-type viral pool, and will probably not carry the (74V), (103N), or (184V) mutations.
Therefore, after developing resistance sequentially to these two treatment regimens, the patient will be carrying three viral strains; the wild-type, one mutant carrying the (74V,103N,184V) RT mutations and another mutant carrying the (41L,210W,215Y) RT mutations and the (30N) Protease mutation,. A resistance mutation to a third antiviral schedule consisting of a permutation of the first two regimens, say abacavir, zidovudine, and nevirapine, would have to arise from one of these parent viruses. The common resistance mutation to this combination with the closest genetic distance to the wild-type carries the mutations (70R, 103N, 184V, 215F) would have a genetic distance of 5. This is also the variant with the shortest genetic distance to the first mutant, with a genetic distance of 4 (the inclusion of the 74V mutation eliminates the resistance to zidovudine, so it must be undone to confer resistance). The variant with the shortest genetic distance to the second mutant would carry the RT mutations (41L,103N,184V,210W,215Y), and would have a genetic distance of either 2 or 3, depending on whether the proteaseresistance mutation renders it inviable.
Genetic Distance Uncertainty. The computation of genetic distance between HIV strains is relatively straightforward, but there are a few ways in which the genetic distance can be over-or under-estimated. The first is the non-uniqueness of the genetic code; multiple genomic sequences can code for the same amino acid. The calculations carried out in the previous section assumed the parent genome contained the sequence with the shortest genetic distance to the mutant; this may lead to an underestimation of genetic distance.
When using genetic distance to estimate mutational barriers, as described in the next section, the existence of viable transitional forms can result in an overestimation of genetic distance. That is, while the genetic distance to the identified resistant strain may be high, an unidentified partially resistant strain with a lower genetic distance may provide a ''stepping-stone'' for the development of the fully resistant strain.

Computation of Risk as a function of Viral Load
Pre-existence. The research of Bonhoeffer and Ribeiro [17], [3] show that emergence of resistant virus strains is most likely caused by preexistent resistant mutants under very general conditions. Bonhoeffer [17] also stated that the preexistence does not mean that there is a stable coexistence of sensitive and resistant virus. The preexistence of resistant mutant is made from mutations between sensitive and resistant mutant. Therefore, in order to quantify the risk of drug-resistance emergence, resistance mutations must be modeled as a stochastic process.
Poisson Modeling. In this section, equations are presented that determine the drug-resistance emergence risk, which is the mutation probability from the current virus pool to a resistant mutant for the new regimen. To accomplish this, a Poisson distribution is used to model the mutation process. The probability of pre-existing resistant genotypes depends upon two key variables: the total number of active virions for each type of virus present in an infected host, v 1 ,v 2 ,:::v n , and the point genetic distance from the current virus strains to the emergent resistant mutant v e : g 1 ,g 2 ,:::g n . Based on the research of [18], [19], [20], [21], the total viral burden can be estimated. The total viral burden of actively circulating virus can be roughly calculated as the viral titer multiplied by the volume distribution of total body extracellular fluid which is 25-30% of body mass [21]. For example, a 100-kg man (roughly 100-liter volume) with a viral RNA load 10000 copies/ml would have (100000| 10000|0:3) circulating virions approximately.
Therefore, assuming a point mutation rate of m, the probability of drug-resistance emergence risk is calculated as follows: where v i (S t ) is the viral load of virus subtype v i present in the patient at the time of introduction of the naive regimen, S t . Consider the simplest case. Assuming that the current virus pool consists of only one kind of virus strain and the genetic distance from the current virus strain to the resistant mutant is either 1, 2 or 3. Fig. 1 shows the relationship between viral load and the risk of resistant mutant emergence.
As Fig. 1 shows, the resistant mutant may pre-exist when the point genetic distance is 1 or 2, but the pre-existence of resistant mutants is very unlikely if the point genetic distance is 3 or bigger. Consider the case where the current virus are all wild-type and the point genetic distance is 2 between the wild-type virus and a mutant resistant to the naive regimen. If the patient is switched to the naive regimen when the viral load is 30000 copies/ml, the probability that the resistant mutant will preexist is 64.11%. However, if the switch is made when the viral load is 2000 copies/ mL, the risk is only 6.86%. Therefore, the task is to create a switch point for the new regimen with the lowest risk.

Model
To model HIV dynamics, a set of ordinary differential equations is used that includes terms describing the mutations among different virus types. This model depicts the interactions between a wild-type virus population sensitive to all antiviral drug regimens and a resistant mutant virus population only sensitive to treatment with some, if any, antiviral drug combinations. The model is in shown in Equation 2 This model includes states x, representing CD4+ T cells that are susceptible to infection; y k , the CD4+ T cells infected by the virus type v k , and v k , the k th type of virus in the patient's virus pool. k~w for wild-type virus and k~r for resistant mutant viruses. The definition of each parameter, and its units, may be found in Table 1.
CD4+ T cells are generated from their source at rate l and disappear at rate d. The target cells are infected by the viral strain v k at rate of b k and the therapy suppresses the infection by strain v j with efficacy j k,i u i , where j k,i is the maximum effectiveness of antiviral regimen u i on virus strain v k .
The infected CD4+ T cells y k are created first by the infection from target cells x with virus v k that are generated at a rate b k . Secondly, there is a contribution due to the activation from longlived reservoirs at rate l k . The infected CD4+ T cells die with a rate of a k .
The virus v k are created from infected CD4+ T cells y k . Virus, v k , die with a rate of v k . These equations are arbitrarily scalable to any number of viral species.
Model identification parameters. To apply the system of equations describing the evolution of viral loads and CD4+T cells to a specific individual, the parameters of the nonlinear differential equations need to be estimated. Using patient data from the AutoVac study [22], a Bayesian estimation technique (specifically, the MCMC, Monte Carlo Markov Chain, method) is used to identify the parameters for this model. In the AutoVac study, HIV patients had viral load measurements taken at a 1-day interval during a series of short treatment interruptions. The data available for estimation is limited to a relatively few values of viral loads after an interruption of medication and reintroduction. In this work, the approach of Huang [23] is used in applying a Bayesian nonlinear mixed-effects model. For the simplified model, there are six parameters to estimate: l, the generation rate of the target cells, d, the death rate of the target cells, b w , the infection rate for the wild type virus, a w , the death rate of the cells infected with the wild type virus, c w , the number of viral particles emitted per infected cell, and v w , the death rate of viral particles. The generation rates for the wild and resistant virus from longterm reservoirs are assumed to be a small, known constant, and j is assumed to be 1. Since the data is so limited, the parameters for the resistant virus are not identified from the data. Instead, they are assumed for the purposes of this paper to be proportional to the parameters for wild-type according to the ratio of viral fitness, with a nominal ratio of 0.5. In practice, this ratio could be estimated from in vitro fitness data available for most common mutation patterns.
Parameter estimates were generated for each of 12 patients. For each patient, treatment was interrupted and after a period of time, the treatment was restarted. This cycle of interrupting and reinstating the treatment is repeated 3 or 4 times. The MCMC procedure produced 200,000 possible combinations of parameters that are consistent with the patients' data. From this result, the marginal probability densities for of the six parameters can be established. Among the 12 patients, 3 of them have no detectable virus after an interruption and another 3 appear to be subject to multiple large exogenous disturbances (likely viral blips) which are not accounted for in our model, and result in poor fits. The results of the other 6 patients by using this identification method are shown in Fig. 2.

Multiple previously failed therapies
In the case of patients who have failed one or more drug regimens previously, the need to preserve the remaining regimens becomes all the more important. Interestingly, the previously failed regimens provide additional control inputs which can be used to reduce the risk of failure for the new regimen at a lower systemic cost than is possible when only one failing regimen is available to use.
Regimen Cycling. The multiple failed regimens allow two options for achieving a transient reduction in viral load. Cycling through the previously failed regimens before returning to the currently failing regimen is one such option. Consider a patient with viral dynamics described by Equation 2 who has developed virus strains v r1 and v r2 resistant to two previous treatment regimens u 1 and u 2 . If the viral strains resistant to those regimens are susceptible to the current regimen, they will have decayed to very low levels, and will take some time to re-emerge. Assuming no cross resistance, the currently dominant viral strain is likely susceptible to the original drug regimen u 1 to which the strain v r1 developed resistance. Strictly speaking, so long as R 0 (v r1 ,u 1 )w R 0 (v r2 ,u 1 ) and R 0 (v r1 ,u 2 )vR 0 (v r2 ,u 2 ), where R 0 is the basic reproductive ratio defined as then a transient viral minimum significantly lower than the steadystate level of viral load may be achieved by this method [24].
Our approach can be formulated as an optimal control problem in two steps. In the first step, the allowable patterns of treatment cycling where either u 1 (t)~1 or u 2 (t)~1 at any time t are searched to find a treatment pattern that minimizes the cost function where P(v e (S t )=0jv i (S t )) is the cost function defined by Equation 1 and S t is the time to introduce a naive regimen. In Fig. 2, S t1 and S t2 represent the time to introduce a naive regimen in current treatment strategy and our proposed treatment strategy respectively. If the genetic distances between the closest strain resistant to regimen u 3 and v w ,v r1 , and v r2 , respectively are all equal, this optimization returns the treatment cycling schedule with the largest decrease in total viral load prior to introducing the naive regimen, as seen in Fig. 3. If a naive regimen was not introduced at S t2 , the viral load would rebound shown as red-dash line in Fig. 3. The treatment schedule may be fixed, or may change as the optimization horizon T M increases, based on the actual values of the parameters in Equation 4 [25]. Since treatment switching on intervals faster than 1 day is not practical, the controls space is discretized at 1 day intervals. This yields an optimization problem with a finite search space, which can be solved in real time using exhaustive search techniques. The second step involves robustly estimating the time at which the minimum in the risk is achieved, and switching to the naive regimen at this point. This study is a subject for future research. The achievable minimum according to this method will be limited by the initial load levels of the various resistant viruses v ri . These in turn are determined by the length of time they have been under suppressive therapy, and their relative prevalence in the viral reservoirs. Using a cycling approach to achieve a risk minimum requires tolerating a relatively high viral load for a short period of time, which may contribute to disease progression or increased resistance.
Permuted Regimen Introduction. A better option for inducing a transient minimum in the case of multiple previously failed regimens is to introduce a permuted antiviral regimen. It is very likely that that no strain exists which is resistant to permutations of the previously failed therapies. While an antiviral regimen consisting of the permuted components of previous regimens will not provide sufficient mutational barrier to be a viable long-term option, they will allow a dramatic transient reduction in the total viral load, and the corresponding risk of preexistent resistance.
In the mathematical formulation, the same cost function as Equation 4 is applied. However, the simplified model is not applicable for this case, because u 1 and u 2 are efficacies of drug cocktails. In this case, u i is the efficacy of an individual drug. Therefore, Equation 2 is modified as follows: In this model, there are multiple virus strains, v i , where i~1,2,3,:::, corresponding to wild-type and all existing resistant-type viruses. The variable l k represents the contribution from long-live reservoir, and u i is the efficacy of the i th individual drug. The optimization is performed to achieve the minimum cost function defined by Equation 4 by exhaustively searching the possible schedules for each individual drug.
When utilizing permuted regimens, the optimal switching strategy involves switching from the failing regimen directly to a permutation of previously failed regimens prior to introduction of the naive regimen. Every previously dominant strain will be susceptible to this permuted regimen, so this will result initially in exponential decline in the total viral load. However, the reduced genetic distance inherent in using components of previously failed regimens means that resistance to the permuted regimen is likely to emerge. The expected achievable risk reduction can be estimated by assuming that the strain resistant to the permuted regimen preexists, with initial conditions related to measures of genetic distance.
Assuming that the two previously dominant resistant strains v r1 and v r2 are resistant to drug combinations a+b+c and A+B+C, respectively, then virus resistant to a permuted drug combination such as A+b+c will pre-exist with expected initial viral load where m is the pointwise mutation rate for HIV and v A {v r1 j j H is the number of point mutations in virus variant v r1 necessary to generate a virus v A with resistance to drug A and v bc {v r2 j j H is the number of point mutations in virus variant v r2 necessary to generate a virus v bc with resistance to drugs b and c, respectively. Note that these are also the Hamming distances applied to the genetic sequences of the respective viruses. Fig. 4 shows the case where only one point mutation separates the dominant resistant strain from a strain resistant to the permuted regimen. Standard treatment introduces a naive regimen at switch point S 1 . By introducing a permuted regimen at S 1 , it is possible to achieve a greater than 2 order-of-magnitude reduction in viral load before introducing the naive regimen at S 2 . The permuted regimen provides insufficient mutational barrier to prevent resistance, so if a switch is not made, the viral load will rebound. The reduction in resistance emergence risk achieved by this intervention depends on the genetic distance of the dominant strain at the switch time to the closest strain with resistance to the naive regimen. Table 2 illustrates the achievable reduction in risk. It is clear that while a genetic distance of 1 rules out any successful intervention, a genetic distance of 2 or three allow a greater than 2 order-of-magnitude reduction in the risk. Especially in the case of a genetic distance of 2, this is a dramatic change in expected outcome.
Although current guidelines [26] suggest changing drugs when the virus load exceeds some threshold values (e.g., 1000-5000 copies/mL), finding the exact time when the viral load reaches the threshold is unlikely, because a patient during therapy is only tested for viral load every 3 or 4 months. The equilibrium value of the viral load can therefore be used as the comparison benchmark. It is worth noting that the genetic distance of 1 between the previous viral strains and a strain resistant to the permuted regimen, as used in this example, represents a worst-case scenario (if the genetic distance is 0, this method cannot be used). The example from the Stanford database in the introduction illustrates a real-world case where this distance could be 2 or higher, yielding even greater reductions in viral load and the corresponding risk of pre-existent resistance.
Frequent sampling for minimum finding. Both of these methods use optimization to find schedules that create a transient minimum in the total viral load. Successful implementation of these methods require accurately finding the time when this minimum occurs and switching to a naive regimen at that point. The exact time point of the achieved minimum may vary considerably from its calculated point due to parameter variation and the stochastic uncertainty in calculating initial values of emerging resistant virus [27,28]. Also, a feedback-free application of the schedule would be disastrous if unanticipated resistance to the permuted regimen is present, as this would result in the patient having uncontrolled virus replication for the duration of the schedule. The simplest method to avoid this is sample the viral load frequently following introduction of the permuted regimen, and switch to the naive regimen either when viral load reduction ceases or when a desired reduction in viral load is achieved.

One Previously failed therapy
The case where a patient has a single previously failed regimen had been the focus of our previous studies [29]. The only strategy that yields effective reduction of future mutation risk in this case involves total treatment interruptions.
Treatment interruptions and optimal scheduling. Our objective is to find a drug-switching schedule that yields the minimum risk, which is calculated as shown in the method in the Analysis section. For patients with only a single previously failed regimen, this can be achieved only through the use of interrupted schedules of treatment. The concept is identical to that driving the regimen cycling approach described in the previous section, except that periods of no treatment are allowed. If the resistant virus generated during the previous therapy has associated fitness cost with respect to the wild-type virus, then these periods of no treatment allow the wild-type virus to re-establish dominance. Reintroduction of the failing regimen will then result in a transient decrease in total viral load before the resistant strain re-establishes dominance.
As in the previous section, this is an optimal control problem in two steps. The optimal schedules consist of an interruption of length T 1 , to sensitize the virus, followed by the re-introduction of the failing regimen, resulting in transient suppression of the sensitized viral load. Before viral load rebound occurs, the naive antiviral regimen is introduced at time T 2 , resulting in a greatly reduced risk of subsequent virological failure. T 1 is again calculated to minimize the cost function of Equation 4, and T 2 is the switching time that achieves this minimum.
An example is used to illustrate how our algorithm works. The simplified model of Equation 2 is employed with the parameters which are identified for Patient 1. Assume that the point genetic distances v w {v e j j H~v r {v e j j H~2 . The optimal time for introduction of the naive regimen is fixed by the initial conditions at time T 1 . The optimal interruption length T 1 changes as a function of parameters; there may be a true optimum, or the optimal time may be infinite, as shown in Fig. 5 [25]. Where there is no minimum, the knee in the curve, where increased interruption time yields only marginal benefit, dictates the switching time. Fig. 6 shows the associated optimal switching schedule, with viral load as a function of time. In this case, the risk of resistance emergence is reduced from almost 55.4% to 0.18% by applying our algorithm.
Evolutionary Risks. Treatment interruption regimens have been associated with a high rate of resistance emergence [30], and are consequently avoided in standard HIV therapy [26]. Careful analysis of their use in this context is needed to avoid the possibility of encouraging development of multi-drug resistant viral strains.

Discussion
We have presented methods to reduce the risk of drug resistance emergence in HIV by manipulating the viral load prior to introduction of a naive antiviral regimen. These methods rely on creating a transient reduction in total viral load prior to introduction of a naive regimen. If the patient has failed multiple previous regimens, this may be accomplished either through regimen cycling or the use of a permuted antiviral regimen. If the patient has failed only a single previous regimen, this may be accomplished through the use of interrupted therapy. The optimal switching regimens are computed using simple model-based openloop optimal control algorithms. In all cases, the models predict achievable order-of-magnitude reduction in the risk of resistance to the naive regimen pre-existing its introduction.
The method proposed in this paper uses predictive models of HIV evolution under conditions of dynamic treatment to find treatment schedules that minimize the probability of certain resistant strains emerging. The application of dynamical systems and control to evolutionary systems will likely find broader application, as problems of drug resistance continue to increase in other areas of medicine.
While risk reduction should be achievable through any of these three methods, regimen cycling and interrupted therapy carry an increased risk of disease progression and/or further development of resistant virus. For these reasons, initial clinical investigations will focus primarily on the method of permuted regimens. Nevertheless, the other two methods may be useful in certain circumstances.
The problem of multi-drug resistance in HIV is extremely widespread, and methods that preserve remaining suppressive antiviral regimens have the potential to significant decrease morbidity and mortality in the HIV-positive population. The necessary next steps for implementing this method are outlined below:

Implementation Issues and Future Works
The methods introduced in this paper have the potential to significantly reduce the incidence of pre-existence related treatment failure. The methods involving permuted regimens could be applied to select groups of patients without significant further work. These patients would have to be in a closely monitored clinical situation, have a history of consistent viral genotyping showing strain patterns amenable to this method, and would have to be available for frequent viral load sampling. For these methods to be more broadly applied, several issues involving measurement, viral load history, and sampling frequency will need to be addressed.
Unmodeled Phenomenon. The model of Equation 2 represents a highly reduced form of the HIV dynamics, and neglects many known factors in HIV infection. Perhaps chief among these is the immune response to HIV, which can change the infection rate b and the death rate of infected cells a significantly. There is some previous work suggesting that the immune response can change dynamically with respect to changes in viral load, and that the consequences of this can be significant, potentially leading to immunological control of the virus under some circumstances [31][32][33]. However, multiple previous experiments summarized in [34] show that, once the chronic infection stage is reached, the immune system is permanently damaged and no longer displays such dramatic dynamic response  to changing viral loads. The excellent fits to the simplified model over several successive treatment interruptions shown in Fig. 2 argues strongly for the sufficiency of the model of Equation 2 to predict the dynamics of HIV during therapy changes. Furthermore, the algorithm proposed in this paper does not actually depend on the exact form of the model, but uses a closedloop sampling method to find the viral load minimum, and is therefore robust to unmodeled phenomenon such as a changing immune response. Measurement issues. The methods described in this paper rely heavily on the frequent and accurate measurement of the HIV viral load, both for risk computation and model identification. Viral load measurements are complicated by the existence of a detection limit of 50 copies/mL. Targeting viral load reduction below the limit of detection is therefore problematic. However, reduction of viral load to this limit of detection is sufficient in most cases to achieve a significant order-of-magnitude reduction in risk.
A second issue related to measurement is the proportion of the measured viral load which is noninfectious. Cells with defective integrated viral genomes may produce replication-incompetent virus particles, and a proportion of the particles produced by all infected cells will also be replication-incompetent. The most direct measurement of this phenomenon estimates that between 5-13% of the total free virus particles are capable of successfully completing the infection process (a much lower percentage actually complete the process, due to high inefficiencies at each intermediate step of replication) [35]. If a model of HIV dynamics is identified against plasma HIV RNA quantifications, which do not differentiate between infectious and noninfectious particles, then the estimate of the infection rate parameter b will implicitly be multiplied by the percentage of plasma virions which are infectious. This will not be a problem so long as all measurements are consistent between the model identification and prediction; estimated values of the number of infected cells will not be affected. The method of action of protease inhibitors results in an increase in the percentage of viruses which are noninfectious; this reduced fraction will be implicitly captured in the estimated value of the drug efficacy parameter j associated with the PI-containing regimen. The proportion of measured virus that is non-infectious will not affect the risk reduction algorithm, as the algorithm attempts to minimize the measured virus prior to switching, which is always proportional to the infectious virus. This phenomenon will slightly alter the calculated risk associated with a given measured viral load, but as the risk only depends logarithmically on the viral load, the change due to this relatively small proportional difference will be negligible.
A third issue, related to the second, is the fact that plasma viral concentrations are not the same as the virus concentrations in the lymphoid tissues, where the majority of the HIV virus resides and the where the majority of virus dynamics occur. HIV virus penetrates into many different tissues in the body, and there is evidence that these tissues are sufficiently compartmentalized to allow for divergent evolution of the virus in different compartments [18][19][20].The conversion factor developed in [21] and used in Equation 1 is a good first order approximation of the relationship between plasma virus level and total viral burden. A recent study in SHIV infected Rhesus Macaques has shown excellent proportional correlation between viral concentrations in plasma and various other tissue types both under treated and untreated conditions [36], indicating that the plasma virus load is a good surrogate measurement of total viral burden, even under conditions of dynamic therapy. The exact ratio between plasma viral load and total viral burden will change from patient to patient, but this will not affect the proposed algorithm, though it will slightly alter the calculated risk associated with a measured viral load, as discussed above.
Choosing drug permutations. The most promising method presented in this paper is the introduction of permuted antiviral regimens prior to the introduction of a naive regimen. In order to safely choose these permuted regimens, it is necessary to know which resistant viruses are present in the patient's viral reservoirs. Consistent viral genotyping following every treatment failure would provide this information; unfortunately, this is rare. It may be possible to estimate the likely distribution of resistant viruses in a patient based on a history of antiviral use and failure, using genetic distance and fitness information from the HIV database.
Finding Minima. All the methods presented in this paper induce a transient crash in the viral load, and depend on being able to switch therapies at or near the minimum of this crash, prior to rebound. In this paper, the method suggested for this is consistent, frequent sampling of the viral load during the transient period. While this should work, it is expensive both in economic  terms and in terms of patient burden. Optimal minimal sampling methods to find the viral load minima with the fewest possible measurements should solve this issue; this research is ongoing [37][38][39].

Author Contributions
Conceived and designed the experiments: RZ RL. Performed the experiments: RL MJP RZ. Analyzed the data: RL MJP JMP RZ. Contributed reagents/materials/analysis tools: JMP. Wrote the paper: RL MJP JMP RZ.