Monte Carlo calculations on HIV-1 reverse transcriptase complexed with the non-nucleoside inhibitor 8-Cl TIBO: contribution of the L100I and Y181C variants to protein stability and biological activity

Marilyn B.Kroeger Smith1,2, Michelle L. Lamb3, Julian Tirado-Rives4, William L. Jorgensen4, Christopher J. Michejda1, Sandra K. Ruby5 and Richard H. Smith, Jr1,5

1 National Cancer Institute-FCRDC, P.O. Box B, Frederick, MD 21702, 3 Department of Pharmaceutical Chemistry, University of California, San Francisco, CA 94143, 4 Department of Chemistry, Yale University, New Haven, CT 06520 and 5 Western Maryland College, Westminster, MD 21157, USA


    Abstract
 Top
 Abstract
 Introduction
 Methods
 Results and discussion
 References
 
A computational model of the non-nucleoside inhibitor 8-Cl TIBO complexed with HIV-1 reverse transcriptase (RT) was constructed in order to determine the binding free energies. Using Monte Carlo simulations, both free energy perturbation and linear response calculations were carried out for the transformation of wild-type RT to two key mutants, Y181C and L100I. The newer linear response method estimates binding free energies based on changes in electrostatic and van der Waals energies and solvent-accessible surface areas. In addition, the change in stability of the protein between the folded and unfolded states was estimated for each of these mutations, which are known to emerge upon treatment with the inhibitor. Results from the calculations revealed that there is a large hydrophobic contribution to protein stability in the native, folded state. The calculated absolute free energies of binding from both the linear response, and also the more rigorous free energy perturbation method, gave excellent agreement with the experimental differences in activity. The success of the relatively rapid linear response method in predicting experimental activites holds promise for estimating the activity of the inhibitors not only against the wild-type RT, but also against key protein variants whose emergence undermines the efficacy of the drugs.

Keywords: HIV-1/inhibitor/reverse transcriptase/mutation/modeling


    Introduction
 Top
 Abstract
 Introduction
 Methods
 Results and discussion
 References
 
A major obstacle in the successful treatment of AIDS patients with HIV-1 reverse transcriptase (RT) inhibitors is the emergence of newly resistant variant strains. The takeover of these mutant strains causes a rapid reduction in the potency of both nucleoside and non-nucleoside inhibitors during the initial course of treatment, although the specific mutations responsible for the decline in efficacy differ between the two classes of drugs (Boyer et al, 1994aGo,bGo). One way to combat this problem is to combine inhibitors in treatment protocols that would induce different resistance profiles, with the ultimate goal of leaving only escape mutants that replicate poorly.

Several key drug-induced variants have been identified following treatment with the non-nucleoside inhibitor 8-chlorotetrahydroimidazo[4,5,1-jk][1,4]benzodiazepin-2(1H)-thione (8-Cl TIBO) (Moermans et al., 1995Go) including Y181C, Y188H, L100I, K103N and V106A. These residues directly surround the hydrophobic binding pocket of the complexed protein (Figure 1Go). From structural studies of the unliganded form of the enzyme (Esnouf et al., 1995Go; Rodgers et al., 1995Go; Hsiou et al, 1996Go), some of these residues also appear to be critical in maintaining the stability of the folded protein.



View larger version (71K):
[in this window]
[in a new window]
 
Fig. 1. Close-up of the binding of 8-Cl TIBO (blue) in the non-nucleoside pocket of RT. Amino acid residues shown in purple are readily mutated upon treatment with the inhibitor.

 
In particular, L100 makes important van der Waals contacts with both Y181 and Y188 (Figure 2Go). Prior to inhibitor binding, the side chains of W229, Y181 and Y188 are proposed to undergo major structural rearrangement in preparation for holding the bound inhibitor in the binding pocket in an energetically favorable orientation (Rodgers et al., 1995Go). Following complexation, the position of L100 is seen to be at the junction of the two regions of the `butterfly-shaped' pocket, contacting both `wings' of the inhibitor (Ding et al., 1995Go) (Figure 3Go). This interaction leads to a substantial increase in stabilization of the inhibitor from that seen in its unbound state (Kroeger Smith et al., 1995Go). Clearly, replacement of leucine with isoleucine could also affect the position of other nearby residues, leading to distortions in the binding pocket. Previous molecular modeling experiments carried out by energy minimization of the protein–inhibitor complex revealed that the L100I change resulted in a ~5 kcal/mol destabilization of the RT–8-Cl TIBO complex (Kroeger Smith et al., 1997Go). These structural changes around the binding pocket could also lead to distortions in the geometry at the nearby polymerase active site, causing a marked reduction in the polymerase activity in the unbound enzyme, as has been suggested previously (Kroeger Smith et al., 1997Go). In fact, measurement of the DNA polymerase activity of the I100 variant showed it to be 40% that of the wild-type enzyme (Kroeger Smith et al., 1997Go).



View larger version (41K):
[in this window]
[in a new window]
 
Fig. 2. Close-up of the future binding pocket area in unliganded RT. Amino acid residue 100 (shown in blue) makes contact with residues Y181 and Y188 (shown in red) which undergo major structural rearrangement prior to non-nucleoside inhibitor binding. L100 also makes contact with K101 and K103 (purple) in the uncomplexed state, residues which mutate upon with the inhibitor.

 


View larger version (12K):
[in this window]
[in a new window]
 
Fig. 3. Diagram of the non-nucleoside binding pocket showing the orientation of 8-Cl TIBO relative to key residues involved in binding.

 
Since the Y181 (and Y188) residues have been implicated in {pi}-stacking stabilization of the bound inhibitors (Ding et al., 1995Go; (Kroeger Smith et al., 1995Go, 1997Go), replacement of tyrosine by cysteine would also be expected to destabilize the drug–protein complex. One would thus expect that the loss of these interactions would result in less favorable interactions between the inhibitor and the enzyme. Indeed, measurement of this interaction energy following minimization for the Y181C substitution in the 8-Cl TIBO complex resulted in a ~1.0 kcal/mol reduction (Smith et al., 1998bGo). Another recent piece of biochemical evidence suggests that another non-nucleoside inhibitor, nevirapine, shows a 500-fold decrease in affinity for the Y181C mutant relative to the wild-type RT (Spence et al., 1996Go). However, the polymerase activity of this mutant enzyme is not reduced at all, suggesting that this side chain is spatially oriented so as not to affect the polymerase site geometry. Recent crystallographic information from the C181 mutant structure (Das et al., 1996Go) supports this hypothesis; very little adjustment in protein backbone or side chain or inhibitor position is seen in the region surrounding the polymerase active site (overall r.m.s. by backbone superposition = 0.621 Å).

As a means of assessing free energy changes during the complexation of non-nucleoside inhibitors with RT presented with the L100I and Y181C mutations, Monte Carlo free energy perturbation and linear response calculations were carried out for the inhibitor 8-Cl TIBO complexed with RT. Free energy perturbation calculations using this Monte Carlo method have previously been carried out for a series of trypsin benzamidine inhibitors (Essex et al., 1997Go), on orthogonal cyclosporin–cyclophilin pairs (Pierce and Jorgensen, 1997Go) and for neurotropic inhibitors of FK506 binding protein (Lamb and Jorgensen, 1998Go). Linear response calculations of this type have been validated in studies of a series of sulfonamide inhibitors of human thrombin (Hertoz-Jones and Jorgensen, 1997Go), FKBP12 inhibitors (Lamb et al., 1999Go) and TIBO inhibitors of HIV-1 reverse transcriptase (Smith et al., 1998aGo). Two new methods, adaptive chemical Monte Carlo molecular dynamics and Poisson–Boltzmann/solvent accessibility calculations (Ericksson et al., 1999Go) have also been used to predict the activity of 8-Cl TIBO derivatives. In addition to the above energy calculations on liganded and unliganded RT, the change in protein stability due to the L100I and Y181C mutations was estimated for each mutant separately by computing the difference in unfolding free energy between model wild-type and mutant unbound proteins. Methodology of this latter type has successfully been applied to proteins such as barnase (Prevost et al., 1991Go; Sun et al., 1996Go). The data from these simulations, and also the resultant implications for the prediction of inhibitor activity and protein stability in the presence of continuous protein mutation, are discussed.


    Methods
 Top
 Abstract
 Introduction
 Methods
 Results and discussion
 References
 
Basis of the method

The application of free energy calculations to assess complexation energy and protein stability following mutation can be illustrated in the thermodynamic cycle shown in Figure 4Go. In this figure, the folded wild-type or mutant proteins are denoted Pfb in their bound form and as Pf in their unbound form. In the unfolded extended states, the uncomplexed solvated protein is denoted Pu and the unsolvated (gas-phase) protein is Pr. The {Delta}G values shown above the horizontal arrows could, in principle, be obtained from experiment, and the {Delta}G values alongside the vertical arrows can be calculated from free energy perturbation simulations as the energy difference between the wild-type and variant proteins. The determination of these values allows one to assess the difference in relative solvation free energy and protein stability through the calculation of mutational free energies.



View larger version (14K):
[in this window]
[in a new window]
 
Fig. 4. Thermodynamic cycles used to describe {Delta}G for the mutation of leucine to isoleucine in RT.

 
The free energy difference between wild-type and mutant systems, for example, {Delta}Gfb(wt->mut), the left-most vertical in Figure 4Go, is exactly described by Equation 1Go (Straatsma, 1996Go):

where kB is the Boltzmann constant, T is the temperature of the system and {Delta}U is the energy difference between wild-type and mutant bound complexes. The < >wt notation represents an ensemble average over configurations of the wild-type complex generated during a Monte Carlo (or molecular dynamics) simulation.

The accuracy of free energy perturbation calculations performed in this way depends upon sampling of all configurations making relevant contributions to the free energy during the simulation, the accuracy of the energy function in describing the interactions within the ensemble and that the ensemble represents configurations of both wild-type and mutant states. In practice, for this final condition to hold, perturbations are performed through a number of intermediate steps (windows), with the free energies computed at each step summed to generate that of the full mutation. While Equation 1Go could in principle be used to compute absolute binding affinities, sampling issues limit its usefulness to protein and ligand mutations. Therefore, relative free binding energies are best obtained through thermodynamic cycles (Figure 4Go).

In addition to free energy perturbation calculations, linear response calculations were carried out to obtain the free energies of binding of 8-Cl TIBO with wild-type and various mutant proteins. The difference in these energies could then be compared with the binding free energies calculated above. In contrast to free energy perturbation calculations, the linear interaction energy or linear response approximation, originally proposed by Aqvist et al. (1994), makes use of ensemble averaging to estimate absolute binding free energies directly, avoiding the intermediate steps of FEP calculations. The form used here for computing the free energy of binding for a given protein and ligand, extended from the original method to better describe ligand solvation (Hertoz-Jones and Jorgensen, 1997Go) is

where the coefficients {alpha}, ß and {gamma} were determined to obtain the best fit to the experimental observable, the free energy of binding. If the surroundings exhibit a linear response to the electrostatic field of the ligand, then by definition, the value of ß will be 1/2. Deviations of linearity may be overcome by assignment of alternative values of ß for special chemical cases (Hansson et al., 1998Go) or by fitting to experimental binding data to optimize this parameter as well.

Construction of the site

Calculations for the inhibitor–RT complex were based on the crystallographic structure for the 8-Cl TIBO (Ding et al., 1995Go) bound in the wild-type RT binding pocket. From the structural coordinates, a model of the inhibitor binding pocket that included amino acid residues approximately within 20 Å of the inhibitor was constructed, as has been documented previously (Kroeger Smith et al., 1995Go). This final model was composed of 149 amino acid residues in eight separable strands, each one terminated with either an acetyl or methylamino group at their N- and C-terminal ends, respectively. Basic (Arg and Lys) and acidic (Asp and Glu) amino acid residues were assumed to be in their charged states. The 133 amino acid residues comprising the site were (flexible) from the p66 domain 94–111, 176–193, 223–237 and from p51 136–138; (fixed) from p66 90–93, 112–113, 157–175, 194–206, 214–222, 238–242, 266–270, 314–322, 347–350, 381, 383 and from p51 133–135, 139–141. Since no structure had been reported for the C181 variant at the time these calculations were begun, this mutant was simulated by replacing the tyrosine residue with cysteine using the Biopolymer module of InsightII (Molecular Simulations). Likewise, the I100 variant was simulated by replacing the leucine residue with isoleucine. Subsequent availability of the 181C mutant coordinates (Das et al., 1996Go) allowed for construction of a true variant model in a fashion analogous to that for the wild-type complex, which was also used in the simulations.

All sites were surrounded by a 6 Å layer of water and minimized with the AMBER force field using the Discover module of InsightII to relieve crystal packing strains. The three-stage minimization protocol for this calculation was the same as that reported previously (Kroeger Smith et al., 1995Go). This minimized site was then used as input for the MC simulations. In these calculations, the protein was partially surrounded by a 26 Å sphere of TIP3P water molecules that was centered at the geometric center of complex. Any water molecules that had serious steric conflicts were deleted, resulting in 1575 water molecules that remained during the simulations. A 1.5 kcal/mol Å half-harmonic restraining force was applied to waters at the surface of the sphere to prevent their migration away from the complex. The perturbations were run at 37°C with 9 Å radius spherical, residue-based cutoffs for all non-bonded interactions. The list of non-bonded interactions was updated every 2x105 configurations during the simulations. For the protein, a partial united atom approach was used in which all hydrogens are explicit except those on saturated carbon atoms. The positions of the backbone atoms and all bond lengths of the protein were fixed throughout the simulations. In a fashion similar to the bound protein, a model of the unbound protein (Rodgers et al., 1995Go) containing exactly the same residues as in the bound state was constructed, minimized and solvated.

The starting conformation for the bound inhibitor was taken from the crystal structure of its respective complex. The 8-Cl TIBO inhibitor has the R absolute configuration for N6, with the tertiary amino nitrogen of the diazepine ring bearing the dimethylallyl substituent. All atoms of the inhibitor were explicitly included in the calculations. Parameters for the inhibitor were assigned from the OPLS all-atom force field (Jorgensen and Tirado-Rives, 1998Go), except for the dihedral torsion parameters for CN6CC dihedral angle, which were determined as described previously (Smith et al., 1998aGo). Partial charges for the inhibitors were obtained by fitting to the electrostatic potential surface (EPS) from 6–31G* ab initio calculations using the ChelpG routine (Breneman and Wiberg, 1990Go; Carlson et al., 1993Go) in the Gaussian 94 program (Frisch et al., 1994Go). In 8-Cl TIBO, the thiazolobenzimidazole ring system and the six atoms associated with the allylic double bond were held internally rigid during the simulations, since variations in the geometry of these portions of the inhibitor are likely to be minimal. Since these regions were linked by geometry elements that can vary, they were still able to move relative to one another and would follow rigid body motions of the inhibitor.

A model of the denatured state of the Leu to Ile transformation was constructed akin to that made for the representative mutation in barnase (Sun et al., 1996Go). The peptide Ace–Gly–Leu–Lys–Ame was assembled, representing residues 99–101 of RT. Ace and Ame denote the acetyl and N-methyl groups used to cap the N- and C-terminal ends of the peptide, respectively. A chloride counterion was placed close to the charged Lys residue in order to neutralize the system. The peptide was constructed in an extended conformation with the {phi} and {Psi} angles of Gly99, Leu100 and Lys101 all equal to 180° and then minimized in vacuo with the AMBER force field using the Discover module of InsightII. To simulate the tyrosine to cysteine mutation, the extended structure of residues 180–182 was constructed in similar fashion and was capped as described above. The minimized structures were then solvated with a 26 Å cap of TIP3P water molecules.

Free energy perturbations

The transformation of wild-type to variant RT was carried out in the same manner for both the folded and unfolded states. For the Leu to Ile mutation, the CD2 atom of leucine was converted to a dummy atom in isoleucine, while a dummy atom in leucine was converted to CG2 in isoleucine (Figure 5Go). For each dummy atom, one force constant of one of the bond angles was changed to zero during the course of the simulation. For the Tyr to Cys mutation, the mutation was carried out in two separate stages; first tyrosine was converted to alanine and then alanine was changed to cysteine (Figure 6Go). The first stage was further divided by first removing the charges on the hydroxyl and the CZ atoms during the course of one FEP simulation. Then the CG, CD1, CD2, CE1, CE2, CZ and OH atoms in Tyr were all transformed to dummy atoms and the bond length from the CG dummy to CB was reduced to 0.65 Å while all the remaining bonds to dummy atoms were shrunk to 0.35 Å. In addition, the CB atom type was changed to the appropriate type for alanine. In the second stage of the transformation, the original CG atom of Tyr was mutated to the SG atom in Cys and a dummy atom in Tyr was converted to the HG atom in Cys; both the HG and SG stoms were uncharged. During this mutation, the bond lengths for the CG–SG and SG–HG bonds were changed to 1.81 and 1.35 Å, respectively. Lastly, these SG and HG atoms were converted to their respective counterparts with the correct charges.



View larger version (10K):
[in this window]
[in a new window]
 
Fig. 5. Atom changes in the mutation of leucine to isoleucine at position 100 in RT for both the folded and unfolded states.

 



View larger version (27K):
[in this window]
[in a new window]
 
Fig. 6. Atom changes in the mutation of tyrosine to cysteine at position 181 in RT for both (a) the folded and (b) the unfolded states.

 
All calculations were performed using the program MCPRO, Version 1.6 (Jorgensen, 1998Go) implemented on R10000 Silicon Graphics workstations. The solution-phase calculations were carried out by first subjecting the peptides, proteins or protein–inhibitor complexes, in separate simulations, to 10x106 configurations of equilibration with only water moving. Next, after 2x106 configurations of equilibration in each window had been completed with the inhibitor and the protein (constrained as described above) allowed to move, averaging was performed for 4x106 configurations. In all, 10 windows were used with double-width sampling for each simulation in which the final structure from each window was used as the starting structure for the next window. The gas-phase calculations were carried out in an analogous fashion, except for the obvious omission of the water equilibration step. In order to determine the complexation, unfolding and solvation free energies of the mutations Leu to Ile or Tyr to Cys, the free energy changes {Delta}Gfb and {Delta}Gf in the folded state and {Delta}Gu and {Delta}Gr in the unfolded state (see Figure 4Go) were calculated in separate simulations.

The relative magnitude of the interactions between the drug and individual residues in the binding pocket of the complexes following the simulations was determined using the Docking module within the InsightII package. Superposition of structures was also carried out using InsightII.

Linear response calculations

MC simulations of the unbound protein (either wild-type or 100I or 181C mutants) consisted of an initial solvent equilibration of 10x106 configurations in which the solute remained fixed, followed by an additional 10x106 configurations of equilibration in which the protein (partially constrained as described above) was also permitted to move. Energy averages were then accumulated over a final 10x106 configurations, during which snapshots of the geometry were saved after each 2x105 configurations. Separate running averages were maintained for the Lennard–Jones (Uvdw) and Coulombic (Uelec) components of the non-bonded interaction energies throughout the simulation. The standard deviation for each of these components was <=1.5% of the mean, indicating that the system had achieved a state of apparent equilibration. Simulations for the bound inhibitor–protein complex (where the ligand remained constant and the protein was mutated) again consisted of three phases: 10x106 configurations of equilibration with water only moving, 4x106 configurations of equilibration in which the inhibitor and protein (partially constrained as described above) were allowed to move and a final 8x106 configurations during which energy averages and geometry snapshots were acquired. Separate running averages were obtained for the Lennard–Jones and Coulombic components for both protein–solvent and inhibitor–protein non-bonded interactions. As in the case of the unbound inhibitor simulations, the standard deviations for each of these components was <=1.5% of the mean. These data, along with a term for the change in the solvent accessible surface area (SASA) of the protein, were then fitted into the linear response equation.

Inhibitor activity

The EC50 values for the variants simulated in this study have been reported previously (Kroeger Smith et al., 1997Go) or are unpublished (R.W.Buckheit, personal communication). These values were converted to estimated experimental free energy values using the equation {Delta}Gbinding = –RTlnKform at 37°C. Here Kform, the equilibrium constant for the formation of the complex, is taken to be the reciprocal of the EC50 value. This estimate was required because Ki data are not available for TIBO complexes with RT. Given the very close similarity in structure of these variants, however, it is reasonable to expect that EC50 and Ki, while not equivalent, should be linearly related (Cheng and Prusoff, 1973Go).


    Results and discussion
 Top
 Abstract
 Introduction
 Methods
 Results and discussion
 References
 
Free energy perturbation simulations: binding free energies

The calculated free energy curves for the mutation of RT residue 100 (Leu to Ile) are shown in Figure 7Go for the bound (with 8-Cl TIBO) and unbound folded states ({Delta}Gfb and {Delta}Gf), as well as the unbound extended states ({Delta}Gu and {Delta}Gr) (see Figure 4Go for an explanation of this notation). The overall free energy differences are given in Table IGo. The curves in Figure 7Go are quite smooth and the {Delta}G values for each window are generally <1.0 kcal/mol, with associated statistical uncertainties ({sigma}) of <0.1 kcal/mol for each step. The error values for the total {Delta}Gs are, on average, 4% of the total energy change.



View larger version (16K):
[in this window]
[in a new window]
 
Fig. 7. Computed free energy curves for the L100I mutation in RT for 8-Cl TIBO.

 

View this table:
[in this window]
[in a new window]
 
Table I. Experimental and computed free energy perturbation valuesa for the L100I and Y181C mutations in RT
 
For the Y181C mutation, the free energies for the {Delta}Gfb and {Delta}Gf, as well as the {Delta}Gu and {Delta}Gr states, are shown in Table IIGo for all four steps of the calculation, while the values for the total mutation are given in Table IGo. The curves for each of the four steps of this simulation (not shown) were also smooth, with associated statistical uncertainties of <0.1 kcal/mol for each window.


View this table:
[in this window]
[in a new window]
 
Table II. Free energy changesa in kcal/mol for the four-step perturbationb of Y to C at residue 181 in RT
 
The change in binding free energy obtained from the FEP calculations, {Delta}{Delta}Gfb to f = {Delta}Gf{Delta}Gfb, is –1.5 kcal/mol for Leu to Ile at amino acid residue 100 when 8-Cl TIBO is complexed with the enzyme. The negative value obtained for {Delta}{Delta}Gfb to f is indicative of the better binding of 8-Cl TIBO to the wild-type protein residue than the mutant isoleucine variant. This value is in excellent agreement with the experimental {Delta}{Delta}G value of –2.1 kcal/mol. For the Y181C change, the sum of the four independently obtained energies ({Delta}Gfb) was subtracted from the {Delta}Gf sum value to give a {Delta}{Delta}Gfb value of +0.32 kcal/mol. This value agrees fairly well with the experimental {Delta}{Delta}G value of –0.4 kcal/mol, keeping in mind that the experimental {Delta}{Delta}G value was determined from EC50 values, instead of the more desirable Ki data.

Given the reasonable match between thermodynamic and experimental results with the 8-Cl TIBO–RT complex for the L100I mutation, the structures from the FEP simulations were analyzed to determine the origin of the –1.5 kcal/mol final difference in energy between the wild-type and mutant proteins. Examination of the position of the inhibitor in the starting wild-type L100 and averaged variant I100 complexes showed an overall displacement of the inhibitor's position of 1.14 Å (Figure 8Go), most likely in response to the `appearance' of the methyl group in isoleucine. In its wild-type orientation, the dimethylallyl portion of the compound has been shown to interact favorably with p66 amino acid residues P95, Y181, Y188, G190 and W229, while the benzodiazepine ring system interacts with residues K101, K103, V106, F227, H235, P236 and Y318 (Ding et al., 1995Go; Das et al., 1996Go; Smith et al., 1998aGo,bGo). At the middle of the pocket, L100, V170 and L234 make contact with both portions of the molecule. During the course of the mutation, adjustment of the inhibitor's position caused shifts in both portions of the molecule. For example, the {gamma} carbon in the dimethylallyl chain (see Figure 3Go) can be seen to move away from Y181 by 0.33 Å and closer to both W2229 and Y188 by ~0.60 Å. Meanwhile, the thiocarbonyl carbon moved closer to Y318 (0.52 Å) and to L100 (1.13 Å), while moving away from K103 by 0.36 Å. This repositioning of the drug is reflected in the measured energy of interaction of TIBO with this residue; in the L100 complex it is –5.4 kcal/mol (totally from van der Waals contributions), whereas in the I100 mutant structure it has increased to –4.4 kcal/mol, indicative of a less favorable interaction upon mutation. In addition, the one observed hydrogen bond in the crystal structure of the complex (Ding et al., 1995Go), TIBO:N1H to the backbone carbonyl oxygen of Lys101 (1.71 Å), was substantially lengthened during the course of the simulation from 1.90 Å in the starting (AMBER minimized) wild-type structure to 2.84 Å in the mutant complex. It has been noted previously (Smith et al., 1998aGo) that the most active inhibitors are those that are strongly hydrogen-bonded in this position, where a strong bond is defined as <2.5 Å. Clearly, the lengthened bond in the mutant structure supports this hypothesis since TIBO has a decreased ability to inhibit I100 RT. No new hydrogen bonds were formed during the course of the simulation.



View larger version (21K):
[in this window]
[in a new window]
 
Fig. 8. Change in position of 8-Cl TIBO and key side chain amino acid residues following free energy pertubation calculations on the mutation L100I in RT.

 
Structural changes for Y181C mutation, where TIBO's activity for wild-type vs mutant protein is essentially zero, were significantly less. The position of the inhibitor was seen to adjust only in the diallyl `tail' portion of the molecule. The greatest movement was in the terminal carbon atoms (~ 0.5 Å) of the tail. This small amount of movement is reflected in the minimal difference in the interaction energy of the drug with the wild-type tyrosine and the mutant 181 cysteine residue in the protein of ~0.01 kcal/mol. In addition, the hydrogen bond length was unchanged following the free energy perturbation.

Free energy perturbation simulations: role of hydration

In order to attempt a separation of intramolecular effects and differences in hydration in the wild-type and mutant unfolded tetrapeptides, the role of hydration can first be examined. The mutational free energy in the vacuum (reference state) {Delta}Gr for the Leu to Ile change is 2.1 kcal/mol. Subtraction of the free energy change in the unfolded solvated state ({Delta}Gu = 3.9 kcal/mol) from that in the reference state results in a solvation free energy difference ({Delta}{Delta}Gr to u = {Delta}Gu{Delta}Gr) of 1.8 kcal/mol. Using 1-butane and isobutane as model compounds, Wolfenden et al. (1981) estimated a {Delta}{Delta}Ghydration value of –0.13 kcal/mol while Ooi et al. (1987) estimated a value of +0.24 kcal/mol, both of which are smaller than our calculated value. This implies that ß-branching results in hindrance of the solvation of the backbone amide groups. For the Y181C change, the calculated solvation free energy difference is larger (6.4 kcal/mol) than that for the leucine to isoleucine mutation, which is expected. An estimation of {Delta}{Delta}Ghydration for this mutation using p-cresol and methanethiol as the model compounds gave a value of 4.9 kcal/mol (Toolpenden et al., 1981Go); this again is smaller by the same amount (~1.6 kcal/mol) than the calculated value in our system. This again suggests that ß-branching plays a role in solvation.

For the Leu to Ile mutation, the total free energy difference between leucine and isoleucine for the model peptide in the gas phase vs for the folded, unbound protein ({Delta}{Delta}Gr to f = {Delta}Gf {Delta}Gr) of 12.6 kcal/mol is significantly larger than the corresponding difference in the solvation free energy {Delta}{Delta}Gr to u (1.8 kcal/mol). This suggests that hydration makes only a small contribution to the difference in stabilization free energy between the leucine and isoleucine side chains in RT. Likewise, the value for {Delta}{Delta}Gr to f for the Y to C change is 14.1 kcal/mol, again much larger than the solvation free energy value of 6.4 kcal/mol. The remaining differences are likely to be a reflection of the intramolecular destabilization of the folded protein caused by the mutations.

The values given in Table IGo for the various free energy differences of the mutations could in principle be utilized to estimate the effect of substitution on the folding stability of RT. For the L100I mutation, the difference between the unbound folded and unfolded solvated states ({Delta}{Delta}Gf to u = {Delta}Gu{Delta}Gf) is 10.8 kcal/mol, while the corresponding value for the Y181C mutation is 7.7 kcal/mol. While the simulations in both systems restricted backbone movement, the lack of motion is a worse approximation in the folded protein because its effect is more severe than in the isolated tetrapeptide (Xu et al., 1998Go). Accordingly, the calculated values are bound to be an overestimate, although the qualitative predictions are probably correct since the majority of single point mutations in proteins are destabilizing (Xu et al., 1998Go). The smaller effect of the Y181C mutation is consistent with its prevalence in resistance to NNIRT inhibition (De Clercq, 1998Go).

Since hydration, as mentioned above, does not greatly favor the stability of the leucine side chain over the isoleucine side chain in the unfolded protein, clearly interactions with water are much more important in the folded state. The results cited above most likely reflect the more `buried' position of tyrosine 181, as opposed to L100, the latter being exposed to solvent molcules surrounding the protein. In spite of the fact that we cannot determine the contribution to {Delta}G from non-covalent interactions from our calculations, the docking energy data show that the measured energy of interaction between the inhibitor and residue 100 or 181 arises solely from van der Waals contributions, which supports the contention that non-covalent interactions are the major driving force in the unfolding free energy of RT. Thus, taken together, these data suggest that the large decrease in stability observed upon mutation to Ile100 or Cys181 from native RT is indicative of the important role that these residues play in maintaining the integrity of the folded protein. The proof for this hypothesis, however, remains to be established.

Linear response simulations

In order to determine if the faster linear response calculations could possibly yield reasonable correlations between the experimental values and the calculated energies, the average interaction energies of the bound and unbound proteins obtained from wild-type and five RT variants complexed with 8-Cl TIBO are shown in Table IIIGo and the solvent-accessible surface area (SASA) values are shown in Table IVGo. The resulting difference in energy between the complexed and uncomplexed inhibitors in the electrostatic, van der Waals and SASA components (Table VGo) were then fitted into the linear response equation. The predicted values for {Delta}Gbinding for wild-type and five mutants using all three components in the linear response equation gave a surprisingly reasonable result with an r.m.s. error of 1.34 kcal/mol (Table VIGo, model 1). When the fit was calculated by eliminating any one parameter (Table VIGo, models 2–4) or setting the {alpha} constant to 0.1 (Table VIGo, model 5), the r.m.s. error did not improve. If the values for the V108I mutation were eliminated from the fit, the r.m.s. error improved to 0.82 kcal/mol if all three components were included (Table VIGo, model 1a). Not only did this reduce the r.m.s. error, it also decreased the difference between the predicted and experimental values compared with model 1. One reason for this may be that the valine to isoleucine is the only mutation where the size of the amino acid residue increased, which may introduce more error into that calculation. If V108I was eliminated from the fit and {gamma} was set to zero (Table VIGo, model 4a), the r.m.s. error increased to 1.20 kcal/mol. Elimination of the energy values for any other mutation aside from V108I did not improve the r.m.s. error no matter how many components were included in the fit (data not shown).


View this table:
[in this window]
[in a new window]
 
Table III. Average interaction energies from bound and unbound protein MC simulations (kcal/mol)a
 

View this table:
[in this window]
[in a new window]
 
Table IV. Calculated average solvent-accessible surface area (Å) for bound and unbound proteina
 

View this table:
[in this window]
[in a new window]
 
Table V. Unweighted changes in energy and solvent-accessible surface area from MC simulations using bound and unbound protein variants
 

View this table:
[in this window]
[in a new window]
 
Table VI. Comparison of experimental and calculated {Delta}Gbinding for various parameter fitting models
 
Given the many possible inherent errors in computing the difference in energy between bound and unbound proteins that have been `mutated' from the wild-type structure by changing only one residue in the complex while keeping the inhibitor the same, the above r.m.s. errors are very reasonable. First, correct assessment of interaction energy terms for residues that are buried in the protein is something that has not yet been clearly determined. Second, the fit was carried out on wild-type protein and the only five `neutral' mutations (all uncharged to uncharged residue changes) for which activity data exist, with three parameters varied in the fitting procedure. In the context of these limitations, this approach may have merit for `quick' assessment of the activity of a particular inhibitor against key mutants in RT.

Calculations were also carried out on a site made from the crystal structure of the mutant 181C protein–8-Cl TIBO complex. The corresponding energy values that were determined from the calculation that correspond to those in Table IIIGo were considerably more negative than those in the `artificially mutated' site, especially for the Uelec value for the bound protein-solv term in the linear response equation (–5274.06 vs –4836.64 kcal/mol). This is not unexpected since the water molecules would be expected to solvate the mutant crystal structure differently. Following the Monte Carlo simulation, backbone superposition of the final position of the inhibitor in the site constructed from mutation of the wild-type protein's crystal structure vs that made from the mutant protein's crystal structure gave an overall r.m.s. value of 0.80 Å (recall that backbone motion occurred only during the minimization procedure, not during the Monte Carlo simulation). The atomic positions of 8-Cl TIBO did not differ by more than 1.4 Å on average between the two structures. Taken together, these results suggest that the mutation of a wild-type crystal structure is a reasonable approach for carrying out these calculations and possibly alleviates the long wait for mutant crystal structure coordinates.

Conclusions

The calculated {Delta}{Delta}Gbinding values obtained from the predicted values from the linear response calculations (Table VIGo, model 1a) for the RT–8-Cl TIBO complex for the difference between wild-type and L100I or Y181C mutant proteins (–0.7 and +0.6 kcal/mol, respectively) agree well with the corresponding FEP value of –1.5 kcal/mol for the L100I variation and +0.3 kcal/mol for the Y181C change. In addition, the predicted results also correspond well to their respective experimental values of –2.1 and –0.4 kcal/mol. Thus, these results imply that linear response calculations of this type may be useful in predicting, at least approximately, the activity of the inhibitors against key mutations in an enzyme. While FEPs are certainly the standard computational method to evaluate such mutations, the success in applying linear response calculations to estimate the activity of inhibitors against key variants is especially promising since they require only one tenth of the calculation time of free energy perturbations.

The results of this study demonstrate that FEP calculations on two key mutants (L100I and Y181C) in the HIV-1 RT–8-Cl TIBO protein complex produce {Delta}{Delta}G values that agree well with the experimentally determined activities of these two mutations. In addition, the energies calculated from individual Monte Carlo linear response simulations on the wild-type and mutant proteins correlate well with their respective experimental {Delta}G values. The fact that these theoretical methods show promise in the prediction of these activities, in particular the more time-efficient linear response calculations, means that the activity of inhibitors against key variants can be ascertained prior to the initiation of lengthy synthetic and crystallographic procedures. Since the emergence of these resistance-engendering variants following treatment with RT non-nucleoside inhibitors is a major obstacle in the treatment of AIDS, the application of computational methodology to this problem will hopefully speed the search for a cure for this deadly disease.


    Notes
 
2 To whom correspondence should be addressed Back


    Acknowledgments
 
The authors acknowledge the Frederick Biomedical Supercomputing Center for access to the Cray-YMP and Silicon Graphs R10000 computers. This research was sponsored in part by the National Cancer Institute, DHHS, under contract with ABL.


    References
 Top
 Abstract
 Introduction
 Methods
 Results and discussion
 References
 
Aqvist,J., Medina,C. and Samuelsson,J.-E. (1994) Protein Engng, 7, 385–391.[Abstract]

Boyer,P.L, Ding,J., Arnold,E. and Hughes,S.H. (1994a) Antimicrob. Agents Chemother., 38, 1909–1914.[Abstract]

Boyer,P.L, Tantillo,C., Jacobo-Molina,A., Nanni,R.G., Ding,J., Arnold,E. and Hughes,S.H. (1994b) Proc. Natl Acad. Sci. USA, 91, 4882–4886.[Abstract]

Breneman,C.M. and Wiberg,K.B. (1990) J. Comput. Chem., 11, 361–373.[ISI]

Carlson,H.A., Nguyen,M., Orozco,M. and Jorgensen,W.L. (1993) J. Comput. Chem., 14, 1240–1249.[ISI]

Cheng,Y. and Prusoff,W.H. (1973) Biochem. Pharmacol., 22, 3099–3108.[ISI][Medline]

Das,K. et al. (1996) J. Mol. Biol., 264, 1085–1100.[ISI][Medline]

De Clercq,E. (1998) Antiviral Res., 38, 153–179.[ISI][Medline]

Ding,J., Das,K., Moereels,H., Koymans,L., Andries,K., Janssen,P.A.J., Hughes,S.H. and Arnold,E. (1995) Nature Struct. Biol., 2, 407–415.[ISI][Medline]

Ericksson,M.A.L., Pitera,J. and Kollman,P.A. (1999) J. Med. Chem., 42, 868–881.[ISI][Medline]

Esnouf,R., Ren,J., Ross,C., Jones,Y., Stammers,D. and Stuart,D. (1995) Struct. Biol., 2, 303–308.

Essex,J.W., Severance,D.L., Tirado-Rives,J. and Jorgensen,W.L. (1997) J. Phys. Chem. B, 36, 9663–9669.

Frisch,M.J. et al. (1994) Gaussian 94. Gaussian, Pittsburgh, PA.

Hansson,T., Marelius,J., Aqvist,J. (1998) Comput.-Aided Mol. Des., 12, 27–35.

Hertoz-Jones,D.K. and Jorgensen,W.L. (1997) J. Med. Chem., 40, 1539–1549.[ISI][Medline]

Hsiou,Y., Ding,J., Das,K., Clark,A.D.,Jr, Hughes,S.H. and Arnold,E. (1996) Structure, 4, 853–860.[ISI][Medline]

Jorgensen,W.L. (1998) MCPRO, Version 1.6. Yale University, New Haven, CT.

Jorgensen,W.L. and Tirado-Rives,J. (1998) J. Am. Chem. Soc., 110, 1657–1666.

Kroeger Smith,M.B. et al. (1995) Protein Sci., 4, 2203–2222.[Abstract/Free Full Text]

Kroeger Smith,M.B., Michejda,C.J, Hughes,S.H., Boyer,P.L, Janssen,P.A.J., Andries,K., Buckheit,R.W.,Jr and Smith,R.H.,Jr (1997) Protein Engng, 10, 1379–1383.[Abstract]

Lamb,M.L and Jorgensen,W.L. (1998) J. Med. Chem., 41, 3928–3939.[ISI][Medline]

Lamb,M.L, Tirado-Rives,J. and Jorgensen,W.L. (1999) Bioorg. Med. Chem., 7, 851–860.[ISI][Medline]

Moermans,M. et al. (1995) In Abstracts of the Fourth International Workshop on HIV-1 Drug Resistance, Sardinia, Italy.

Ooi,T., Oobatake,M., Nemethy, G and Scherage,H.A. (1987) Proc. Natl Acad. Sci. USA, 84, 3086–3090.[Abstract]

Pierce,A.C. and Jorgensen,W.L. (1997) Angew. Chem., Int. Ed. Engl., 36, 1466–1469.[ISI]

Prevost,M., Wodak,S.J., Tidor,B. and Karplus,M. (1991) Proc. Natl Acad. Sci. USA, 88, 10880–10884.[Abstract]

Rodgers,D.W., Gamblin,S.J., Harris,B.A., Ray,S., Culp,J.S., Hellmig,B., Woolf,D.J., Debouck,C. and Harrison,S.C. (1995) Proc. Natl. Acad. Sci. USA, 92, 1222–1226.[Abstract]

Smith,R.H.,Jr, Jorgensen,W.L., Tirado-Rives,J., Lamb,M.L., Janssen,P.A.J., Michejda,C.J. and Kroeger Smith,M.B. (1998a) J. Med. Chem., 41, 26, 5272–5286.[ISI][Medline]

Smith,R.H.,Jr, Michejda,C.J., Hughes,S.H., Arnold,E.A., Janssen,P.A.J. and Kroeger Smith,M.B. (1998b) J. Mol. Struct. (Theochem), 423, 67–77.

Spence,R.A., Anderson,K.S. and Johnson,K.A. (1996) Biochemistry, 35, 1054–1063.[ISI][Medline]

Straatsma,T.P. (1996) Rev. Comput. Chem., 9, 81–127.

Sun,Y.-C., Veenstra,D.L. and Kollman,P.A. (1996) Protein Engng, 9, 273–281.[Abstract]

Toolpenden,R., Andersson,L., Cullis,P.M., Southgate,C.C.B. (1981) Biochemistry, 20, 849–855.[ISI][Medline]

Wolfenden,R., Andersson,L., Cullis,P.M. and Southgate,C.C.B. (1981) Biochemistry, 20, 849–855.[ISI][Medline]

Xu,J., Baase,W.A., Baldwin,E. and Matthews B.W. (1998) Protein Sci., 7, 158–177.[Abstract/Free Full Text]

Received September 7, 1999; revised January 18, 2000; accepted February 29, 2000.





This Article
Abstract
FREE Full Text (PDF)
Alert me when this article is cited
Alert me if a correction is posted
Services
Email this article to a friend
Similar articles in this journal
Similar articles in ISI Web of Science
Similar articles in PubMed
Alert me to new issues of the journal
Add to My Personal Archive
Download to citation manager
Search for citing articles in:
ISI Web of Science (15)
Request Permissions
Google Scholar
Articles by Smith, M. B.K.
Articles by Smith, R. H.
PubMed
PubMed Citation
Articles by Smith, M. B.K.
Articles by Smith, R. H., Jr