The Evolution of Catalytic Function in the HIV-1 Protease

J. Mol. Biol. (2011) 408, 792–805

doi:10.1016/j.jmb.2011.02.031 Contents lists available at www.sciencedirect.com

Journal of Molecular Biology j o u r n a l h o m e p a g e : h t t p : / / e e s . e l s e v i e r. c o m . j m b

The Evolution of Catalytic Function in the HIV-1 Protease Manoj Kumar Singh 1 , Kristina Streu 2 , Andrew J. McCrone 3 and Brian N. Dominy 1 ⁎ 1

Department of Chemistry, Clemson University, Clemson, SC 29634, USA Department of Chemistry, Illinois Wesleyan University, Bloomington, Il 61701, USA 3 Department of Chemistry, University of Strathclyde, Glasgow G1 1XQ, Scotland, UK 2

Received 17 September 2010; received in revised form 19 January 2011; accepted 14 February 2011 Available online 2 March 2011 Edited by M. Levitt Keywords: HIV-1 protease; activity; enzymatic evolution; MM-PBSA; molecular dynamics

The evolution of species is a complex phenomenon based on the optimization of a multidimensional function referred to as fitness. At the level of biomolecular evolution, the fitness function can be reduced to include physiochemical properties relevant to the biological function of a particular molecule. In this work, questions involving the physical-chemical mechanisms underlying the evolution of HIV-1 protease are addressed through molecular simulation and subsequent analysis of thermodynamic properties related to the activity of the enzyme. Specifically, the impact of 40 single amino acid mutations on the binding affinity toward the matrix/ capsid (MA/CA) substrate and corresponding transition state intermediate has been characterized using a molecular mechanics Poisson–Boltzmann surface area approach. We demonstrate that this approach is capable of extracting statistically significant information relevant to experimentally determined catalytic activity. Further, no correlation was observed between the effect of mutations on substrate and transition state binding, suggesting independent evolutionary pathways toward optimizing substrate specificity and catalytic activity. In addition, a detailed analysis of calculated binding affinity data suggests that ground-state destabilization (reduced binding affinity for the substrate) could be a contributing factor in the evolutionary optimization of HIV-1 protease. A numerical model is developed to demonstrate that ground-state destabilization is a valid mechanism for activity optimization given the high concentrations of substrate experienced by the functional enzyme in vivo. © 2011 Elsevier Ltd. All rights reserved.

Introduction The evolution of organisms can be described as occurring in the context of a fitness landscape that maps genetic information onto phenotypic properties and their interaction with the environment. One class of phenotypic property is the catalytic efficien-

*Corresponding author. E-mail address: [email protected]. Abbreviation used: MA/CA, matrix/capsid.

cy of enzymes expressed within the species. Enzymes act to enhance the rate of biochemical reactions and consequently are also crucial in kinetically controlling the concentrations of metabolic products within the organism. The catalytic activity of enzymes is therefore inextricably linked to the survival of the organism and the propagation of the species. The functional properties of enzymes evolve along with the organism in order to optimize the specific metabolic functions of these organisms, thereby aiding in the ability of the organism to adapt to its environment. The essential contribution of enzymes to the fitness of an organism links the

0022-2836/$ - see front matter © 2011 Elsevier Ltd. All rights reserved.

Evolution of HIV-1 Protease Catalytic Function

evolution of these molecules to the broader evolution of species. In considering the HIV virus, the constituent enzyme HIV-1 protease is known to play an essential role in the lifecycle of the pathogen. The enzyme catalyzes the hydrolysis of specific peptide bonds releasing constituent proteins from the Gag and Pol polyproteins during viral maturation to an infectious form.1 Inhibition of this enzyme is a key treatment strategy for halting the viral lifecycle and preventing the progression of disease.2 Treatment of AIDS using protease inhibitors has further demonstrated that the evolution of HIV-1 protease is closely linked to the survival of the HIV virus. While the essential role of the protease makes it an attractive target for drug therapy, the rapid replication rate (108 –10 9 virions day – 1 ) of the virus, complemented by the high error rate of reverse transcriptase (about 1 in 10,000 base pairs),3 results in the rapid accumulation and selection of drugresistant variants of the HIV-1 protease. This drug resistance process, which is an evolutionary event, results in the selection of viral enzymes that avoid strong interactions with drugs, but still maintain sufficient catalytic efficiency for the survival of the virus. Drug resistance continues to be seen as one of the most challenging problems resulting from molecular evolution. Proper enzymatic function relies on a number of physical characteristics including fold stability and folding rate, oligomer stability and rate of oligomer formation, as well as factors coupled directly to the catalytic event, such as the affinity between the enzyme and the transition state. These properties can be considered constraints that determine the viability of specific evolutionary pathways and outcomes of an enzyme. Any mutation that has a significantly adverse affect on these properties is less likely to be naturally selected. In order to develop predictive models of evolutionary pathways and outcomes, including those related to drug resistance, the factors responsible for these evolutionary constraints should be correctly identified and characterized. In the present study, we focus on the characterization of factors governing the catalytic events associated with HIV-1 protease. Specifically, we have examined the effect of changes in substrate binding upon mutation on the catalytic activity of HIV-1 protease. It is well established that the change in transition state binding upon mutations significantly affect the catalytic activity of an enzyme.4 The direct role of substrate binding toward catalytic activity in vivo remains an important question. The characterization of factors governing the effect of mutations on the catalytic activity of HIV-1 protease can be considered as a step toward understanding the molecular evolutionary process. From the perspective of molecular modeling, this process involves modeling the reaction catalyzed by

793 the enzyme and the effect of a mutation on the catalytic process. The modeling of the reaction itself suggests the need for quantum mechanical models that properly evaluate the energy changes occurring during the bond breaking and forming process. However, modeling large molecules over long (nanoseconds or longer) time scales typically requires the use of classical potentials. Considering these two competing factors, QM/MM models have been historically applied to the study of enzyme catalysis, successfully addressing the chemical mechanisms associated with enzyme catalysis.5 However, the extensive computational time required for these QM/MM calculations constitutes a significant challenge for studying a large number of mutations, which is often necessary while characterizing broader evolutionary aspects of enzyme catalysis. The use of classical molecular mechanics models to assess enzyme catalysis can fill an important role in studying the evolution of enzyme function. Such models tend to be computationally tractable for studying the effect of a large number of mutations on enzyme activity. One of the classical approaches for studying enzyme catalysis involves the consideration of a ground state model of the transition state, often described as a transition state analogue.6 These transition state analogue models can be used in the context of a classical molecular mechanics protocol to evaluate the effect of mutations on transition state binding. The transition state binding effects can then be used to describe the effects on enzyme activity through transition state theory.7 This reduces the problem of estimating changes in activity to estimating the changes in binding affinity upon mutation, a task suited to classical molecular mechanics approaches. By using computationally tractable classical models, it is possible to characterize larger portions of the evolutionary landscape with the aim of identifying common physical mechanisms underlying the process of enzyme evolution. The calculated activity obtained for different mutants in the current study were compared to activity data determined using a bacteriophage lambda-based genetic screen system that closely captures the in vivo conditions associated with HIV1 protease function.8–10 A statistically significant correlation was observed between the calculated and experimental ranks of enzyme activity for HIV1 protease. In addition to the comparison with experiment, analysis of 40 single amino acid mutants of HIV-1 protease indicate that most mutations, as expected, destabilize transition state binding resulting in reduced activity. In addition, however, most mutations appear to enhance binding to the ground state, also resulting in reduced activity as characterized through νmax. The reduction of enzyme-catalyzed reaction rates through

Evolution of HIV-1 Protease Catalytic Function

794 ground state stabilization is consistent with the Michaelis–Menten mechanism under conditions where the substrate concentration is not significantly smaller than KM. A numerical analysis of HIV-1 protease catalysis in the immature viral particle environment suggests that the protease is acting and

evolving under such conditions. Further analysis demonstrated that the impact of mutations on ground state binding is correlated only weakly to the impact on transition state binding. This decoupling of the effect of mutations might have implications for the evolvability of enzyme catalysis.

Theoretical Background of Protease Activity Calculation The complete equation describing the peptide hydrolysis reaction catalyzed by the HIV1 protease can be written as: 1 Folded Monomer ðFÞ W Unfolded Monomer ðU Þ 2o Folded dimer k1
 Active Conformor of Enzyme ðEÞ W Enzyme − Substrate complexðESÞ k−1 þ SubstrateðSÞ 3

ð1Þ

kcat

Enzyme − Product complex ðEPÞ Y Enzyme + Product Y 5 4 The above reaction can be represented in terms of the free energy change with respect to the reaction coordinate as shown in Fig. 1. Assuming that the activation energy barrier of the first two steps is sufficiently low and the concentration of the active conformer of the enzyme [E]0, which can form a complex with the substrate, is governed primarily by the free energy difference between the folded, unfolded and dimer form of the enzyme, one can focus on steps 3, 4 and 5 of reaction (1) for the study of enzyme catalysis. For these catalytic steps, the Michaelis–Menten formulation provides a widely accepted quantitative description of enzyme catalysis.11,12 The rate of the reaction (ν), according to the Michaelis–Menten formulation, can be given as: kcat ½E0 ½S KM + ½S kcat + k − 1 = k1

v= KM

where kcat, k– 1 and k1 are the rate constants of the corresponding reaction steps.

Fig. 1. Energy profile of a peptide hydrolysis reaction catalyzed by the HIV-1 protease.

ð2Þ

Evolution of HIV-1 Protease Catalytic Function

795

A mutation in the HIV-1 protease potentially can have an effect on transition state binding, substrate binding and/or the stability of the active conformer of the protein. In the first case, where a mutation affects transition state binding, the activation energy barrier (ΔG‡cat) and consequently kcat will be affected, which will affect the rate of the corresponding reaction. In the second case, where a mutation affects substrate binding, the effect on catalytic activity is highly dependent on the concentration of the substrate. If the concentration of the substrate is low compared to the Michaelis constant (KM), the concentration of the enzyme substrate complex will be governed by the reversible step 3 of reaction (1), and the change in substrate binding affinity will affect the forward and reverse reaction rate constants (k– 1 and kcat) by the same multiplicative factor:   −DDGS exp RT Therefore, according to Eq. (2), this will have no affect on the reaction rate. However, in the case of a high concentration of substrate relative to KM, the enzyme binding sites will be saturated by the substrate and therefore the reversible step 3 of reaction (1), describing substrate binding to the enzyme, becomes less significant. In this case, the effective rate equation can be written as: v = kcat ½E0

ð3Þ

Therefore, in the case of a high concentration of substrate, a change in substrate binding affinity upon mutation will directly affect ΔG‡cat and consequently the rate of the reaction. According to Eqs. (2) and (3), the change in the concentration of enzyme caused by a change in folding free energy or dimer stability will also result in a change in the catalytic activity. In the present work, the binding affinity is calculated using a single trajectory MM-PBSA method.13–15 The use of this technique facilitates the analysis of a large number of mutations while focusing primarily on the affect of mutations on non-bonded interactions between the enzyme and the substrate/transition state analogue. The effects of mutations on monomer and dimer stability are not captured directly by the current method. The calculated binding affinity toward the transition state analogue (and substrate) does not provide any direct information about the activation energy of the catalyzed reaction; rather, it provides information about the reduction in the activation energy relative to the uncatalyzed reaction. For a general reaction, the activation energy of an uncatalyzed reaction can be related to the free energy of the substrate state (GS) and transition state (GT) as: z GTuncat − GS uncat = DG uncat

ð4Þ

and for the catalyzed reaction: GTcat − GScat = DG

z

cat

ð5Þ

The difference in the activation energies of the catalyzed and unanalyzed reaction is equivalent to the difference between the substrate and transition state binding free energies:

 DGT − DGS = GTcat − GTuncat − GS cat − GS uncat ð6Þ z z = DG cat − DG uncat Therefore, the binding affinity toward the substrate and transition state analogue can be directly related to the catalytic activity of the enzyme. The same approach can be extended to the study of the effect of a mutation on the catalytic activity of the enzyme, assuming the mutation does not significantly affect the fold stability of the active enzyme. For an enzymatic reaction taking place at rate vmax, and, therefore, at the effective rate constant kcat, the effect of the mutation on the catalytic activity can be described as:     z z vðwtÞ = vðmutÞ = exp − DG catðwtÞ − DG catðmutÞ = RT       ð7Þ = exp − DGTðwtÞ − DGS ðwtÞ − DGTðmutÞ − DGS ðmutÞ = RT However, for a reaction taking place at a reduced concentration of substrate, as discussed earlier, assuming the mutation does not affect the stability of the folded protein, the affect of a mutation on substrate binding affinity is neglected and the reaction rate is affected primarily by changes to the transition state binding

796

Evolution of HIV-1 Protease Catalytic Function

affinity. Under these conditions, it could be appropriate to calculate the effect of a mutation on the reaction rate using the following expression:     mðwtÞ =mðmutÞ ¼ exp − DGTðwtÞ −DGTðmutÞ =RT ð8Þ As a result, the effects of mutations on enzyme activity can be ranked by changes in the ground state/ transition state energy gap (at higher concentrations of substrate) or changes in the transition state binding affinity (at lower concentrations of substrate). In this study, non-parametric correlations are determined between these energy changes and the percentage change in enzyme activity.

Results and Discussion Rank correlation to the experimental data The rank correlation between the calculated activity and the experimental activity were determined for 12 systems (G17E, M36V, N37S, G40E, K43R, M46V, M46T, G48E, L63P, L63V, N88D and the wild type structure), available from a random mutagenesis activity analysis of the HIV-1 protease.8 The calculated change in activation energy based on the change in substrate or ground state (GS) binding and transition state (TS) binding (using Eq. 7) show a strong and statistically significant Spearman correlation (r = 0.81, r2 = 0.66, p = 0.001) with the experimental activity (Fig. 2). These results suggest that the affinity between the enzyme and both the transition state and the ground state substrate are factors for activity manipulation during the course of evolution in HIV-1 protease, which is consistent with the numerical model developed considering the physiological environment experienced by HIV-1 protease (see below). The MM-PBSA energies for the wild type and each mutant are given in the Supplementary Data.

The correlation between rank orders of experimental and calculated activities is not reduced significantly (r = 0.74, r2 = 0.55, p = 0.005, Fig. 3) if only the change in transition state binding free energy is considered for the activity estimation (using Eq. (8)). Mathematically, this is because the variance in TS binding due to the mutation (29.05 kcal mol– 1) is an order of magnitude larger than the variance in the GS binding (2.37 kcal mol– 1) for the mutants used in the experimental comparison. Physically, the larger variance in TS binding upon mutation is consistent with its larger average binding free energy. Stronger binding affinity towards the transition state relative to the substrate is expected because of the strong electrostatic interaction between transition state stabilizing residues in the enzyme and the transition state.16 The change in transition state binding affinity appears to be a major cause of the change in activity resulting from mutation. However, as mentioned above, the concentrations of substrate experienced in vivo by HIV-1 protease along with the correlation described in Fig. 2 suggest that substrate binding is also likely to contribute toward the catalytic activity of the enzyme. In the present study, we have chosen a nonparametric Spearman correlation17 of activities over

Fig. 2. The rank correlation between experimental and predicted rank of activity considering substrate binding along with transition state binding for activity calculation.

Evolution of HIV-1 Protease Catalytic Function

797

Fig. 3. The rank correlation between experimental and predicted rank of activity considering only transition state binding for activity calculation.

a Pearson correlation17 of activation energies for several different reasons. First, the experiment that very closely represents the environment of the HIV1 virus in vivo does not provide any direct information about the activation energies. The experiment gives information about the percentage change in activity upon mutation relative to the wild type activity.8–10 No functional relationship has been established between the experimental measure of activity (the counting of plaques as a growth score) and the activation energy. Second, our focus in the present study is to characterize evolutionary constraints and not to predict evolutionary outcomes. Therefore, the qualitative relationship between the predicted and experimentally determined activities, which validates the numerical model for HIV-1 protease catalytic activity, is sufficient within the scope of the present study. In addition to these reasons, the ability of the MM-PBSA method to estimate absolute binding affinity quantitatively has been questioned.18 However, the MM-PBSA binding affinities were found to be in excellent qualitative agreement with experiment in several earlier studies.19,20 All of these reasons suggested that a non-parametric rank correlation was most appropriate for the current study. Further analysis of the binding free energy components can provide more insight into the role of electrostatic interactions. On comparing the components of the free energy of binding for the substrate and transition state, the relatively strong affinity for the transition state compared to the substrate arises primarily through electrostatic interactions (Supplementary Data Table S1). Therefore, the electrostatic affinity toward the transition state distinguishes its binding mechanism from that

of the substrate. Looking further into the components of the binding free energy, the change in calculated activity considering both the substrate and the transition state binding (Eq. (7)) is mostly described by the change in the electrostatic component of transition state binding affinity (Table 1). The significant contribution of electrostatic interactions toward the stabilization of transition states in enzymes has been demonstrated numerous times.21–23 In addition, these results are consistent with those of a recent study where the electrostatic interactions with the transition state were shown to

Table 1. MM-PBSA free energy decomposition describing the effects of 11 single amino acid mutants on substrate and transition state binding affinity as well as the change in activation energy (ΔΔGactive = ΔGt – ΔGs)

Free energy term Total free energy Total enthalpy (ΔH) van der Waals (ΔGvdW) Coulomb (ΔGcoul) Polar solvation (ΔGpol-solv) Electrostatic (ΔGelec) Chain entropy (TΔS)

Average change in substrate binding (ΔGs)

Average change in transition state binding (ΔGt)

Average change in ΔΔGactive

–6.27(1.57)

5.49(5.39)

11.77(5.44)

–1.91(1.92)

4.02(6.81)

5.93(8.01)

–0.06(2.35)

0.24(2.87)

0.31(3.03)

–2.95(4.00)

3.60(7.36)

6.55(9.32)

0.95(3.77)

0.21(2.69)

–0.74(4.29)

–1.99(2.12)

3.81(6.32)

5.81(6.16)

4.35(2.27)

– 1.47(2.39)

-5.83(3.58)

Numbers in parentheses are standard deviations.

798 be key factors in predicting drug-resistant mutations in HIV-1 protease.24 In the present study, a simple and computationally tractable classical model approach has been shown to capture a significant portion of the information relevant to changes in enzyme activity. Using the endpoint free energy calculation method MM-PBSA, it is possible to qualitatively predict the relative activity of 11 mutants and the wild type HIV-1 protease. Earlier work demonstrated the utility of classical methods for assessing the effects of mutations on the function of enzymes.7,25 For instance, it was shown recently that classical techniques for evaluating TS binding are capable of predicting the emergence of specific mutations in HIV-1 protease in response to drug treatment.24 The present study takes this approach to characterizing enzyme activity using classical methods a step further by characterizing the relative activities of a statistically significant number of mutants. The availability of experimental data suggests HIV-1 protease as a favorable model system for the study of enzyme evolution but the low stability of the enzyme presents challenges. HIV-1 protease has a stability of ~ 1.05 kcal mol– 1 monomer,26 suggesting that a greater fraction of random single amino acid mutations could significantly disrupt the structure. This is consistent with the experimental observation that two-thirds of random mutations in HIV-1 protease result in N90% loss of activity.8 The typical limitations of fully atomic molecular mechanics studies to nanosecond timescale simulations limit the ability to detect changes in binding affinity that result, in part, from significant conformational changes in mutant enzymes. These limitations were exemplified by analyzing mutations in sequence positions 50–57, which are known to disrupt the structure of HIV-1 protease. The Spearman correlation between the activities determined from experiment and those determined from calculations decreases significantly (r2 = 0.19, p = 0.16 using Eq. (7) and r2 = 0.21, p = 0.13 using Eq. (8)) when considering mutations in positions 50–57 (I50V, F53L, I54V, R57S and R57G). Mutations at positions 50, 53 and 54 have been shown to affect the stability of the active form of HIV-1 protease.27,28 In addition, the side chain of R57 has been shown to form several important contacts within the protease,29–31 consistent with this residue being important for fold stability. Substitutions at this position with a smaller side chain, which cannot compensate for the lost hydrogen bonds and salt bridges, can have an adverse affect on the stability of the protease. This might explain, in part, the existence of R57K as the only significant polymorphism of position R57 in HIV-1 protease.32 This evidence suggests that some mutations act to disrupt the stability of the well characterized folded conformation of HIV-1 protease leading to reduced

Evolution of HIV-1 Protease Catalytic Function

activity. These instances present challenges to traditional computational techniques for characterizing the effect of random mutations on binding free energy and catalytic activity. This part of the analysis suggests that the fold stability could act as a significant constraint on the evolution of HIV-1 protease. Statistical analysis of effect of 40 mutants on enzyme activity Some approximations are clearly inherent to any modeling approach but statistical analysis of a broader dataset can provide information that is more robust than an individual measurement. In order to more broadly assess the physical mechanisms underlying the effect of mutations on enzyme activity, the study was expanded to include 29 additional mutations that were found in a population of untreated patients based on an analysis of the HIV-1 protease sequence database from the Stanford University HIV drug resistance database.32 Therefore, in examining these mutations, naturally occurring polymorphisms are studied that are also less likely to result in a significant destabilization of the native fold and the subsequent loss of all catalytic

Table 2. Impact of 40 HIV-1 protease mutations on substrate and transition state binding

Free energy term Total (ΔG) van der Waals (ΔGvdW) Coulomb (ΔGcoul) Polar solvation (ΔGpol-solv) Electrostatic (ΔGelec) Chain Entropy (TΔS)

Total (ΔG) van der Waals (ΔGvdW) Coulomb (ΔGcoul) Polar Solvation (ΔGpol-solv) Electrostatic (ΔGelec) Chain entropy (TΔS)

Mutations stabilizing substrate binding (%)

Mutations destabilizing transition state binding (%)

100 42.5

80 60

80

45

50

75

90

75

95

87.5

Average change in substrate binding (ΔGs)

Average change in transition state binding (ΔGt))

–6.79(2.69) –0.01(3.10)

4.37(5.17) 0.89(2.83)

–3.06(3.49)

0.09(5.59)

0.20(4.00)

1.74(2.59)

–2.85(1.98)

1.84(4.41)

4.09(2.11)

Numbers in parentheses are standard deviations.

–1.7(1.86)

Evolution of HIV-1 Protease Catalytic Function

799

Fig. 4. (a) The change in binding affinity toward the substrate for different mutants relative to the wild type. (b) The change in binding affinity toward the transition state for different mutants relative to the wild type.

activity. In other words, these mutations are likely to be both relevant to the natural evolution of the enzyme and amenable to the analysis used here. More insight into the binding mechanisms of the substrate and transition state analogue can be gained through analysis of the free energy components. As discussed earlier, the binding enthalpy to the transition state is more favorable than binding to the ground state primarily due to the electrostatic component (Table 2). This is followed by some compensation in –TdS, because the stronger MMPBSA energy is compensated by a weaker –TdS contribution to binding affinity (Table 2). The mutations, in general, tend to weaken binding affinity to the transition state (Fig. 4a), which appears to happen through a combination of energy

terms. The most significant contributor appears to be the electrostatic component of the binding free energy (75% of mutants show weakened affinity in this term with an average change of 1.84 kcal mol– 1 for all 40 mutants) and the binding entropy contribution, –TdS (87.5% of mutants show weakened affinity with an average change of 1.7 kcal mol– 1) as shown in Table 2. The van der Waals term appears to play some role, with 57% of mutants showing weaker binding through this term resulting in an average change of 0.8 kcal mol– 1. Mutations also result in an increase in binding affinity toward the substrate (Fig. 4b). The increase in binding affinity for the matrix/capsid (MA/CA) upon random mutation suggests that HIV-1 protease has evolved to bind less tightly to this substrate.

800 This observation is further supported by the fact that enzymes, in general, have not evolved to bind the ground state (substrate) very tightly.16,33,34 Further, this could be a consequence of the fact that HIV-1 protease is a promiscuous enzyme that is known to catalyze the hydrolysis of multiple sequences within the Gag-Pol polyprotein.35 Weak binding toward any one substrate might be a consequence of frustration in evolving promiscuous activity. This weaker binding toward the substrates might also impact the catalytic activity of the protease because this enzyme is likely acting in a high substrate concentration environment where weakened binding to the substrate (ground state destabilization) can theoretically result in a lower activation barrier and an increased catalytic rate.4 An increase in the electrostatic component of the binding free energy is the major contributor (an average of –2.85 kcal mol– 1 for all 40 mutants) toward the increase in substrate binding affinity for most (90%) of the mutations. The entropic term also plays a significant role in discriminating the effect of mutations on substrate and transition state binding; 95% of mutations show improved binding affinity toward the substrate with an average of –4.09 kcal mol– 1 through the entropic term (Table 2) for all 40 mutants. As discussed earlier, these factors affecting ground state stabilization, along with transition state destabilization, could have direct implications for protease evolution by affecting the reaction rate. Relation between change in substrate binding and transition state binding upon mutation The correlation between the calculated change in substrate and transition state binding for the 11 mutants that have been compared to experiment, as well as the larger set of 40 mutants, was determined. If the effects of mutation on substrate binding and transition state binding were highly correlated with some significant slope, mutations would have an attenuated effect on the previously mentioned energy gap and consequently a weaker effect on kcat and νmax. The lack of correlation (Supplementary Data Fig. S2) indicates the evolution of the affinity to the ground state substrate and transition states are independent, possibly resulting in faster evolution of kcat or νmax. The weak correlation between a change in substrate binding and transition state binding upon mutation could have implications in evolving activity and specificity independently. When considering moderate or low concentrations of substrate, the catalytic rate is determined by kcat/KM. Under these conditions, the effective rate constant is insensitive tod changes in substrate binding affinity while remaining sensitive toward changes in transition state binding affinity. Alternatively, the specificity of the enzyme towards particular substrates can be determined primarily through sub-

Evolution of HIV-1 Protease Catalytic Function

strate affinity.36 Therefore, the ability to evolve ground state binding affinity and affinity towards the transition state independently through different mutations could lead to the efficient and independent evolution of enzyme activity and enzyme specificity. Numerical model for the HIV-1 protease catalytic process The affinity calculations described here suggest that mutations in the wild type HIV-1 protease generally reduce the affinity toward the transition state or alternatively enhance the affinity toward the substrate. Regarding this observation from the evolutionary perspective of mutations that lead to the wild type, the calculations suggest that in the last stages of the evolution of HIV-1 protease, mutations might strengthen binding to the transition state or reduce the affinity toward the substrate. As described earlier, the strengthened binding toward the transition state can lead to enhanced activity by directly lowering the activation energy barrier. However, based on the Michaelis–Menten mechanism, the effect of weakened substrate affinity on enzyme activity is dependent on the local concentrations of the substrate. In order to correctly rationalize the catalytic mechanism of HIV-1 protease under physiological conditions, a numerical solution to the Michaelis–Menten equation is derived using experimentally determined parameters associated with the activity of HIV-1 protease. The numerical model suggests that the concentration of the substrate in the immature virus particle (before activation of the protease) is significantly higher than KM and supports the theory that ground state destabilization is a viable evolutionary mechanism for manipulating the activity of this enzyme. Reaction (1) can be rewritten as: HIVPRopen WHIVPRclosed WHIVPRclosed : Substrate WHIVPRclosed : TSWHIVPRcopen + Product

ð9Þ

In this model, a mechanism was proposed that incorporates the characteristic flap motion of HIV-1 protease. The flaps within the HIV-1 protease are known to adopt different conformations depending on whether the enzyme is bound or unbound. In the current kinetic model, the opening and closing of the flaps in the apo state is assumed to establish an equilibrium before binding the substrate. The remaining mechanism follows the traditional Michaelis–Menten model where a steady-state concentration is assumed for the enzyme/transition state complex. The equilibrium constant of the flap (Kflap) was estimated to be 7.9 × 10- 6, considering the free energy change of 7 kcal mol– 1 upon flap closing

Evolution of HIV-1 Protease Catalytic Function

based on a potential of mean force calculation.37 Assuming the transition state barrier for substrate binding and unbinding to be 10 kcal mol– 1 and 5 kcal mol– 1, respectively, the rate constant for the substrate binding (k1) and unbinding (k– 1) was calculated to be 11,303.96 Ms – 1 and 2.57 s – 1 , respectively. The calculated rate constants k1 and k2 are found to be consistent with those estimated from NMR experiments, which suggests that k1 should be b20,000 Ms– 1 and k– 1 should be a maximum of 10 s– 1.38 Therefore, in the current analysis, the simplified model of the substrate binding mechanism provides a reasonably good estimate of the rate constants for substrate binding. The calculated dissociation constant (Kd) using the present model was calculated to be 0.22 mM, which is in excellent agreement with another experimental observation that suggests Kd should be ~ 0.2 mM.38 Considering the activation energy of the reaction to be 18 kcal mol– 1,39 kcat is estimated to be 0.48 s– 1, also in agreement with experiment.40 Therefore, KM is 0.27 mM, which is consistent with the suggested range of 0.01–0.5 mM.40 Considering these calculated parameters, Fig. 5 demonstrates the relationship between the substrate concentration, the change in substrate binding affinity and the turnover rate for a standard enzyme concentration of 1 M. The concentration of the enzyme will vary as the virus matures; 41 however, the value of the enzyme concentration will affect the rate of cleavage, but will not affect the relationship toward the ground state (substrate) binding affinity. Figure 5 indicates that at low concentrations of substrate, as expected, changes in the ground state substrate binding affinity do not change the rate of the reaction. However, within a 100 nm diameter

801 immature viral particle42 comprised of 1500 gag and gag-pol precursor molecules43 (each of which contains at least 12 substrate cleavage sites) the concentration of only the MA/CA site described here is estimated to be ~ 4.7 mM (N KM). Figure 5 shows that any change in the substrate affinity at this concentration of substrate can have a significant effect on the catalytic rate. This simplistic model suggests that ground state destabilization is a plausible physical mechanism contributing to the evolution of HIV-1 protease activity within the environment of the immature HIV-1 viral particle.

Conclusions We have attempted to further characterize the role of substrate binding during the evolutionary process of the important disease target HIV-1 protease. The significant correlation between the predicted activity and the experimental activity along with the numerical model developed to describe the activity of HIV-1 protease suggests that substrate binding could be a determining factor in manipulating enzyme activity during the evolution of HIV-1 protease. Our analysis suggests that the ground state destabilization in the wild type HIV-1 protease can lead to enhanced catalytic activity in a high substrate concentration environment within the immature viral capsid. Further, we observed that changes in substrate binding and transition state binding due to mutation are not correlated. As discussed earlier, this could have implications in rapidly optimizing νmax and in more effectively evolving both enzyme specificity and catalytic activity. Finally, we have demonstrated the potential of classical force fields in predicting enzyme activity by estimating the relative activities of a statistically significant number of HIV-1 protease mutants.

Materials and Methods Modeling initial structures and molecular dynamics simulation

Fig. 5. Relationship between substrate concentration, change in substrate binding affinity and catalytic rate of the enzyme (on the z-axis).

The amino acid sequence of HIV-1 wild type protease was used as described.32 The initial coordinates of HIV-1 protease consistent with the wild type sequence were taken from the crystal structure of an inhibitor bound complex at 1.80 Å resolution (PDB ID 1HXW.pdb).44 The structures of the MA/CA (Val-Ser-Gln-Asn-Tyr/Pro-IleVal-Gln) substrate and tetra-coordinated water molecule were derived from another crystal structure of a catalytic aspartate mutated (D25N) HIV-1 protease (PDB ID 1KJ4. pdb) bound to MA/CA.45 The substrate and tetracoordinated water molecule were placed into the binding site of 1HXW by superimposing the backbone of its active site (residues with atoms within 4.5 Å of the ligand) on the

802 corresponding atoms of 1KJ4 using CHARMM.46 Another water molecule was placed at the catalytic center based on interatomic distances at the catalytic site determined earlier (Supplementary Data Fig. S1).47 The model of the MA/CA transition state analogue was created in compliance with the previously described tetrahedral anionic transition state.48 This ground state model of the transition state structure has been used successfully in the computational studies of HIV-1 protease.24 The catalytic site of the transition state analogue structure was built using the information from the crystal structure of the tetrahedral intermediate bound to HIV-1 protease49 and the distances between heavy atoms were derived from the crystal structure (Supplementary Data Fig. S1). The protonation state of Asp was used as shown in the Supplementary Data Fig. S1. Amino acid point mutations were made to the wild type protease structure using the protocol developed by Feyfant et al.50 This algorithm uses MODELLER51 to predict the positions of the atoms in the new amino acid side chain using a knowledge-based scoring function and molecular dynamics simulation methods. The mutations were made in the original inhibitor-bound crystal structure (1HXW.pdb) before replacing the inhibitor with substrate in the binding pocket. The atomic partial charges within the transition state analogue model were adopted from similar atoms in the Charmm22 parameter52 and from an earlier study.48 The bond, angle and dihedral parameters for the transition state analogue were also assigned from similar chemical systems present in the Charmm22 parameter set. All energy minimizations and molecular dynamics simulations were done with the molecular mechanics package CHARMM,46 using the 22nd version of the protein force field.52 A dielectric constant of 1 was used and the Lennard–Jones (L–J) interaction was truncated beyond 10 Å using a switch function. The Particle-Mesh Ewald (PME) method53 was used to treat the electrostatic interaction with a 10 Å cutoff. A B-spline order of 4 and real space Gaussian width of 0.34 Å– 1 was used in the PME calculation. The system was solvated in a truncated octahedron box of pre-equilibrated TIP3 water molecules54 and the distance between any atom of the protein to the nearest face of the box was maintained to at least 10 Å. The system was neutralized and the salt concentration was maintained close to 150 mM by adding KCl. Water molecules with oxygen atoms closer than 2.6 Å to any atom of the solute were deleted. The bonds involving hydrogen atoms were constrained using SHAKE in CHARMM.55 All molecular dynamics simulations were performed in an NPT ensemble. The temperature of the system was kept at 300 K on average using the Nose–Hoover thermostat. Constant pressure was maintained at 1 atm (101,325 Pa) by the Langevin piston method.56 Migration of the solute protein outside the primary solvent box was discouraged during the molecular dynamics simulation by a weak (0.5 kcal mol– 1) center-of-mass translational restraint using the MMFP module of CHARMM.57 The systems (wild type and mutant HIV1 protease structures bound with substrate or the transition state analogue) were first minimized for 500 steps using the steepest descent (SD) algorithm, keeping solute atoms fixed in order to initially optimize the solute interactions with water. In the next step, a distance-based harmonic

Evolution of HIV-1 Protease Catalytic Function restraint was applied between atoms at the catalytic site of the systems in order to restrain the distances to values given in Supplementary Data Fig. S1, with a high force constant (1000 kcal mol– 1 Å– 2). The system was then minimized with the distance-based harmonic restraints for 5000 steps using the Adopted basis Newton–Raphson (ABNR) algorithm keeping only both chains of the protein fixed, allowing the substrate/transition state analogue and water molecule(s) to optimize with respect to the binding pocket under the influence of the restraint. In the following steps, the distance-based restraints were removed from the transition state analogue-bound structures but kept for the substrate-bound structures. The systems were then minimized for 500 steps using the SD algorithm with harmonically restrained solute using a force constant of 50 kcal mol– 1 Å– 2 followed by 100 steps of SD minimization with a force constant of 10 kcal mol– 1 Å– 2 on the solute atoms. The force constant associated with distance-based restraint in the substrate-bound structure was reduced to 200 kcal mol– 1 Å– 2 before starting molecular dynamics, whereas no distance-based restraint was present for transition state analogue-bound structure at this stage. The systems were gradually heated from 100 K to 300 K during 100 ps of simulation, with a 1 fs time-step and a center of mass translation restraint. All atom-based restraints were then removed from the substrate-bound structures, leaving only center-of-mass restraints on the systems (both on substrate-bound and transition statebound structures). The molecular dynamics simulation in an NPT ensemble at a temperature of 300 K was used for equilibration (2–7 ns) and production (2 ns), using a 2 fs time step. Most of the structures were equilibrated within the initial 2 ns equilibration simulation after heating, however, for a few mutant systems the equilibration time was extended to 7 ns. MM-PBSA and entropy calculations Binding affinities were calculated using the MM-PBSA method,13–15 which is a relatively fast technique used to estimate binding affinity that has been shown to be in excellent qualitative agreement with experiment.19,20 The computational efficiency of the MM-PBSA technique permits the characterization of the effect of a large number of amino acid mutations on the activity of HIV-1 protease. According to the MM-PBSA protocol, the binding affinity of a ligand to its receptor can be estimated by calculating the average change in gas phase energy ΔUgas, the average change in the solvation free energy ΔGsolv and the average change in the conformational entropy contribution to the binding –TΔS.13–15,18 In the present study, the free energy changes were estimated from single trajectories of protein– ligand complexes (single trajectory MM-PBSA).13,15,18 The flap motion in HIV1 protease is a relatively long time scale motion and therefore unbound protease cannot be sampled accurately with relatively short (nanosecond time scale) simulations in explicit solvent.18,58 This is the primary reason for choosing a single simulation MM-PBSA method over the frequently used three simulation MM-PBSA calculation in which separate simulation trajectories of protein, ligand and complex are used for free energy estimation. The single trajectory MM-PBSA calculation method has been used successfully with the HIV1 protease in earlier studies.58,59

Evolution of HIV-1 Protease Catalytic Function The gas phase energy (Ugas) in the single trajectory MMPBSA calculation includes Coulombic energy (Gcoul) and van der Waals energy (GvdW). The gas phase energies are calculated from molecular mechanical energy function using Charmm22 parameters.52 An internal dielectric constant of 4 and no cutoff was used while evaluating the MM energies: DGbind = Gcomplex − Greceptor − Gligand gas solv GMM − PBSA = hU i + hG i − ThSi The polar part of the solvation energies (Gsolv-pol) was calculated using the Poisson–Boltzmann (PB) equation. The PBEQ module of Charmm was used to calculate the electrostatic solvation free energy. The dielectric interface defined from the molecular (contact and reentrant) surface and a grid spacing of 0.4 Ǻ was used. An internal dielectric constant of 4 and solvent dielectric constant of 80 was used in the PB calculation. The non-polar part of solvation (GSA) energy was calculated using the solvent-accessible surface area (SASA) method, where SASA was calculated using a probe radius of 1.4 Å: DGnp = gSASA + b The values of 0.00542 kcal mol– 1 for surface tension coefficient γ and 0.92 kcal mol– 1 for offset b were used as suggested.60 The change in entropy estimation was done with the VIBRAN module of CHARMM with the snapshots extracted from the production run every 40 ps.61 Normal mode analysis of the protein, ligand and protein–ligand complexes was used order to estimate the vibrational contribution to the change in entropy upon ligand binding. A distance-dependent dielectric with definition of 4r was used for the normal mode calculations, where r is the pairwise distance between atoms. The electrostatic interaction was truncated beyond 12 Å using the shift function. The loss of transitional and rotational degrees of freedom upon protein–ligand binding was estimated from the mass and moment of inertia of the protein,62 ligand and protein–ligand complex evaluated within the VIBRAN module of CHARMM. Rank correlation of calculated and experimental activities In calculating the Spearman correlation, the lowest rank was assigned to the mutant with the lowest experimental activity and higher rank was assigned to mutants with greater experimental activity. Similarly, the lowest rank was assigned to the mutant with the lowest calculated activity and higher rank was assigned to mutants with larger calculated activity. The Spearman correlation (ρ) between the experimental rank (xi) and the calculated rank (yi) was calculated as: P i ðxi − xÞðyi − yÞ ffi U = qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi P 2 2 i ðxi −xÞ ðyi −yÞ The p-values associated with the correlation coefficients were determined using the two-tailed test based on the t-distribution.17

803

Acknowledgements This work was supported by a National Science Foundation Career award (MCB-0953783).

Supplementary Data Supplementary data associated with this article can be found, in the online version, at doi:10.1016/j. jmb.2011.02.031

References 1. Brik, A. & Wong, C. H. (2003). HIV-1 protease: mechanism and drug discovery. Org. Biomol. Chem. 1, 5–14. 2. Wlodawer, A. & Vondrasek, J. (1998). Inhibitors of HIV-1 protease: a major success of structure-assisted drug design. Annu. Rev. Biophys. Biomol. Struct. 27, 249–284. 3. Coffin, J. M. (1995). HIV population dynamics in vivo: implications for genetic variation, pathogenesis, and therapy. Science, 267, 483–489. 4. Garcia-Viloca, M., Gao, J., Karplus, M. & Truhlar, D. G. (2004). How enzymes work: analysis by modern rate theory and computer simulations. Science, 303, 186–195. 5. Acevedo, O. & Jorgensen, W. L. (2009). Advances in quantum and molecular mechanical (QM/MM) simulations for organic and enzymatic reactions. Acc. Chem. Res. 43, 142–151. 6. Eksterowicz, J. E. & Houk, K. N. (1993). Transitionstate modeling with empirical force fields. Chem. Rev. 93, 2439–2461. 7. Pan, Y., Gao, D., Yang, W., Cho, H. & Zhan, C. (2007). Free energy perturbation (FEP) simulation on the transition states of cocaine hydrolysis catalyzed by human butyrylcholinesterase and its mutants. J. Am. Chem. Soc. 129, 13537–13543. 8. Parera, M., Fernandez, G., Clotet, B. & Martinez, M. A. (2007). HIV-1 protease catalytic efficiency effects caused by random single amino acid substitutions. Mol. Biol. Evol. 24, 382–387. 9. Sices, H. J. & Kristie, T. M. (1998). A genetic screen for the isolation and characterization of site-specific proteases. Proc. Natl Acad. Sci. USA, 95, 2828–2833. 10. Martinez, M. A., Cabana, M., Parera, M., Gutierrez, A., Este, J. A. & Clotet, B. (2000). A bacteriophage lambda-based genetic screen for characterization of the activity and phenotype of the human immunodeficiency virus type 1 protease. Antimicrob. Agents Chemother. 44, 1132–1139. 11. Briggs, G. E. & Haldane, J. B. S. (1925). A note on the kinetics of enzyme action. Biochem. J. 19, 339. 12. Michaelis, L. & Menten, M. (1913). Die Kinetik der Invertinwirkung. Biochem. Z. 49, 333–369. 13. Gilson, M., Given, J., Bush, B. & McCammon, J. (1997). The statistical-thermodynamic basis for computation of binding affinities: a critical review. Biophys. J. 72, 1047–1069.

Evolution of HIV-1 Protease Catalytic Function

804 14. Kollman, P. A., Massova, I., Reyes, C., Kuhn, B., Huo, S., Chong, L. et al. (2000). Calculating structures and free energies of complex molecules: combining molecular mechanics and continuum models. Acc. Chem. Res. 33, 889–897. 15. Srinivasan, J., Cheatham, T. E., Cieplak, P., Kollman, P. A. & Case, D. A. (1998). Continuum solvent studies of the stability of DNA, RNA, and Phosphoramidateâ DNA helices. J. Am. Chem. Soc. 120, 9401–9409. 16. Villa, J. & Warshel, A. (2001). Energetics and dynamics of enzymatic reactions. J. Phys. Chem. B, 105, 7887–7907. 17. Devore, J. L. (1991). Probability and Statistics for Engineering and the Sciences, 3rd edit. Duxbury Press, Belmont, CA. 18. Lee, M. S. & Olson, M. A. (2006). Calculation of absolute protein-ligand binding affinity using path and endpoint approaches. Biophys. J. 90, 864–877. 19. Rastelli, G., Del Rio, A., Degliesposti, G. & Sgobba, M. (2010). Fast and accurate predictions of binding free energies using MM-PBSA and MM-GBSA. J. Comput. Chem. 31, 797–810. 20. Thorsteinsdottir, H. B., Schwede, T., Zoete, V. & Meuwly, M. (2006). How inaccuracies in protein structure models affect estimates of protein-ligand interactions: computational analysis of HIV-I protease inhibitor binding. Proteins, 65, 407–423. 21. Bruice, T. C. (2006). Computational approaches: reaction trajectories, structures, and atomic motions. enzyme reactions and proficiency. Chem. Rev. 106, 3119–3139. 22. Gao, J., Ma, S., Major, D. T., Nam, K., Pu, J. Z. & Truhlar, D. G. (2006). Mechanisms and free energies of enzymatic reactions. Chem. Rev. 106, 3188–3209. 23. Warshel, A., Sharma, P. K., Kato, M., Xiang, Y., Liu, H. & Olsson, M. H. M. (2006). Electrostatic basis for enzyme catalysis. Chem. Rev. 106, 3210–3235. 24. Ishikita, H. & Warshel, A. (2008). Predicting drugresistant mutations of HIV protease13. Angew Chem. Int. Ed. 47, 697–700. 25. Yang, W., Pan, Y., Zheng, F., Cho, H., Tai, H. & Zhan, C. (2009). Free-energy perturbation simulation on transition states and redesign of butyrylcholinesterase. Biophys. J. 96, 1931–1938. 26. Noel, A. F., Bilsel, O., Kundu, A., Wu, Y., Zitzewitz, J. A. & Matthews, C. R. (2009). The folding freeenergy surface of HIV-1 protease: insights into the thermodynamic basis for resistance to inhibitors. J. Mol. Biol. 387, 1002–1016. 27. Liu, F., Boross, P. I., Wang, Y., Tozser, J., Louis, J. M., Harrison, R. W. & Weber, I. T. (2005). Kinetic, stability, and structural changes in high-resolution crystal structures of HIV-1 protease with drug-resistant mutations L24I, I50V, and G73S. J. Mol. Biol. 354, 789–800. 28. Liu, F., Kovalevsky, A. Y., Tie, Y., Ghosh, A. K., Harrison, R. W. & Weber, I. T. (2008). Effect of flap mutations on structure of HIV-1 protease and inhibition by saquinavir and darunavir. J. Mol. Biol. 381, 102–115. 29. Freedberg, D. I., Ishima, R., Jacob, J., Wang, Y., Kustanovich, I., Louis, J. M. & Torchia, D. A. (2002). Rapid structural fluctuations of the free HIV protease flaps in solution: relationship to crystal structures and

30.

31.

32.

33. 34.

35.

36.

37.

38.

39.

40.

41.

42. 43.

44.

comparison with predictions of dynamics calculations. Protein Sci. 11, 221–232. Martin, P., Vickrey, J. F., Proteasa, G., Jimenez, Y. L., Wawrzak, Z., Winters, M. A. et al. (2005). “Wideopen” 1.3 Å structure of a multidrug-resistant HIV-1 protease as a drug target. Structure, 13, 1887–1895. Meiselbach, H., Horn, A., Harrer, T. & Sticht, H. (2007). Insights into amprenavir resistance in E35D HIV-1 protease mutation from molecular dynamics and binding free-energy calculations. J. Mol. Mod. 13, 297–304. Rhee, S., Gonzales, M. J., Kantor, R., Betts, B. J., Ravela, J. & Shafer, R. W. (2003). Human immunodeficiency virus reverse transcriptase and protease sequence database. Nucleic Acids Res. 31, 298–303. Dinner, A. R., Blackburn, G. M. & Karplus, M. (2001). Uracil-DNA glycosylase acts by substrate autocatalysis. Nature, 413, 752–755. Fromme, J. C., Bruner, S. D., Yang, W., Karplus, M. & Verdine, G. L. (2003). Product-assisted catalysis in base-excision DNA repair. Nat. Struct. Biol. 10, 204–211. Kohl, N. E., Emini, E. A., Schleif, W. A., Davis, L. J., Heimbach, J. C., Dixon, R. A. et al. (1988). Active human immunodeficiency virus protease is required for viral infectivity. Proc. Natl Acad. Sci. USA, 85, 4686–4690. Lopes, A., Busch, M. S. A. & Simonson, T. (2010). Computational design of protein-ligand binding: modifying the specificity of asparaginyl-tRNA synthetase. J. Comput. Chem. 31, 1273–1286. Rick, S. W., Erickson, J. W. & Burt, S. K. (1998). Reaction path and free energy calculations of the transition between alternate conformations of HIV-1 protease. Proteins: Struct. Funct. Bioinf. 32, 7–16. Katoh, E., Louis, J. M., Yamazaki, T., Gronenborn, A. M., Torchia, D. A. & Ishima, R. (2003). A solution NMR study of the binding kinetics and the internal dynamics of an HIV-1 protease-substrate complex. Protein Sci. 12, 1376–1385. Piana, S., Bucher, D., Carloni, P. & Rothlisberger, U. (2004). Reaction mechanism of HIV-1 protease by hybrid Car-Parrinello/classical MD simulations. J. Phys. Chem. B, 108, 11139–11149. Tözsér, J., Bláha, I., Copeland, T. D., Wondrak, E. M. & Oroszlan, S. (1991). Comparison of the HIV-1 and HIV-2 proteinases using oligopeptide substrates representing cleavage sites in Gag and Gag-Pol polyproteins. FEBS Lett. 281, 77–80. Darke, P. L., Jordan, S. P., Hall, D. L., Zugay, J. A., Shafer, J. A. & Kuo, L. C. (1994). Dissociation and association of the HIV-1 protease dimer subunits: equilibria and rates. Biochemistry, 33, 98–105. ÖHagen, Å., Luftig, R. B., Reicin, A. S., Yin, L. E. I., Ikuta, K., Kimura, T. et al. (1997). The Morphology of the immature HIV-1 virion. Virology, 228, 112–114. Layne, S. P., Merges, M. J., Dembo, M., Spouge, J. L., Conley, S. R., Moore, J. P. et al. (1992). Factors underlying spontaneous inactivation and susceptibility to neutralization of human immunodeficiency virus. Virology, 189, 695–714. Kempf, D. J., Marsh, K. C., Denissen, J. F., McDonald, E., Vasavanonda, S., Flentge, C. A. et al. (1995). ABT538 is a potent inhibitor of human immunodeficiency

Evolution of HIV-1 Protease Catalytic Function

45.

46.

47.

48.

49.

50. 51. 52.

virus protease and has high oral bioavailability in humans. Proc. Natl Acad. Sci. USA, 92, 2484–2488. Prabu-Jeyabalan, M., Nalivaika, E. & Schiffer, C. A. (2002). Substrate shape determines specificity of recognition for HIV-1 protease: analysis of crystal structures of six substrate complexes. Structure, 10, 369–381. Brooks, B. R., Brooks, C. L., 3rd, Mackerell, A. D., Jr, Nilsson, L., Petrella, R. J., Roux, B. et al. (2009). CHARMM: the biomolecular simulation program. J. Comput. Chem. 30, 1545–1614. Okimoto, N., Tsukui, T., Hata, M., Hoshino, T. & Tsuda, M. (1999). Hydrolysis mechanism of the phenylalanine−proline peptide bond specific to HIV1 protease: investigation by the ab initio molecular orbital method. J. Am. Chem. Soc. 121, 7349–7354. Bjelic, S. & Aqvist, J. (2004). Computational prediction of structure, substrate binding mode, mechanism, and rate for a malaria protease with a novel type of active site. Biochemistry, 43, 14521–14528. Kovalevsky, A. Y., Chumanevich, A. A., Liu, F., Louis, J. M. & Weber, I. T. (2007). Caught in the act: the 1.5 Å resolution crystal structures of the HIV-1 protease and the I54V mutant reveal a tetrahedral reaction intermediate. Biochemistry, 46, 14854–14864. Feyfant, E., Sali, A. & Fiser, A. (2007). Modeling mutations in protein structures. Protein Sci. 16, 2030–2041. Sali, A. & Blundell, T. L. (1993). Comparative protein modelling by satisfaction of spatial restraints. J. Mol. Biol. 234, 779–815. MacKerell, A. D., Bashford, D., Bellott, M., Dunbrack, R. L., Evanseck, J. D., Field, M. J. et al. (1998). All-atom empirical potential for molecular modeling and dynamics studies of proteins. J. Phys. Chem. B, 102, 3586–3616.

805 53. Darden, T., York, D. & Pedersen, L. (1993). Particle mesh Ewald: an N.log(N) method for Ewald sums in large systems. J. Chem. Phys. 98, 10089–10092. 54. Jorgensen, W. L., Chandrasekhar, J., Madura, J. D., Impey, R. W. & Klein, M. L. (1983). Comparison of simple potential functions for simulating liquid water. J. Chem. Phys. 79, 926–935. 55. van Gunsteren, W. F. & Berendsen, H. J. C. (1977). Algorithms for macromolecular dynamics and constraint dynamics. Mol. Phys. 34, 1311–1327. 56. Feller, S. E., Zhang, Y., Pastor, R. W. & Brooks, B. R. (1995). Constant pressure molecular dynamics simulation: the Langevin piston method. J. Chem. Phys. 103, 4613–4621. 57. Beglov, D. & Roux, B. (1994). Finite representation of an infinite bulk system: solvent boundary potential for computer simulations. J. Chem. Phys. 100, 9050–9063. 58. Hou, T. & Yu, R. (2007). Molecular dynamics and free energy studies on the wild-type and double mutant HIV-1 protease complexed with amprenavir and two amprenavir-related inhibitors: mechanism for binding and drug resistance. J. Med. Chem. 50, 1177–1188. 59. Stoica, I., Sadiq, S. K. & Coveney, P. V. (2008). Rapid and accurate prediction of binding free energies for saquinavir-bound HIV-1 proteases. J. Am. Chem. Soc. 130, 2639–2648. 60. Sitkoff, D., Sharp, K. A. & Honig, B. (1994). Correlating solvation free energies and surface tensions of hydrocarbon solutes. Biophys. Chem. 51, 397–409. 61. Brooks, B. R., Jane, D. & Karplus, M. (1995). Harmonic analysis of large systems. I. Methodology. J. Comput. Chem. 16, 1522–1542. 62. McQuarrie, D. A. (1973). Statistical Thermodynamics. University Science Books, Mill Valley, CA.