Drug-Class Specific Impact of Antivirals on the Reproductive Capacity of HIV

Predictive markers linking drug efficacy to clinical outcome are a key component in the drug discovery and development process. In HIV infection, two different measures, viral load decay and phenotypic assays, are used to assess drug efficacy in vivo and in vitro. For the newly introduced class of integrase inhibitors, a huge discrepancy between these two measures of efficacy was observed. Hence, a thorough understanding of the relation between these two measures of drug efficacy is imperative for guiding future drug discovery and development activities in HIV. In this article, we developed a novel viral dynamics model, which allows for a mechanistic integration of the mode of action of all approved drugs and drugs in late clinical trials. Subsequently, we established a link between in vivo and in vitro measures of drug efficacy, and extract important determinants of drug efficacy in vivo. The analysis is based on a new quantity—the reproductive capacity—that represents in mathematical terms the in vivo analog of the read-out of a phenotypic assay. Our results suggest a drug-class specific impact of antivirals on the total amount of viral replication. Moreover, we showed that the (drug-)target half life, dominated by immune-system related clearance processes, is a key characteristic that affects both the emergence of resistance as well as the in vitro–in vivo correlation of efficacy measures in HIV treatment. We found that protease- and maturation inhibitors, due to their target half-life, decrease the total amount of viral replication and the emergence of resistance most efficiently.


Introduction
Since 1996, human immunodeficiency virus (HIV) infection is treated with a combination therapy, known as highly active anti-retroviral therapy (HAART) [1,2], which has substantially improved the clinical management of HIV [3]. Despite the success of HAART, eradication of HIV can currently not be achieved [4,5], most likely due to the persistence of virus in very long lived, latently infected cells [6,7]. For HIV-infected individuals, life-long therapy is therefore required to prevent progression to the acquired immunodeficiency syndrome (AIDS) and death.
During therapy, plasma viral load (HIV RNA per mL blood plasma) is recommended by the National Institute of Health as a marker of therapy success [8], whereas measurement of the CD4 cell count is the most important clinical marker of disease progression [9]. The in vivo potency of novel antivirals is usually assessed by viral load decline in small clinical trials of monotherapy, e.g., [10,11], and later evaluated utilizing the novel agent in combination with an optimized background therapy, e.g., [12]. The in vitro potency of antivirals is typically assessed by using phenotypic/single-round infectivity assays [13][14][15][16], which measure the number of offspring after one round of virus replication.
Investigation of novel drug targets for the treatment of HIV infection resulted in the development of new drug classes. In 2003 and 2007, the fusion inhibitor (FI) enfuvirtide [17], the CCR5antagonist maraviroc [18] and the integrase inhibitor raltegravir [19] were approved for the treatment of HIV infection. Many more drugs are in late clinical development [20]. With the introduction of new drug classes, in particular integrase inhibitors, a huge discrepancy between the efficacy measured in vitro, using phenotypic/single-round infectivity assays, and in vivo, using viral load decline, was observed [14,21]. Although integrase inhibitors cause a steep initial decline of plasma viral load [21][22][23][24][25][26], the in vitro efficacy is amongst the lowest [14].
Mathematical modelling of viral dynamics has lead to many insights into the pathogenesis and treatment of HIV. It is a valuable tool to interpret the time course of virological markers (e.g. viral load) during HIV treatment [27][28][29][30][31] and contributes much to our current understanding of the in vivo dynamics of HIV. Sedaghat et al. [32,33] used a mathematical modelling approach to analyze the rapid decay of plasma viral load after application of integrase inhibitors. They infer that this characteristic viral decay is a result of the inhibited stage within the viral life cycle rather than superior in vivo potency.
Consequently, viral load decay may be misleading for assessing the potency of integrase inhibitors (and other novel inhibitors) in comparison to existing drug classes. However, an alternative, more appropriate measure of drug efficacy, which allows to directly compare drugs from different drug classes is still missing.
The objectives of this article are (i) to develop a novel, generic measure of drug potency that facilitates comparison across different drug classes; (ii) to develop a novel mathematical model of the viral replication cycle that incorporates the action of established and novel drugs in a mechanistic way; and (iii) to analyze determinants of drug efficacy critical for drug discovery and development. The proposed measure of drug efficacy, termed reproductive capacity, extends the established in vivo marker, plasma viral load, by incorporating additional infectious viral stages, and the in vitro phenotypic/single-round infectivity assays by taking into account host specific defense mechanisms. This enables us to understand the observed discrepancies between in vitro and in vivo efficacy for integrase inhibitors, and to elucidate and quantify the role of immune-system related clearance mechanisms in drug action. The results presented herein are of particular value to categorize different molecular targets in the HIV life cycle and are expected to be of significance for guiding future HIV drug discovery and development.

Development of a detailed model of viral life cycle and action of anti-retroviral drugs
We derived a detailed virus-target cell interaction model as depicted in Fig. 1. The model incorporates the mechanisms of action of all currently approved drugs and some drugs in late clinical development.
Target cells are produced by the immune system with some constant rate l T . An infectious virus V I reversibly binds (with effective rate constants k on and k off ) to a target cell TU, forming a complex V I : TU. After binding, the virus irreversibly fuses (with rate constant k fus ) with the target cell and the viral capsid containing the viral genomic RNA is released; this state is denoted by T RNA . During reverse transcription (with effective rate constant k rev ), genomic viral RNA is irreversibly transformed into a more stable DNA. Viral DNA and viral proteins form the pre-integration complex (PIC), denoted by T 1 . In the next step, viral DNA of the PIC is irreversibly integrated into the DNA of the target cell (with rate constant k T ), forming the provirus T 2 . After integration, the infected cell cannot return to an uninfected stage. From the proviral DNA, viral proteins are amplified and new viruses are released (with effective rate constant b N N T ½1=(cells : day)). Only a given percentage pw0 of the released viruses are correctly assembled immature viruses V IM , while the remaining percentage (1{p) are defective virions V D that might e.g. lack the (gag-pol-polyprotein contained) enzymes. During the final step, the viral protease, which is packed into the correctly assembled, immature virions V IM , is responsible for the maturation of the virus.

Author Summary
To guide drug discovery and development, measures of drug efficacy that are linked to clinical outcome are of key importance. In HIV treatment, decay of plasma viral load is typically used as an in vivo measure of drug efficacy, whereas phenotypic assays are used to assess drug efficacy in vitro. The recent development of novel HIV drugs resulted in a huge discrepancy between viral load decay and in vitro predictions of drug efficacy. We used a mathematical modelling approach to resolve this discrepancy by introducing a new quantity, the reproductive capacity, that allows a transfer of the in vitro drug efficacy measure into the in vivo context, enabling a direct comparison. We developed a novel model of viral dynamics that incorporates the mechanism of action of all established and novel antivirals. Based on the model, we analyzed the ability of the viral infection to replicate under different drug treatments, and estimated classspecific times until virological failure. We conclude that the half life of the targeted viral stage is an important classspecific attribute that impacts on the overall success of a drug in vivo. Our findings have direct implication for the drug discovery and development process. The maturation of HIV virions has been shown to be dependent on the highly ordered cascade of cleavages, governed by differences in the inherent processing rates at each cleavage site [34,35]. We assume that a fraction (1{q) of the released virus matures abnormally, contributing to the pool of defective virions V D . Successful maturation eventually leads to new infectious virus particles V I (with rate constant k mat and probability q).
Depending on the stage of the life cycle, the host organism has different abilities to clear the virus. It was assumed that infectious, immature and defective virions V I , V IM , and V D , respectively, are cleared with rate constant CL by the host. The uninfected target cells TU, the T RNA stage and the early infected stage T 1 are assumed to be cleared with rate constant d T , since none of these stages express viral proteins, while the virus-producing late infected cell T 2 is assumed to be cleared with rate constant d T2 &d T . In addition to cell death, the target cell may fend-off the viral infection by degrading the viral RNA or parts of the PIC, rendering the cell uninfected. RNA is very unstable with a half life ranging from seconds to a maximum of two hours [36,37]. Therefore, through degradation or, e.g., by hypermutation through APOBEC3G [38], the viral RNA can be cleared with rate constant d RNA . The cell might also destroy essential components of the PIC (with rate constant d PIC,T ) to clear the virus.
The system of ordinary differential equations (ODEs) describing the rate of change of the different viral species and target cells in the detailed model (depicted in Fig. 1) is given in Supplementary Text S1, Eqs. (S1)-(S8). As typically done in kinetic studies, complex aspects of the viral dynamics are subsumed by 'lumped' parameters in the model. For instance, the rate constant of the reverse transcription k rev contains all the steps necessary to transform the viral RNA into a double stranded DNA. The mechanisms of action of the seven drug classes are based on interfering with the viral life cycle at different stages. We assumed that the effect of a drug on the targeted process is specified by some parameter e(t)[½0,1, i.e., assuming some underlying averaged drug concentration C~b C C, see [39], some fifty percent inhibitory concentration IC 50 , and some drug specific Hill coefficient n, see [14]. For the purpose of the study, this rough approximation is sufficient, however, it is possible to also use time-varying drug concentration C~C(t) resulting from some pharmacokinetic model, or to use more mechanistic effects models [40,41]. The actions of the different drug classes within the viral life cycle are shown in Fig. 1. CCR5 antagonists inhibit the association of HIV with the CCR5 receptor in CCR5-tropic virus. They thus affect the association constant k on . Fusion inhibitors (FI) inhibit the process of HIV fusion, affecting k fus . Activated nucleoside reverse transcriptase inhibitors (NRTI) compete with endogenous deoxynucleoside triphosphates for prolongation of the growing DNA chain, while non-nucleoside reverse transcriptase inhibitors (NNRTI) allosterically inhibit the function of the reverse transcriptase. The effects of both drug classes result in a reduced rate at which the RNA is reversely transcribed into DNA. Integration inhibitors affect the integration of viral DNA into the host genome catalytically [42][43][44][45]. In the proposed model, this alters the transition rate constant k T from early infected cells T 1 to the late infected cells T 2 . Protease inhibitors (PI) bind to the catalytic pocket of the viral protease enzyme, which is responsible for the processing of the viral precursor polyproteins and thus the maturation of viral particles. In the proposed model ( Fig. 1), PIs therefore inhibit maturation by decreasing the maturation constant k mat . Maturation inhibitors (MI) bind to the substrate of the viral protease (Gagpolyprotein) [46] at a specific site. This binding perturbs the ordered sequence of cleavages that is necessary for proper maturation [47,48], resulting in defective virus morphology [49]. In the proposed model (Fig. 1), MIs therefore decrease the probability q that immature virus matures normally, increasing the proportion of abnormally matured, defective viruses V D .

Impact of antiviral drugs on relative abundance of infectious viral stages
We used the detailed virus-target cell interaction model to predict the effect of the different drug classes on the distinct stages of the viral life cycle. In order to enable a direct comparison between the different drug classes, we artificially eliminated the feedback by keeping the uninfected target cell TU and the infective virions V I that 'enter' the infection cycle constant (the two leftmost species in Fig. 1), resulting in 'downstream' quasi-steady state numbers T 1,ss , T 2,ss , V IM,ss , V I,ss , and V D,ss . For a given drug class and inhibition of the targeted molecular process e, the effect of the drug on the life cycle was quantified by the four ratios as shown in Fig. 2. As expected, the drugs perturb the ratios of viral states that encompass their site of action within the viral life cycle. In the present example, all states that lie downstream of the drugs' target site are affected, while the states that lie upstream are usually not affected. The exception are InIs, which increase the abundance of the preceding stage T 1 ( Fig. 2A), while decreasing the number of the subsequent infectious stage T 2 (Fig. 2B). Interestingly, the effect on the ratios is not always a linear function of drug efficacy. PIs and MIs also show a different behavior ( Fig. 2D): PIs affect the ratio of infectious-to-defective virions by decreasing the maturation rate k mat , which lowers the number of infective virions V I , but also lowers the number of virions that mature abnormally (contributing to V D ). MIs increase the proportion of virus that matures abnormally and decrease the proportion of virus that matures normally, thus decreasing V I and increasing V D , without affecting k mat .

Development of a simplified two stage virus dynamics model
The detailed model ( Fig. 1) contains parameters that are difficult to measure and currently not available. We therefore reduced the detailed model based on reasonable quasi-steady state assumptions to obtain a simplified model of virus-target cell interaction dynamics that is parameterizable in terms of established and validated parameter values (see Supplementary Text S1). In particular, we have eliminated the intermediate stages of the cellvirus complex TU : V I , the infected cells prior to reverse transcription T RNA and the immature virus V IM in the original model ( Fig. 1). As a consequence, we derived a lumped infection rate constant b, which describes the infection of a susceptible cell towards the stage, where the viral RNA has been successfully transformed into DNA. We also derived a virus clearance CL T that is associated with the loss of virus during the intermediate stages before reverse transcription and the release rate constant of infectious virus N.
The infection rate constant is given by where k fus denotes the fusion rate constant, K D the dissociation constant of the virus-target cell complex, and r rev,w denotes the probability that reverse transcription is successfully completed (see Supplementary Text S1). The lumped virus clearance (loss of virus by, e.g., genome destruction) in the intermediate stages is given by the parameter The number of released, infectious viruses is given by where p and q are the probabilities that the released virus is correctly assembled and matures normally, and r PR,w is the probability that the released virus matures before being cleared by the immune system (see Supplementary Text S1). The lumped model can be parameterized in terms of six unknown parameters (b, b N N,l T ,d T ,d T 2 ,CL), which equals the number of estimated parameters using standard models [28]. For the remaining parameters, we have provided values from the literature (see Supplementary Text S1).
In the following, we considered two types of target cells (T-Cells and a longer lived cell population, which we refer to as macrophages) and finally incorporated the viral mutation process (resulting from erroneous reverse transcription) into the overall model. Whether the longer lived cell population consists solely of macrophages in vivo remains unknown. There is, however, some evidence that the kinetic characteristics of the longer lived cell population are similar to those of the macrophage population [33]. The proposed simplified twostage virus dynamics model is shown in Fig. 3. It comprises T-cells, macrophages, free non-infectious virus (TU,MU,V NI , respectively), free infectious virus of mutant strain i,V I (i), and four types of infected cells belonging to mutant strain i: infected T-cells and macrophages prior to proviral genomic integration (T 1 (i) and M 1 (i), respectively) and infected T-cells and macrophages after proviral genomic integration (T 2 (i) and M 2 (i), respectively). The rates of change of the different species in the reduced two-stage HIV model are given by the following system of ODEs:  (5)). The parameters CL T (i) and CL M (i) denote the clearance of mutant virus i through unsuccessful infection of T-cells and macrophages respectively (see Eq. (4)) and the parameters b T (i) and b M (i) denote the successful infection rate constants of mutant virus i for T-cells and macrophages respectively. The parameter p k?i denotes the probability to mutate from strain k to strain i (to be defined below).
The model enabled us to mechanistically incorporate the action of all drugs that are approved or in late clinical trial. The impact of a compound on a corresponding (lumped) parameter in the model is specified by g: The same quantities are defined for macrophages by replacing the subscript T by M; see Supplementary Text S1 for details. The overall viral dynamics model comprises a complete mutagenic graph. In HIV infection, genomic mutation occurs during the reverse transcription process [50]. The reverse transcriptase of HIV lacks a proof reading mechanism in contrast to host polymerase enzymatic reactions. However, viral proteins from newly mutated viral genomes are only produced after integration of the viral genome into the host cell DNA. The proteins required for the stable integration of the newly mutated viral genome originate from the founder virus. Therefore, phenotypically, drug resistance of new mutants will only be observed after integration, i.e., in the infectious stages T 2 and M 2 . In total, the model includes 2 L different viral strains i that contain point mutations in any pattern of the modelled L possible mutations. For two distinct mutations L~2, the mutagenic graph is shown in Fig. 4A. Each mutant i can mutate into every other mutant k in one step. The probability p k?i to mutate from a strain k into another strain i can be directly derived from the mutagenic pathways in Fig. 4A, i.e., where m denotes the mutation probability per base and reverse transcription process (m&2:16 : 10 {5 [50]), h(i,k) denotes the hamming distance between strain k and strain i, and L is the total number of different positions that are considered in our model. The phenotype of each mutant strain i is modelled by introducing a selective disadvantage s(i), which denotes the loss of functionality (e.g., in the activity of some viral enzyme that is affected by the mutation) relative to the wild type, and a strain specific inhibitory activity (g(i,j)) of treatment j against the mutant strain i. For example, the strain specific infection rate i under a certain treatment j is given by denotes the infection rate constant of the wild type wt in the absence of drug w (given in Table 1). Since some viral strains are present only in very low copy numbers, we used a hybrid stochastic deterministic approach [51] to model the overall virus dynamics model (see Materials and Methods section for details).

Reproductive capacity for predicting drug-specific impact on viral replication
The production of infectious offspring is crucial for the survival of a viral population. The phenotypic single-round infectivity assay measures the amount of infectious offspring after one round of replication. For a given drug, the assay quantifies the drug's efficacy by measuring the reduction in viral offspring relative to the drugfree situation. We defined a new quantity-termed the reproductive capacity R cap -, which transfers the principle of the phenotypic single-round infectivity assay into a mathematical term. Its definition involves the quasi-species distribution and the basic reproductive numbers of all pathogenic sub-stages. The reproductive capacity characterizes the fitness of a given state of the infection from the perspective of a potential treatment j by quantifying the expected total number of offspring under the treatment j.
The basic reproductive number R 0 is a well characterized quantity in epidemiology that denotes the expected number of  secondary infections caused by a single infected cell/virus [52]. If R 0 w1 then the infection will spread, while for R 0 v1 the infection will die out. The strain associated reproductive number R 0 (i,j) characterizes the fitness of a viral strain i in a pharmacologically modified environment, specified by a drug treatment j. We used the 'survival function' approach [53] to calculate the reproductive numbers for mutant strains i under treatment j. In our context, the survival function is of particular value, since it captures the possible event of mutation for all infective classes. Based on the two-stage virus dynamics model, the basic reproductive number R V (i,j) of a single virus of strain i under treatment j is given by with constants Since infected cells are also pathogens, which can lead to a rebound of the disease even in the absence of any virus, we also determined their basic reproductive numbers under a given treatment j. The basic reproductive numbers R T 1 (i,j) and R M 1 (i,j) of the infectious stages T 1 and M 1 , associated with the viral strain i, are given by Finally, the reproductive numbers R T 2 (i,j) and R M 2 (i,j) of the infectious stages T 2 and M 2 , associated with the viral strain i, are given by We defined the reproductive capacity R cap (j) of the entire quasispecies ensemble under treatment j as the weighted sum of the basic reproductive numbers of all pathogenic stages of mutant strain i, i.e., free virus, infected T-cells and infected macrophages, where the weights are the abundance of the corresponding pathogenic stage: The reproductive capacity R cap (j) can be interpreted as the expected total number of infectious offspring that the infection produces in one round of replication under a certain treatment j, given the current state of the infection.
Relation to viral load and phenotypic/single-round infectivity assay. The viral load considers the total concentration of free virus, consisting of non-infectious virus V NI and infectious virus V I (i), belonging to all mutant strains i. In contrast to the reproductive capacity, viral load does not assess the ability of distinct viral strains i to replicate (in terms of R V ). In mathematical terms, the viral load is given by The in vitro reproductive capacity, corresponding to the read-out of the phenotypic assay R pA (j) (under treatment j) is conceptionally similar to Eq. (20). However, in comparison to the above defined in vivo measure, the in vitro measure would not take into account: (i) the clearance of any infective stage by the immune system (relating to the parameters CL,CL T (i,j),CL M (i,j),d T ,d M ,d T 2 , and d M 2 ), and (ii) the abundance of the different infected cell types (e.g., Tcells and macrophages). The assay measures one round of replication, denoted byR RT T 2 , starting from a late stage infected cellT T 2 . Mathematically expressed, the primary output is given by Drug-class specific decay of viral load and reproductive capacity Application of drugs/drug classes changes the total size and the composition of the viral population. The impact of this change is typically evaluated in terms of the decay of the viral load over time. We used the reproductive capacity R cap (j) to also evaluate viral replication under various hypothetical treatments j. In Fig. 5, we predicted the impact of the different drug classes on the decay of the plasma viral load and the reproductive capacity R cap (w), i.e., the fitness of the whole virus population, evaluated in the absence of drugs. As typically done, we assumed 100% drug efficacy g.
In terms of the plasma viral load decay (Fig. 5A), we observe a faster initial decay for InIs in comparison to all other compound classes, in agreement with clinical data [21] and theoretical analysis [32,33]. The onset of viral load decay is delayed for all other compound classes as observed clinically [12,27], see also Figure S1. In agreement with clinical data [21], in the case of InI treatment, the second phase of viral decay starts earlier after treatment initiation and exhibits &70% less viremia in comparison to other drug classes, but shows the same decay. Notably, the change of the ratio of infective virus-to-total virus (see Fig. 5, inset) upon PI or MI administration is not reflected by the total viral decay in the blood plasma.
Most noticeable, the reproductive capacity (Fig. 5B) discriminates between RTIs, FIs and CCR5-antagonists vs. InI vs. PIs and MIs. It can be seen, that protease and maturation inhibitors reduce R cap most efficiently initially and shift it to an overall lower level. Integrase inhibitors cause a slightly faster initial decay in R cap , in comparison to RTIs, FIs and CCR5-antagonists, which consistent with the rapid decay in viral load (Fig. 5A). However, in contrast to viral load decay, the initial fast decay of R cap levels off and the second phase decay is flatter for InIs in comparison to RTIs, FIs, CCR5-antagonists, PIs and MIs. The effect of NRTIs, NNRTIs, CCR5 inhibitors and FIs on R cap is comparable (Fig. 5B). Remarkably, these inhibitors induce an initial increase in R cap (see next section for details), followed by a slow first phase decay, followed by a second phase decay that is parallel to the decay of R cap in the case of PI-and MI-treatment, sustaining overall higher levels of R cap in comparison to PIs and MIs. In the next section, we further elucidate the reasons for these classspecific differences.

Immune-system related clearance is critical determinant of drug-class specific decay
In view of the analysis performed in Fig. 5B, R cap is directly correlated to the overall abundance of viral infectives (V I ,T 1 ,T 2 ,M 1 ,M 2 ).
PIs and MIs primarily act on infectious virus V I (see Fig. 5, inset), by reducing the proportionality factor (N T =CL and N M =CL) that determines the abundance of infectious virus in the first-and second phase decay (see Eq. (10)). The infectious virus V I is rapidly cleared by the immune system [54]. Therefore, application of highly efficient PIs and MIs leads to a rapid reduction of infectious virus V I , as illustrated in Fig. 6D and Fig. 5 (inset). This reduction is also reflected by the initial drop of R cap in Fig. 5B. In the case of PI and MI treatment, infected T-cells are quickly becoming the most abundant infectious compartment (Fig. 6D) and subsequently dominate the decay characteristics of R cap in Fig. 5B. In the final phase, late infected macrophages (M 2 ) are becoming the most abundant compartment and thus dominate the decay of R cap in the final phase.
Integrase inhibitors prevent the integration of the viral genome and thus prevent the transition of early infected cells (pre-integration, T 1 and M 1 ) to late infected cells (post-integration, T 2 and M 2 ), see Fig. 3. By inhibiting the transition from early to late infectious cells, integrase inhibitors increase the decay of late infected T 2 -cells (see Fig. 6C). In the case of InI treatment, infectious virus V I is initially proportional to T 2 , explaining the observed more rapid first-phase decline in R cap in Fig. 5B. However, blocking the transition from T 1 to T 2 can also slow the decay of the T 1 -compartment, which might become more abundant than V I after the initial decay. In the final phase both T 1 and V I become proportional to-and remain more abundant than M 2 , which explains the overall higher levels of R cap in the final phase (see Fig. 5B).
The effects of NRTIs, NNRTIs, CCR5 inhibitors and FIs on R cap are comparable (Fig. 5B), as they primarily act on preintegrative early infected cells (T 1 and M 1 ). The difference between entry inhibitors and reverse transcriptase inhibitors is marginal, because the clearance of virus by infection is negligible compared to the clearance by the immune system (CL T vCL and CL M vCL). A positive result of entry inhibitors (FI/CCR5) and RTIs (NRTIs/ NNRTIs) is an increased number of uninfected cells, which also results in an initial increase in the reproductive capacity R cap (see Fig. 5B). During treatment with NRTIs, NNRTIs, CCR5 inhibitors and FIs, infective virus V I is the most abundant compartment. The decay in the first phase is proportional to the decay of the late infected cells, T 2 . Once the abundance of T 2 falls below M 2 , the decay of V I and thus R cap in Fig. 5B is proportional to the decay of late infected macrophages M 2 .

The pattern of virological removal influences the time to virological rebound after treatment application
In the following, we predict how the distinct viral dynamics after drug application affect drug efficacy in vivo. The long-term in vivo efficacy of an antiviral drug depends on many different factors,  Table 1. Black diamonds indicate median viral load data from [27] (PI monotherapy), numerically available in [70]. Black squares and black bullets indicate median viral load data from [21] (NRTI + background therapy and InI+background therapy, respectively). The horizontal dashed black line indicates the limit of detection of current assays (50 copies of HIV RNA per mL). Inset: Protease-and maturation inhibitors (PI and MI) change the ratio of infectious to total virus (V I : V tot ). B: The evolution of the reproductive capacity (evaluated at the drug free state R cap (w)) after treatment with different drug classes. Model parameters are as indicated in Table 1. The initial infection was assumed to consist of wild type only. Drug efficacy g was assumed to be 100%. Total body virus has been converted to plasma viral load by assuming that the virus distributes into the plasma (V plas:~3 :1 liters, which surrounds 2% of infected cells) and the interstitial space (V int:~9 :6 liters [71], which surrounds 98% of infected cells). The volume of distribution with reference to the plasma concentration has been calculated using the well-known formula vol. distr~K int::plas: : V int: zV plas: , see e.g. [72], where K int::plas:~9 8%=2%~50. Finally, we assume that on average each virus contains 2 viral RNAs (which are measured [viral RNA/mL] plasma). doi:10.1371/journal.pcbi.1000720.g005 including the ability of the virus to adapt to the pharmacological challenge by developing resistance mutations. The ability to develop drug resistance is strongly dependent on the induced pattern of resistance mutations against a particular drug, but also on the velocity at which replication competent compartments (V I ,T 1 ,T 2 ,M 1 ,M 2 ) are removed from the body. Since antiretroviral drug classes target different stages in the viral life cycle, they are likely to induce different patterns by which viral compartments are removed from the body (see Fig. 6) and might therefore exhibit different long-term in vivo efficacies.
To illustrate the sole impact of virological removal (V I ,T 1 ,T 2 ,M 1 ,M 2 ) on resistance development and therefore on drug efficacy, we have intentionally chosen a simplistic, unified mutational landscape and considered the time to viral rebound as a long-term measure of efficacy. We denoted virological rebound, if the viral load reaches 90% of the pre-treatment viral load. We assumed that the drugs inhibited their targeted (lumped) parameter (see Eqs. (7)-(10)) by 90% in the wild type (g~0:9), by 45% in a one-mutation strain (g~0:45) and are entirely inefficient in the double-mutant (g~0). Drug-specific and more realistic mutational landscapes are possible, but in view of the current analysis (elucidating the impact of class-specific virological removal), they would blur the results.
In Table 2, the time to virological rebound for the different drug classes based on the above simplistic mutation model is reported.
The virus generally rebounds to 90% of pre-treatment levels after 1-2 month of monotherapy, which is in the same order of magnitude as clinically observed rebound times [55][56][57]. Although inhibition g was assumed to be identical across all drug classes, the  The time to virological rebound depends on both the cost of resistance ('selective disadvantage', s) and the choice of drugs. Each table entry shows the time to virological rebound in [days] in an ensemble of 1000 hybrid stochastic deterministic simulations, where we assumed that the efficacy of the drugs against the wild type was 90%. The drug was 45% effective against an onemutation strain and completely inefficient against the double-mutant. The fraction of non-infectious viruses (1{p : q : r PR,w ) was set to one-third and the initial population was assumed to be all wild type. The viral load was said to be rebounded, if the viral load reached 90% of the pre-treatment viral load. doi:10.1371/journal.pcbi.1000720.t002 times to virological rebound differed. In particular, when resistance confers a marked loss in fitness (i.e. selective disadvantage = 30%), PIs show the longest time to virologically rebound, and the InIs the shortest. For integrase inhibitors, the difference between the decay of plasma viral load and their predicted long-term efficacy is quite pronounced. Their comparably shorter times to virological rebound are in strong contrast to their steep initial decrease of plasma viral load (see Fig. 5A), but consistent with the decay pattern of the reproductive capacity (Fig. 5B). For the EIs, RTIs, PIs and MIs, the predicted time to virological rebound is also much more consistent with the decay characteristics of the reproductive capacity (Fig. 5B) than with the decay pattern of total viral load (see Fig. 5A).

Discussion
In clinical studies, the first approved integrase inhibitor, raltegravir, induced an extremely rapid decline in viral load when applied both as monotherapy [10] and in combination with an optimized NRTI background therapy [21][22][23][24]. While it was initially speculated that the observed decline might be a result of superior potency of raltegravir, it is now emerging that the viral decline in InI-based therapy could be a class-specific phenomenon [25,26]. Moreover, superior potency of InIs (in terms of g) was not confirmed by single-round infectivity assays [14]. The mechanisms underlying the decay dynamics are still not clear [58] and controversially discussed [21,32].
In [32], a two stage model of the viral replication cycle is presented, which explains the differences between the decay of viral load between RTIs and InIs based on the stage at which the drugs affect the dynamics of viral replication. The model explicitly distinguishes two viral stages, early-stage infected cells and late-stage infected target cells, which are specifically defined for a pair of drugs under examination. The authors further conclude that the viral dynamics produced by drugs from different anti-retroviral classes should not be directly compared to infer drug potency [33]. An alternative measure, as it is imperative for guiding drug discovery and prioritizing drug candidates in later development stages, is still lacking.
All currently approved antivirals exert their effect by inhibiting the replication of HIV. The extent at which replication is inhibited, is therefore a unifying indicator for drug efficacy across all drug classes. Replication assays, e.g., phenotypic assays [15] or replication capacity assays [59], analyze drug efficacy in terms of viral replication in vitro. The replicative fitness of HIV in vivo, however, depends on the interaction of a multitude of viral and host factors. Replication assays represent the dynamics of HIV under the assay conditions, which lack many host factors, in particular the immune responses to the infection. However, since it is particularly useful to compare compounds in terms of replication inhibition, we adopt the dynamic approach of replication assays to define the reproductive capacity R cap . In silico, we are able to consider the host response to the viral infection and can thus extrapolate the replication approach from in vitro to in vivo. In [60], the authors used a similar approach to compare the effect of distinct antiviral classes utilizing age-structured models.
We derived a single detailed model of the viral replication cycle and deduced a reduced two stage model, which incorporates the action of all approved HIV drugs. Our two-stage model allows to predict the action of any number of drugs simultaneously, including common HAART cocktails, potentially belonging to different drug classes. In contrast, in [32], the stages of the two-stage model of viral replication are not specified a priori and have to be determined by the two drugs that are analyzed and compared.
Based on the proposed detailed and reduced model, we identify the following effects of currently approved drugs: EI and RTIs decrease the infection rate and thus the number of new infections. The impact on the release of new virus (and virus decline) is therefore delayed by the viral life cycle. MIs and PIs do not interfere with the total amount of virus that is being released, but rather shift the ratio of infective to total virus, V I : V tot (see Fig. 5, inset), which is not directly reflected by total plasma viral load. Since the kinetics of the free virus are rapid [54], this has an immediate impact on the number of new infections. Subsequently, this impact on the number of new infections affects the total viral release (and thus total plasma virus load) in a similar manner as EIs and RTIs, creating a 'shoulder' phase. Hence, we obtain In our model, EIs, RTIs, PIs and MIs produce an identical decay of plasma viral load (see Fig. 5A), when assuming 100% inhibition, respectively. In particular, the onset of viral load decay is similarly delayed ('shoulder phase') with these inhibitors (see Figure S1), in agreement with clinical observations [12,27]. Previously discussed theoretical differences in the viral response between RTIs and PIs (see Eq. (5.7) vs. Eq. (5.16) in [61]) yield similar dynamics when more recent (higher) estimates of viral clearance are used [54].
In contrast to other inhibitor classes, InIs decrease the amount of late infected cells (T 2 ,M 2 ) (see Fig. 2), which has an immediate impact on total virus release, i.e., total virus release~b N N : The impact of InIs on viral load decay is immediate and not delayed by the viral replication cycle as in the case of all other compounds [12,27]. Thus, the onset of observed total viral decay is faster for InIs than for other compounds, irrespective of their potency (which was set equal for all compounds in Fig. 5A). Furthermore, the decay of viral load in the first phase is steeper for InIs in comparison to other inhibitor classes (see Fig. 5A). The viral load decline in the first phase is proportional to the decay of the late infected T-cells T 2 (see Fig. 6). Sedaghat et al. [32] derived analytical solutions for the viral decay dynamics after InI and RTI treatment (see Eqs. (9) and (10) in [32]), which demonstrate that the viral decay after InI treatment is determined by the death rate of late infected cells (d T 2 ), while in the case of RTI treatment, the decay is determined by the ''flushing-out'' of the early infected cells (T 1 ) and the death rate of the late infected cells d T 2 , leading to overall faster viral declines in the case of InI treatment in the first phase. The long-term in vivo efficacy of an antiviral drug depends on many different factors, particularly the ability of the virus to adapt to the pharmacological challenge by developing resistance mutations. The ability to develop drug resistance is strongly dependent on the induced pattern of resistance mutations against a particular drug, but might also be influenced by the velocity at which replication competent compartments are removed from the body. However, viral load decay focusses on only one single variable, namely the total output of virus, whereas other infectious stages (e.g. T 1 ,T 2 ,M 1 ,M 2 ) remain 'hidden'. Furthermore, the ratio of infective virus-to-total virus (V I =V tot ) is not resolved, which might underestimate the long-term efficacy of PIs and MIs that target this ratio (see Table 2 in relation to Fig. 5A). In the section 'The pattern of virological removal influences the time to virological rebound after treatment application', we have compared the impact of drug-class specific removal patterns on the long-term efficacy of antivirals (in terms of resistance development). We showed that although inhibition g was assumed to be identical across all drug classes, the times to virological rebound (used as a measure of long-term efficacy) differed, with PIs showing the longest time to virologically rebound, and InIs the shortest.
The reproductive capacity has been monitored over time in Fig. 5B to depict class-specific long-term efficacy of antivirals based on the hosts' ability to clear the targeted infectant in the viral life cycle. The main conclusion is that the long-term efficacy is larger for compounds that target viral life-stages that are cleared at a fast rate. It is generally assumed that the free virus is cleared at the fastest rate [27,54]. Since MIs and PIs reduce the production of infective virus V I (see Fig. 2), they reduce the virus' ability to produce offspring faster than all other drug classes. Furthermore, since resistance development is correlated with the extent of replication, we infer that PIs and MIs, based on their viral target, are the most efficient drug classes in terms of reducing the probability of resistance development. This assumption correlates well with the observed rebound times in Table 2 and is also supported by the fact that the introduction of PIs marked the success of HAART [1].
During drug discovery, the pre-clinical-and the clinical development process, in vitro surrogate measures or in vivo drug efficacy measures are important to prioritize drug candidates.
The mechanistic mode of action of a compound at its target site can be elucidated by cell free assays that use purified viral protein, e.g. reverse transcriptase for RTIs. The influence of viral mutation, the immune system and pharmacokinetics are absent in this type of assay. However, it is possible to deduce the pharmacodynamic mode (e.g. Eq. (1), see also [41]) and thus the parameter e from these types of assays, which denotes the extent of inhibition of the molecular process by the compound. Mathematical models of HIV dynamics use a minimal number of parameters, making them suitable for parameter fitting and comparison with clinical data. The parameters used in the models are often lumped, summarizing many viral processes. For example, binding, fusion and reverse transcription are part of the infection rate b (see Eq. (3)). Inhibition of lumped model parameters (denoted by g) might therefore differ from inhibition of the molecular process e, which is measured by cell-free in vitro assays. We have provided equations (Eqs. (S24) and (S31), Supplementary Text S1) that enable the use of pharmacodynamic information e, derived from cell free assays (inhibition of the targeted molecular process), in a (lumped) mathematical model of HIV dynamics (utilizing g).
The presented model can be extended to incorporate drugspecific escape pathways [62,63] or realistic time-varying drug pharmacokinetics [41]. If in vivo pharmacokinetic data is available (in terms of time-varying concentrations C(t) in Eq. (1)), then extrapolation from in vitro to in vivo is possible and the mechanistic understanding of drug effects, its parametrization and extrapolation is facilitated. For RTIs and PIs, we found a nonlinear relationship between e and g (see Eqs. (S24) and (S31), Supplementary Text S1). Utilization of Eqs. (S24) and (S31) allows to simulate drug effects based on their mechanistic understanding in a lumped model, that can be compared with clinical data.
The model can also be extended to include latently infected cells (very long lived infected cells). We did not consider them in this study, since they are expected to contribute little to the dynamics analyzed herein (the first and the second decay phase).
The reproductive capacity is a useful concept to analyze and monitor drug efficacy in silico. In its current form, the reproductive capacity requires detailed knowledge about (i) the composition of the viral population, and (ii) the fitness of the different viral strains under a given treatment (reproductive numbers, Eqs. (12) and (16)-(19)).
The fitness of certain viral strains can be assessed in vitro, e.g., by phenotypic assays. We model strain specific fitness i under treatment j, in terms of two parameters: the selective disadvantage s(i), which denotes the loss in replication of mutant i, relative to the wild type; and the efficacy of treatment j against mutant i in terms of the parameter g(i,j). The selective disadvantage can, e.g., be estimated by performing a phenotypic assay with a certain mutant virus i in the absence of drug and then comparing it to the assay with the wild type. The parameter g(i,j) is already being assessed in practice (e.g., [15]), usually in terms of a fold increase in IC 50 .
Acquisition of detailed knowledge about the composition of the viral population might, due to recent advances in sequencing technology [64][65][66][67], become feasible in the future. However, novel sequencing technology requires large amounts of viral RNA, which cannot be derived when the viral load is below the limits of detection.

Realization of hybrid simulations
The overall virus dynamics in our model comprise different viral strains with copy numbers that can vary over several orders of magnitude. For this reason we have chosen a hybrid (stochastic deterministic) setting for numerical simulation. This approach (i) takes stochastic fluctuations in the slow reaction processes into account; and (ii) reduces the computational costs for the simulation of the fast (deterministic) system dynamics. We used the direct hybrid method proposed in [51]. Elementary reactions were treated stochastically whenever their propensity function or the quantity of at least one of their reactants was below a certain threshold (for all numerical simulations this threshold was set to 5). For the numerical integration of the deterministic part of the system, we implemented a solver in C++ that is based on numerical differentiation formulas [68] and uses strategies for error control and step size control comparable to ode15s in Matlab [69]. To generate the data for Table 2, we performed 1000 hybrid simulations for each condition. With realization start (t = 0), the effects of the drug treatment were simulated until the viral population size reached 90% of its pre-treatment value, i.e., virological rebound occurred. During a simulation, the stochastic partitioning of the reaction system was dynamically updated and stochastic reaction events were realized accordingly. Every numerical calculation was computed with a relative error tolerance of 10 {6 and an absolute error tolerance of 10 {9 . The hybrid simulations for Table 2 were performed on two Intel Quad-Core Xeon E5345 processors with 2.33 GHz and 32 GB RAM, which took nearly 46 hours in total or approximately 6 seconds per simulation, respectively.

Supporting Information
Text S1 This file contains the derivation of the simplified model (Fig. 3) from the detailed model (Fig. 1). Found at: doi:10.1371/journal.pcbi.1000720.s001 (0.30 MB PDF) Figure S1 Delay in the onset of viral load decay, exemplified for PI treatment. Simulation results (red line) using the novel two stage virus dynamics model and simulating 100% effective PI treatment are shown together with median clinical data (black diamonds) from PI (RTV) monotherapy. Found at: doi:10.1371/journal.pcbi.1000720.s002 (0.91 MB EPS)