Precision of free energies calculated by molecular dynamics simulations of peptides in solution

Perspective: Alchemical free energy calculations for drug discovery

Computational techniques see widespread use in pharmaceutical drug discovery, but typically prove unreliable in predicting trends in protein-ligand binding. Alchemical free energy calculations seek to change that by providing rigorous binding free energies from molecular simulations. Given adequate sampling and an accurate enough force field, these techniques yield accurate free energy estimates. Recent innovations in alchemical techniques have sparked a resurgence of interest in these calculations. Still, many obstacles stand in the way of their routine application in a drug discovery context, including the one we focus on here, sampling. Sampling of binding modes poses a particular challenge as binding modes are often separated by large energy barriers, leading to slow transitions. Binding modes are difficult to predict, and in some cases multiple binding modes may contribute to binding. In view of these hurdles, we present a framework for dealing carefully with uncertainty in binding mode or conformation in the context of free energy calculations. With careful sampling, free energy techniques show considerable promise for aiding drug discovery.

Computational ligand-based rational design: Role of conformational sampling and force fields in model development

A significant number of drug discovery efforts are based on natural products or high throughput screens from which compounds showing potential therapeutic effects are identified without knowledge of the target molecule or its 3D structure. In such cases computational ligand-based drug design (LBDD) can accelerate the drug discovery processes. LBDD is a general approach to elucidate the relationship of a compound's structure and physicochemical attributes to its biological activity. The resulting structure-activity relationship (SAR) may then act as the basis for the prediction of compounds with improved biological attributes. LBDD methods range from pharmacophore models identifying essential features of ligands responsible for their activity, quantitative structure-activity relationships (QSAR) yielding quantitative estimates of activities based on physiochemical properties, and to similarity searching, which explores compounds with similar properties as well as various combinations of the above. A number of recent LBDD approaches involve the use of multiple conformations of the ligands being studied. One of the basic components to generate multiple conformations in LBDD is molecular mechanics (MM), which apply an empirical energy function to relate conformation to energies and forces. The collection of conformations for ligands is then combined with functional data using methods ranging from regression analysis to neural networks, from which the SAR is determined. Accordingly, for effective application of LBDD for SAR determinations it is important that the compounds be accurately modelled such that the appropriate range of conformations accessible to the ligands is identified. Such accurate modelling is largely based on use of the appropriate empirical force field for the molecules being investigated and the approaches used to generate the conformations. The present chapter includes a brief overview of currently used SAR methods in LBDD followed by a more detailed presentation of issues and limitations associated with empirical energy functions and conformational sampling methods.

Absolute free energy calculations by thermodynamic integration in four spatial dimensions

An optimized technique for calculating the excess chemical potential of small molecules in dense liquids and the binding affinity of molecular ligands to biomolecules is reported. In this method, a molecular species is coupled to the system of interest via a nonphysical fourth spatial dimension w through which insertion or extraction can be carried out [R. Pomes, E. Eisenmesser, C. B. Post et al., J. Chem. Phys. 111, 3387 (1999)]. Molecular simulations are used to compute the potential of mean force (PMF) acting on the solute molecule in the fourth dimension. The excess chemical potential of that molecule is obtained as the difference in the PMF between fully coupled and fully decoupled systems. The simplicity, efficiency, and generality of the method are demonstrated for the calculation of the hydration free energies of water and methanol as well as sodium, cesium, and chloride ions. A significant advantage over other methods is that the 4D-PMF approach provides a single effective and general route for decoupling all nonbonded interactions (i.e., both Lennard-Jones and Coulombic) at once for both neutral and charged solutes. Direct calculation of the mean force from thermodynamic integration is shown to be more computationally efficient than calculating the PMF from umbrella sampling. Statistical error analysis suggests a simple strategy for optimizing sampling. The detailed analysis of systematic errors arising from the truncation of Coulombic interactions in a solvent droplet of finite size leads to straightforward corrections to ionic hydration free energies.

Unexpected relative aqueous solubilities of a phosphotyrosine analogue and two phosphonate derivatives

Phosphotyrosine (pTyr) is an essential component of biological signaling, often being a determinant of protein-protein interactions. Accordingly, a number of drug discovery efforts targeting signal transduction pathways have included phosphotyrosine and analogues as essential components of the lead compounds. Toward the goal of improved biological efficacy, the phosphonate and difluoro phosphonate analogues of pTyr have been employed in inhibitor design because of their stability to hydrolysis and enhanced binding affinity in certain cases. To quantitate the contribution of aqueous solubility of pTyr, phosphonomethyl phenylalanine (Pmp), and difluorophosphonomethyl phenylalanine (F(2)Pmp) to their relative binding affinities, free energy perturbation calculations were undertaken on the mimetics phenol phosphate (PP), benzyl phosphonate (BP), and difluorobenzyl phosphonate (F(2)BP), including development of empirical force field parameters compatible with the CHARMM all-atom force fields. Notably, it is shown that the most favorably solvated compound of the series is BP, followed by PP, with F(2)BP the least favorably solvated for both the mono- and dianionic forms of the compounds. The molecular origin of this ordering is shown to be due to changes in charge distribution, in the comparatively larger size of the fluorine atoms, as well as in differences of local solvation between PP and BP. The implications of the differences in aqueous solubility toward the relative binding potencies of pTyr-, Pmp-, and F(2)Pmp-containing peptide ligands are discussed. Our results indicate that one general principle explaining the efficacy of selective fluorination to enhance binding affinities may lie in the ability of fluorine atoms to increase the hydrophobicity of a ligand while maintaining its capability to form hydrogen bonds.

Modeling protein-small molecule interactions: structure and thermodynamics of noble gases binding in a cavity in mutant phage T4 lysozyme L99A

The complexes of phage T4 lysozyme L99A with noble gases have been studied by molecular dynamics simulation. In a long simulation of the complex with one Xe atom, the structure was found to undergo global conformation change involving a reversible opening and closing of the entrance to the substrate-binding site, during which the conformations of the N and C-terminal domains varied little. The distributions of Xe positions sampled in dynamics simulations were refined in terms of anisotropic Gaussian distributions via least-squares minimization of the difference between Fourier transforms. In addition, molecular transformation simulations have been applied in order to calculate the binding free energies of Xe, Kr and Ar relative to a standard state at a pressure of 1 bar. A single bound Xe is found to assume an equilibrium distribution over three adjacent preferred sites, while in a two-Xe complex, the two Xe atoms preferentially occupy two of these. The positions of the three sites agree closely with the positions of bound Xe determined in the refined crystal structure of a complex formed at a pressure of 8 bar Xe, and the calculated affinities agree well with the observed partial occupancies. At a pressure of 8 bar, a mixture of one-Xe and two-Xe complexes is present, and similarly for complexes with Kr and Ar, with single occupancy relatively more prevalent with Kr and Ar. (Binding of a third Xe atom is found to be quite unfavorable.) A comparison with simulation results for the binding of benzene to the same site leads to the conclusion that binding of Xe within cavities in proteins is common because of several favorable factors: (1) Xe has a large atomic polarizability; (2) Xe can be applied at a relatively high pressure, i.e. high chemical potential; (3) an unfavorable entropic term related to the need to orient the ligand in the binding site is absent. Finally, it is found that the model's binding energy of a water molecule in the cavity is insufficient to overcome the unfavorable binding entropy.

Non-Boltzmann thermodynamic integration (NBTI) for macromolecular systems: relative free energy of binding of trypsin to benzamidine and benzylamine

The relative free energies of binding of trypsin to two amine inhibitors, benzamidine (BZD) and benzylamine (BZA), were calculated using non-Boltzmann thermodynamic integration (NBTI). Comparison of the simulations with the crystal structures of both complexes, trypsin-BZD and trypsin-BZA, shows that NBTI simulations better sample conformational space relative to thermodynamic integration (TI) simulations. The relative binding free energy calculated using NBTI was much closer to the experimentally determined value than that obtained using TI. The error in the TI simulation was found to be primarily due to incorrect sampling of BZA's conformation in the binding pocket. In contrast, NBTI produces a smooth mutation from BZD to BZA using a surrogate potential, resulting in a much closer agreement between the inhibitors' conformations and the omit electron density maps. This superior agreement between experiment and simulation, of both relative binding free energy differences and conformational sampling, demonstrates NBTI's usefulness for free energy calculations in macromolecular simulations.

Protein-protein recognition: an experimental and computational study of the R89K mutation in Raf and its effect on Ras binding

Binding of the protein Raf to the active form of Ras promotes activation of the MAP kinase signaling pathway, triggering cell growth and differentiation. Raf/Arg89 in the center of the binding interface plays an important role determining Ras-Raf binding affinity. We have investigated experimentally and computationally the Raf-R89K mutation, which abolishes signaling in vivo. The binding to [gamma-35S]GTP-Ras of a fusion protein between the Raf-binding domain (RBD) of Raf and GST was reduced at least 175-fold by the mutation, corresponding to a standard binding free energy decrease of at least 3.0 kcal/mol. To compute this free energy and obtain insights into the microscopic interactions favoring binding, we performed alchemical simulations of the RBD, both complexed to Ras and free in solution, in which residue 89 is gradually mutated from Arg into Lys. The simulations give a standard binding free energy decrease of 2.9+/-1.9 kcal/mol, in agreement with experiment. The use of numerous runs with three different force fields allows insights into the sources of uncertainty in the free energy and its components. The binding decreases partly because of a 7 kcal/mol higher cost to desolvate Lys upon binding, compared to Arg, due to better solvent interactions with the more concentrated Lys charge in the unbound state. This effect is expected to be general, contributing to the lower propensity of Lys to participate in protein-protein interfaces. Large contributions to the free energy change also arise from electrostatic interactions with groups up to 8 A away, namely residues 37-41 in the conserved effector domain of Ras (including 4 kcal/mol from Ser39 which loses a bifurcated hydrogen bond to Arg89), the conserved Lys84 and Lys87 of Raf, and 2-3 specific water molecules. This analysis will provide insights into the large experimental database of Ras-Raf mutations.

Conformation of the Ras-binding domain of Raf studied by molecular dynamics and free energy simulations

Recognition of Ras by its downstream target Raf is mediated by a Ras-recognition region in the Ras-binding domain (RBD) of Raf. Residues 78-89 in this region occupy two different conformations in the ensemble of NMR solution structures of the RBD: a fully alpha-helical one, and one where 87-90 form a type IV beta-turn. Molecular dynamics simulations of the RBD in solution were performed to explore the stability of these and other possible conformations of both the wild-type RBD and the R89K mutant, which does not bind Ras. The simulations sample a fully helical conformation for residues 78-89 similar to the NMR helical structures, a conformation where 85-89 form a 3(10)-helical turn, and a conformation where 87-90 form a type I beta-turn, whose free energies are all within 0.3 kcal/mol of each other. NOE patterns and H(alpha) chemical shifts from the simulations are in reasonable agreement with experiment. The NMR turn structure is calculated to be 3 kcal/mol higher than the three above conformations. In a simulation with the same implicit solvent model used in the NMR structure generation, the turn conformation relaxes into the fully helical conformation, illustrating possible structural artifacts introduced by the implicit solvent model. With the Raf R89K mutant, simulations sample a fully helical and a turn conformation, the turn being 0.9 kcal/mol more stable. Thus, the mutation affects the population of RBD conformations, and this is expected to affect Ras binding. For example, if the fully helical conformation of residues 78-89 is required for binding, its free energy increase in R89K will increase the binding free energy by about 0.6 kcal/mol.

Energetic decomposition of the alpha-helix-coil equilibrium of a dynamic model system

Using molecular dynamics simulations to calculate free energies of molecular transformation, we have computed helix-coil transition free energies for alanine oligomers up to 14 residues long. The simulations have been done on the model in vacuo with dielectric constant, epsilon = 1, 5, 25, and infinity and on the model in solution with explicit representation of water molecules and with partial charges on the oligomer set to zero. (The analogous simulations of the solvated model with full charges on the oligomer were reported elsewhere [L. Wang et al. (1995) Proceedings of the National Academy of Science USA 92, 10924-10928]). In vacuo, both entropic and electrostatic contributions oppose formation of a 3-residue helical nucleus in the helix initiation step. The entropy change opposing helix growth is found to be 3 e.u., van der Waals interactions favor helix growth by 1.9 kcal/mol, and electrostatic interactions favor helix growth by 3 kcal/mol (for epsilon = 1; all these values are per residue). In water, helix stability is slightly greater for the zero-charge model than for the full-charge model, i.e., the polypeptide's electrostatic interactions, which include hydrogen bonds, slightly destabilize the helix. The helix stabilizing contribution of the hydrophobic effect was found to be identical to that of the van der Waals interactions in vacuo (i.e., 1.9 kcal/mol per residue). The zero-charge model has nearly identical helix stability in vacuo and in water, the almost identical free energies of transfer of helix and coil state of the zero-charge oligomer from vacuum to water are found to be small. Thus, the results of this systematic variation of the force field afford a meaningful decomposition of the free energies for helix initiation and growth.

Free energy simulations of axial contacts in sickle-cell hemoglobin

Molecular dynamics simulations have been used to investigate the thermodynamic stability of axial contacts in sickle-cell hemoglobin (HbS). Free energy changes were evaluated for the point mutation beta 121 Glu --> Gln in the axial contact region of HbS crystals. The calculations predict a free energy change of-3.6 kcal/mol per contact for the mutation, which is in qualitative agreement with experimental observations of aggravated sickling found in the double mutant Hb D Los Angeles (beta 6 Glu --> Val. beta 121 Glu --> Gln) relative to HbS (beta 6 Glu --> Val). The beta 121 Glu is sequestered in a salt link with beta 17 Lys located on the same polypeptide chain, making the Glu interactions with its surroundings similar in aggregates and individual hemoglobins. Due to this cancellation of the large electrostatic Glu contributions, the weak nonspecific interactions between the Gln and the neighboring polypeptide chain are the main contributing factor to the enhanced aggregation of Hb D Los Angeles relative to HbS. Together with the previous study of the lateral contact [K. Kuczera et al. (1990) Proceedings of the National Academy of Science USA, Vol. 87, pp, 8481-8485], the present results provide a more complete picture of the forces driving the sickling aggregation. A comparison of different treatments of internal flexibility in free energy simulations and analysis of rate of convergence of the different calculated properties has also been performed.

Contributions of asparagine at alpha 97 to the cooperative oxygenation process of hemoglobin

According to the X-ray crystallographic results from human deoxyhemoglobin, beta 99Asp at the alpha 1 Beta 2 interface forms hydrogen bonds with alpha 42Tyr and alpha 97Asn. To clarify the structural and functional roles of the hydrogen bond between alpha 97Asn and beta 99Asp, we have engineered a recombinant hemoglobin in which alpha 97Asn is replaced by Ala, and have investigated its oxygen-binding properties, and have used proton nuclear magnetic resonance spectroscopy to determine the structural consequences of the mutation. Recombinant Hb (alpha 97Asn-->Ala) shows a milder alteration of functional properties compared to the severely impaired beta 99 mutants of the human abnormal hemoglobins. The addition of inositol hexaphosphate, an allosteric effector, causes recovery of the functional properties of recombinant Hb (alpha 97 Asn-->Ala) almost to the level of human normal adult hemoglobin without this allosteric effector. r Hb (alpha 97 Asn-->Ala) shows very similar tertiary structure around the heme pockets and quaternary structure in the alpha 1 beta 2 interface compared to those of human normal adult hemoglobin. The proton nuclear magnetic resonance spectrum of the deoxy form of this recombinant hemoglobin shows the existence of an altered hydrogen bond which is believed to be between alpha 42Tyr and beta 99Asp at the alpha 1 beta 2 interface. Thus, the present results suggest that the intersubunit hydrogen bond between alpha 97 Asn and beta 99Asp at the alpha 1 beta 2 interface is not as crucial as the one between alpha 42Tyr and beta 99Asp in the deoxy quaternary structure. Preliminary molecular dynamics simulations have been used to calculate the contributions of specific interactions of several amino acid residues in r Hb (alpha 97Asn-->Ala) to the free energy of cooperativity of this recombinant hemoglobin. The results of these calculations are consistent with the experimental results.

Proline cis-trans isomerization in staphylococcal nuclease: multi-substrate free energy perturbation calculations

Staphylococcal nuclease A exists in two folded forms that differ in the isomerization state of the Lys 116-Pro 117 peptide bond. The dominant form (90% occupancy) adopts a cis peptide bond, which is observed in the crystal structure. NMR studies show that the relatively small difference in free energy between the cis and trans forms (delta Gcis-->trans approximately 1.2 kcal/mol) results from large and nearly compensating differences in enthalpy and entropy (delta Hcis-->trans approximately delta TScis-->trans approximately 10 kcal/mol). There is evidence from X-ray crystal structures that the structural differences between the cis and the trans forms of nuclease are confined to the conformation of residues 112-117, a solvated protein loop. Here, we obtain a thermodynamic and structural description of the conformational equilibrium of this protein loop through an exhaustive conformational search that identified several substates followed by free energy simulations between the substrates. By partitioning the search into conformational substates, we overcame the multiple minima problem in this particular case and obtained precise and reproducible free energy values. The protein and water environment was implicitly modeled by appropriately chosen nonbonded terms between the explicitly treated loop and the rest of the protein. These simulations correctly predicted a small free energy difference between the cis and trans forms composed of larger, compensating differences in enthalpy and entropy. The structural predictions of these simulations were qualitatively consistent with known X-ray structures of nuclease variants and yield a model of the unknown minor trans conformation.

Free energy simulations: the meaning of the individual contributions from a component analysis

A theoretical analysis is made of the decomposition into contributions from individual interactions of the free energy calculated by thermodynamic integration. It is demonstrated that such a decomposition, often referred to as "component analysis," is meaningful, even though it is a function of the integration path. Moreover, it is shown that the path dependence can be used to determine the relation of the contribution of a given interaction to the state of the system. To illustrate these conclusions, a simple transformation (Cl- to Br- in aqueous solution) is analyzed by use of the Reference Interaction Site Model-Hypernetted Chain Closure integral equation approach; it avoids the calculational difficulties of macromolecular simulation while retaining their conceptual complexity. The difference in the solvation free energy between chloride and bromide is calculated, and the contributions of the Lennard-Jones and electrostatic terms in the potential function are analyzed by the use of suitably chosen integration paths. The model is also used to examine the path dependence of individual contributions to the double free energy differences (delta delta G or delta delta A) that are often employed in free energy simulations of biological systems. The alchemical path, as contrasted with the experimental path, is shown to be appropriate for interpreting the effects of mutations on ligand binding and protein stability. The formulation is used to obtain a better understanding of the success of the Poisson-Boltzmann continuum approach for determining the solvation properties of polar and ionic systems.

Helix-capping interaction in lambda Cro protein: a free energy simulation analysis

The stability mutant Tyr-26-->Asp was studied in the Cro protein from bacteriophage lambda using free energy molecular dynamics simulations. The mutant was calculated to be more stable than the wild type by 3.0 +/- 1.7 kcal/mol/monomer, in reasonable agreement with experiment (1.4 kcal/mol/monomer). Moreover, the aspartic acid in the mutant was found to form a capping interaction with the amino terminus of the third alpha-helix of Cro. The simulations were analyzed to understand better the source of the stability of this helix-capping interaction and to examine the results in light of previous explanations of stabilizing helix caps--namely, a model of local unsatisfied hydrogen bonds at the helix termini and the helix macrodipole model. Analysis of the simulations shows that the stabilizing effect of this charged helical cap is due both to favorable hydrogen bonds with backbone NH groups at the helix terminus and to favorable electrostatic interactions (but not hydrogen bonds) with their carbonyls (effectively the next row of local dipoles in the helix). However, electrostatic interactions are weak or negligible with backbone dipolar groups in the helix further away from the terminus. Moreover, the importance of other local electrostatic interactions with polar side chains near the helix terminus, which are neglected in most treatments of this effect, are shown to be important. Thus, the results support a model that is intermediate between the two previous explanations: both unsatisfied hydrogen bonds at the helix terminus and other, local preoriented dipolar groups stabilize the helix cap. These findings suggest that similar interactions with preoriented dipolar groups may be important for cooperativity in other charge-dipole interactions and may be employed to advantage for molecular design.

Molecular dynamics study of the binding of phenylalanine stereoisomers to thermolysin

The stereospecificity in binding of phenylalanine as inhibitor in the active site of the thermolysin, has been investigated by means of molecular dynamics simulations using free energy integration techniques. The difference in the free energy of binding was found to be 2.0 +/- 1.8 kJ/mol in favour of the D-form. This agrees with the experimental value, 2.8 kJ/mol. The result was obtained using a standard empirical force field (that of GROMOS). A different force field with 30% bigger charges (more like ab initio charges) was also tried. This resulted in less fluctuations and a more precise binding, but in a free energy difference that was clearly larger than the experimental one. The phenylalanine backbone is located close to the zinc atom and the ring stays in the hydrophobic pocket in both the cases. The two stereoisomers differ mainly in the orientation of the backbone plane with respect to the active site and the rotational state of the dihedral around the C alpha-C beta bond.

A free-energy simulation study of a bacterial collagenase inhibitor

The cis-trans conformational isomerisation of the N-methylated peptide bond of a bacterial collagenases peptide inhibitor (HS-CH2-CH2-CO-Pro-NMe-Ala) has been investigated by molecular dynamics simulations, energy minimisations and free-energy simulations in presence of solvent molecules. The free-energy difference between the cis and trans forms obtained by the thermodynamics integration method is equal to 0.95 kcal/mol in favour of the trans form, in accord with the experimental result. The main contribution to this free-energy difference comes from solute-water electrostatic interactions. Interestingly, we show that the number of interactions between water molecules and the oxygen atoms of the inhibitor is larger in the trans form than in the cis. Thus the organisation of water molecules around the inhibitor appears crucial in determining the population of the cis and trans conformers.

Computational chemistry and molecular modeling of electron-transfer proteins

The methods of computational chemistry and molecular modeling are becoming more and more accessible to biochemists with the advent of fast, inexpensive graphics workstations and well-tested computer programs. The state of the art in small molecules allows chemists to use these programs as "black boxes" and be confident of the results at an amazingly high level of precision. This is not the case, however, for biological macro-molecules at this time. Therefore, it is necessary before using the programs listed in Section I that we familiarize ourselves with their theoretical basis and limitations. It is also important that they be used in the context of the accumulated literature on their use. The survey given in this chapter is intended as an introduction to these tools and as a source for initiating the discovery process with the literature cited.

Helix propensity of Ala and Val: a free energy perturbation study

Difference in helix propensity between Ala and Val was examined by a molecular dynamics/free energy perturbation method. Simulations were based on a simple two state model of helix-coil transition. Val10 in an Ala based 17mer peptide (Y1K2A3A4A5A6K7A8A9V10A11K12A13A14A15A16K17) used as the alpha-helical model, or in an extended trimer (A9V10A11) used as the random coil model, was perturbed to Ala by a slow growth method. The computed delta delta G (-1.27 +/- 0.92 kcal/mol) reproduced semiquantitatively the experimental delta delta G (-0.7 kcal/mol; Padmanabhan et al., Nature 344 (1990) 268). Inclusion of intraperturbed contributions was essential. Free energy component analysis showed that intra-perturbed interaction within Va10 contributed positive value to delta delta G and that interaction between Lys7 and Val10 contributed negative value. The latter dominated the former, which resulted in the larger helix propensity of Ala.

Molecular dynamics study of structure and stability of a model coiled coil

This paper employs methods used earlier to study helix propensity in a model alpha-helix. The methods are extended to simulations of a motif structure of the alpha-helical coiled coil, i.e., a structure with a simple amino acid sequence, containing only alanine, leucine, and valine, with leucine and valine forming hydrophobic contacts in the helix interface (positions "d" and "a"). Dynamic simulations of the model coiled-coil structure reproduce characteristic features of the coiled-coil motif seen in experimental studies. Free energy simulations were used to assess the change in stability of the model when a leucine pair or a valine pair in the helix interface was replaced with an alanine pair. A leucine pair at position d was found to contribute 3.4 kcal/mol to the stability of the coiled coil relative to an alanine pair, and a valine pair at position a was found to contribute 0.8 kcal/mol relative to an alanine pair. The value for the leucine pair agrees with reports in two experimental studies with molecules having different amino sequence. The value for the valine pair is reasonable given the smaller size of the valine side chain and the intrinsic low helix propensity of valine. No experimental value was available for comparison.

Binding of an antiviral agent to a sensitive and a resistant human rhinovirus. Computer simulation studies with sampling of amino acid side-chain conformation. I. Mapping the rotamers of residue 188 of viral protein 1

The mutation of valine 188 to leucine in the viral protein 1 of human rhinovirus 14 renders the virus resistant to certain antiviral compounds. Thermodynamic-cycle perturbation theory provides a means of calculating the difference in the binding free energies of an antiviral compound to the wild-type virus and to the mutant virus. In calculating the relevant free-energy differences in molecular dynamics simulations, it is important to sample the multiple rotational isomers of residue 188 correctly. In general, these rotamers will not be fully sampled during a single molecular dynamics simulation. However, the contributions of all the rotamers to the free-energy differences associated with mutation of residue 188 may be considered explicitly once they have been identified and their relative free energies determined. Therefore, we describe here the mapping of the rotamers of residue 188 by steric-bump search and energy minimization techniques, and by the computation of potentials of mean force (p.m.f.s.) using umbrella sampling. The usefulness, validity and efficiency of these methods of examining rotameric states is discussed. Adiabatic mapping by energy minimization was found to be unreliable for this residue due to the small magnitude of its interactions with the surrounding protein atoms. Ambiguities in the adiabatic maps were resolved by computing p.m.f.s. The p.m.f. for valine 188 in the unliganded wild-type virus shows a minimum corresponding to the crystallographically observed conformation of valine 188. The p.m.f.s. for valine 188 in the liganded virus and for leucine 188 in the unliganded mutant virus suggest that the experimentally observed conformations may be interpreted as averages of a number of conformations corresponding to those at the minima in the p.m.f.s. The calculations suggest also that the conformation of leucine 188 may change when the ligand binds. The use of the calculated p.m.f.s. to compute the difference in the free energy of binding of an antiviral compound to the wild-type and mutant rhinoviruses is described in the accompanying article.

Binding of an antiviral agent to a sensitive and a resistant human rhinovirus. Computer simulation studies with sampling of amino acid side-chain conformations. II. Calculation of free-energy differences by thermodynamic integration

Thermodynamic-cycle perturbation theory and molecular dynamics simulations were used to calculate the difference in the free energy of binding of the antiviral compound WIN53338 to the wild-type human rhinovirus 14 and to a drug-resistant mutant of the virus in which valine 188 of the viral protein 1 is mutated to leucine. Because of the difficulty of achieving adequate sampling of all of the rotational isomers of amino acid side-chains in molecular dynamics simulations, an explicit treatment of the effects of the existence of multiple rotational isomers of residue 188 on the calculated free energies was used. The rotamers of residue 188 were first mapped by steric and energetic techniques as described in the accompanying article. Thermodynamic integration was then carried out during simulations of the virus, both with and without the antiviral compound bound, by mutating residue 188 while restraining its side-chain to one conformation. The contributions of the other rotamers of residue 188 to the free-energy changes for this mutation were then added to those calculated by thermodynamic integration as correction factors. Binding of WIN53338 to the wild-type virus was calculated to be favored over binding to the mutant virus by 1.7(+/- 3.0) kcal/mol. This is consistent with experimental data which, if differences in activity are assumed to be due to differences in binding, indicate that the binding affinity of WIN53338 for the wild-type virus is at least 0.15 to 1.7 kcal/mol greater than for the mutant virus. Thermodynamic integration was also performed in the conventional manner without restraints and was found to give less accurate results.