Calculation of the relative change in binding free energy of a protein-inhibitor complex.

By means of a thermodynamic perturbation method implemented with molecular dynamics, the relative free energy of binding was calculated for the enzyme thermolysin complexed with a pair of phosphonamidate and phosphonate ester inhibitors. The calculated difference in free energy of binding was 4.21 +/- 0.54 kilocalories per mole. This compares well with the experimental value of 4.1 kilocalories per mole. The method is general and can be used to determine a change or "mutation" in any system that can be suitably represented. It is likely to prove useful for protein and drug design.

Free energy calculations by computer simulation.

A fundamental problem in chemistry and biochemistry is understanding the role of solvation in determining molecular properties. Recent advances in statistical mechanical theory and molecular dynamics methodology can be used to solve this problem with the aid of supercomputers. By using these advances the free energies of solvation of all the chemical classes of amino acid side chains, four nucleic acid bases and other organic molecules can be calculated. The effect of a site-specific mutation on the stability of trypsin is predicted. The results are in good agreement with available experiments.

Structures of the (+)- and (-)-trans-7,8-dihydroxy-anti-9,10-epoxy-7,8,9,10-tetrahydrobenzo(a)pyre ne adducts to guanine-N2 in a duplex dodecamer.

The structures of the mirror image (+)- and (-)-trans-anti-adducts of 7,8-dihydroxy-9,10-epoxy-7,8,9,10-tetrahydrobenzo(a)pyrene to guanine N2 have been of great interest because the high biological activity of 7,8-dihydroxy-9,10-epoxy-7,8,9,10-tetrahydrobenzo(a)pyrene in mammalian mutagenesis and tumorigenesis has been attributed to the predominant (+)-trans-anti-adduct. We have carried out new potential energy minimization studies, involving wide-scale conformational searches on small modified DNA subunits, followed by energy-minimized build-up techniques, to generate atomic resolution views of these adducts. These energy-minimized duplex dodecamers were then subjected to 100-ps molecular dynamic simulations with solvent and salt to yield animated molecular structures. The most favored computed structure for the (+)-adduct places the pyrenyl moiety in the B-DNA minor groove, with its long axis directed toward the 5' end of the modified strand, and with a pronounced bend in the helix axis. In the (-)-adduct, there are 2 favored structures. One places the pyrenyl moiety in the minor groove, whereas the other positions it in the major groove; in both cases, the pyrenyl long axis is directed more toward the 3' end of the modified strand, and with much less helix axis bend. Structures with intercalation character computed for these adducts are less preferred. The favored computed structures agree with spectroscopic data on the (+)- and (-)-trans-anti-adducts, whereas recent experimental evidence suggests that cis-adducts assume intercalation-type structures. Perhaps the conformational distinctions elucidated for the (+)- and (-)-trans anti-adducts play a role in their differential tumorigenic properties in mammalian systems.

Systematic variation in bond lengths of peptides.

The lengths of the carbonyl as well as of the adjacent C-N and C-C bonds in peptides are shown to vary systematically with the central C-N bond length. Results of ab initio calculations on N-methylacetamide and its Li+, Na+ and Mg2+ complexes are also discussed.

Protein-ligand dynamics. A 96 picosecond simulation of a myoglobin-xenon complex.

A 96 picosecond dynamics trajectory of myoglobin with five xenon-probe ligands in internal cavities is examined to study the effect of protein motions on ligand motion and internal cavity fluctuations. Average structural and energetic properties indicate that the simulation is well behaved. The average protein volume is similar to the volume of the X-ray model and the main-chain atom root-mean-square deviation between the X-ray model and the average dynamical structure is 1.25 A. The protein volume oscillates 3 to 4% around the volume of the X-ray structure. These fluctuations lead to changes in the internal free volume and in the size, shape and location of atom-sized cavity features. Transient cavities produced in the simulation have a crucial role in the movement of two of the ligands. One of the ligands escapes to the protein surface, whilst a second ligand travels through the protein interior. Complex gating processes involving several protein residues are responsible for producing the necessary pores through which the ligand passes between transient cavities or packing defects.

Computational studies of the interaction of myoglobin and xenon.

Computational studies are used to investigate the energies of xenon binding to myoglobin and to describe pathways through the protein interior for a metmyoglobin-xenon complex. Empirical energy calculations indicate a favorable enthalpic contribution of 0.6 to 4.2 kcal/mol to xenon binding for four experimentally determined xenon sites. These calculated enthalpies help to explain the different xenon occupancies observed experimentally. A fifth site, modeled in place of the iron co-ordinated water molecule in the distal cavity, is also predicted to bind xenon. The largest contribution to the binding energy is from van der Waals' interactions with smaller contributions from polarization and protein strain terms. Ligand trajectory calculations as well as a new geometric algorithm define a connecting network of channel-like pathways through the static protein structure. One or two pathways appear to lead most easily from each major internal cavity to the protein surface. The importance of these channels in protein dynamics and their implications as routes for ligand motion are discussed.

Molecular mechanics simulations on covalent complexes between anthramycin and B DNA.

We present molecular mechanics simulation of the covalent interactions of the potent antitumor antibiotic belonging to the pyrrolo[1,4]benzodiazepine class, anthramycin, with six deoxydecanucleotides, d(GCGCGCGCGC)2, d(G10) X d(C10), d(GCGCGTGCGC) X d(GCGCACGCGC), d(GCGCGAGCGC) X d(GCGCTCGCGC), d(GGGGGAGGGG) X d(CCCCTCCCCC), and d(GGGGGTGGGG) X d(CCCCACCCCC), in their minor grooves. The complexes are characterized by both a network of hydrogen bonds between the drug and the polynucleotide and good packing interactions. The DNA double helix in these complexes shows very minimal distortion, and interactions of the drug with the decanucleotides seem to be not very sensitive to the sequence variation around the site of complex formation. The conformational features in the complexes obtained are generally consistent with the experimentally derived conclusions by recent NMR and 2-D NOE studies.

Conformations of complexes between mitomycin and decanucleotides. 2. Application of the model to mitomycin C derivatives. Extension to covalent binding with adenine.

Molecular mechanics simulation of the interactions of important mitomycin C analogues monocovalently bound to DNA models are presented. These analogues included substituents such as p-hydroxyphenyl, 2-mercaptoethyl, and dimethylamidinium on N7 of mitomycin C and the DNA models consisted of d(GCGCGCGCGC)2 and d(GCGCATGCGC)2. The excellent fits and strong binding affinities of these highly potent analogues support the usefulness of the model. The binding of a mitomycin-related N-phenylpyrrole with a carbamoyloxy substituent to 06 of guanine was studied. Finally, a reactive mitomycin intermediate proposed by Moore was shown to interact with DNA in a way consistent with the formation of a covalent adduct.

Modeling complex molecular interactions involving proteins and DNA.

We have presented a perspective of progress in three areas of simulations of complex molecules: the development of force fields for molecular simulation; the application of computer graphics, molecular mechanics and molecular dynamics in simulations of DNA and DNA-drug complexes and the application of computer graphics, molecular mechanics and quantum mechanics in studies of enzyme substrate interactions. It is our perspective that improvements are being made in force fields, and these will allow a more accurate simulation of structures and energies of complex molecules. In the area of DNA molecular mechanics and dynamics, it is clear that the use of computer graphics model building combined with NMR NOE data is a potentially very powerful tool in accurately determining structures of drug-DNA complexes using molecular mechanics and dynamics. Finally, we are in a position to reasonably simulate structures and (qualitatively) energies for complete reaction pathways of enzymes using a combination of computer graphics, molecular mechanics and quantum mechanics. More accurate energies and pathways are sure to follow, using the combined molecular mechanics/quantum mechanics optimization developed by Singh and the free energy perturbation methods pioneered in Groningen and Houston.

On loop folding in nucleic acid hairpin-type structures.

In a series of studies, combining NMR, optical melting and T-jump experiments, it was found that DNA hairpins display a maximum stability when the loop part of the molecule comprises four or five nucleotide residues. This is in contrast with the current notion based on RNA hairpin studies, from which it had been established that a maximum hairpin stability is obtained for six or seven residues in the loop. Here we present a structural model to rationalize these observations. This model is based on the notion that to a major extent base stacking interactions determine the stability of nucleic acid conformations. The model predicts that loop folding in RNA is characterized by an extension of the base stacking at the 5'-side of the double helix by five or six bases; the remaining gap can then easily be closed by two nucleotides. Conversely, loop folding in DNA is characterized by extending base stacking at the 3'-side of the double helical stem by two or three residues; again bridging of the remaining gap can then be achieved by one or two nucleotides. As an example of loop folding in RNA the anticodon loop of yeast tRNAPhe is discussed. For the DNA hairpin formed by d(ATCCTAT4TAGGAT) it is shown that the loop structure obtained from molecular mechanics calculations obeys the above worded loop folding principles.

A combined 2D-NMR and molecular dynamics analysis of the structure of the actinomycin D: d(ATGCAT)2 complex.

We present a comparative analysis of an NMR experiment and molecular and harmonic dynamics simulations of an actinomycin D: d(ATGCAT)2 complex. A comparison of NOE measurements and 1/R6 weighted proton-proton distances confirm the general correctness of the Actinomycin D-DNA model proposed by Sobell. There are, however, some substantial differences between the proton-proton distances inferred from the NOE results and the molecular and harmonic dynamics simulations. The remaining discrepancies could either come from contributions of other conformations to the average properties of the complex or from uncertainties in the NMR distance analysis. An analysis of the molecular dynamics helix properties, sugar puckers, hydrogen bonding, rms fluctuations and torsional properties are qualitatively consistent with those from previous simulations, but the presence of an intercalated drug leads to some new structural and dynamical features.

Ab initio quantum mechanics-based free energy perturbation method for calculating relative solvation free energies.

A free energy perturbation (FEP) method was developed that uses ab initio quantum mechanics (QM) for treating the solute molecules and molecular mechanics (MM) for treating the surroundings. Like our earlier results using AM1 semi empirical QMs, the ab initio QM/MM-based FEP method was shown to accurately calculate relative solvation free energies for a diverse set of small molecules that differ significantly in structure, aromaticity, hydrogen bonding potential, and electron density. Accuracy was similar to or better than conventional FEP methods. The QM/MM-based methods eliminate the need for time-consuming development of MM force field parameters, which are frequently required for drug-like molecules containing structural motifs not adequately described by MM. Future automation of the method and parallelization of the code for Linux 128/256/512 clusters is expected to enhance the speed and increase its use for drug design and lead optimization.

Probing the salt bridge in the dihydrofolate reductase-methotrexate complex by using the coordinate-coupled free-energy perturbation method.

The importance of the ionic interaction due to the formation of the salt bridge between the Asp-27 and the pteridine ring in Escherichia coli dihydrofolate reductase-methotrexate complex has been studied by using the free-energy perturbation method. The calculation suggests that the ion-pair contribution to the binding energy is insignificant, as the enzyme surroundings do not stabilize the salt bridge to the extent of the desolvation of the charged groups. The activation barrier for the proton exchange between the pteridine ring and the Asp-27 is calculated to be 20.1 kcal/mol (1 cal = 4.184 J) by using the coordinate-coupled perturbation method, implying that this may be a channel to the proton exchange from the pteridine ring to the solvent. The Gibbs-energy difference of binding between the Asn-27 and Ser-27 is calculated to be 3.2 kcal/mol and is mainly due to the electrostatic interactions.

A free-energy perturbation study of the binding of methotrexate to mutants of dihydrofolate reductase.

The importance of hydrophobic residues to the binding of methotrexate in the active site of dihydrofolate reductase (EC 1.5.1.3) was examined by a free-energy perturbation method. The replacement of a strictly conserved residue, Phe-31, by tyrosine or valine costs 1.8 and 5.1 kcal/mol, respectively, to the binding of the drug (1 cal = 4.184 J). In the case of the Phe31----Tyr mutation, the loss of the binding energy is due to the desolvation of the phenolic group; in the case of Phe31----Val mutation, it is mainly due to the loss of the interaction with the drug. The replacement of Leu-54 by glycine decreases the binding energy by 4.0 kcal/mol. A calculation on the mutation of Phe-31 to serine shows that the alteration could reduce the binding energy of methotrexate by 9.7 kcal/mol. The calculations clearly show that the hydrophobic interactions are as important as the hydrophilic ones in the binding of methotrexate.

Development of a quantum mechanics-based free-energy perturbation method: use in the calculation of relative solvation free energies.

Free-energy perturbation (FEP) is considered the most accurate computational method for calculating relative solvation and binding free-energy differences. Despite some success in applying FEP methods to both drug design and lead optimization, FEP calculations are rarely used in the pharmaceutical industry. One factor limiting the use of FEP is its low throughput, which is attributed in part to the dependence of conventional methods on the user's ability to develop accurate molecular mechanics (MM) force field parameters for individual drug candidates and the time required to complete the process. In an attempt to find an FEP method that could eventually be automated, we developed a method that uses quantum mechanics (QM) for treating the solute, MM for treating the solute surroundings, and the FEP method for computing free-energy differences. The thread technique was used in all transformations and proved to be essential for the successful completion of the calculations. Relative solvation free energies for 10 structurally diverse molecular pairs were calculated, and the results were in close agreement with both the calculated results generated by conventional FEP methods and the experimentally derived values. While considerably more CPU demanding than conventional FEP methods, this method (QM/MM-based FEP) alleviates the need for development of molecule-specific MM force field parameters and therefore may enable future automation of FEP-based calculations. Moreover, calculation accuracy should be improved over conventional methods, especially for calculations reliant on MM parameters derived in the absence of experimental data.

Molecular dynamics simulations of d(C-G-C-G-A) X d(T-C-G-C-G) with and without "hydrated" counterions.

We present the results of molecular dynamics simulations on d(C-G-C-G-A) X d(T-C-G-C-G) with fully charged phosphates with and without inclusion of counterions. The average structures found in the two simulations are similar, but the simulation with counterions does give an average helix repeat, tilt, and twist in better agreement with those found in the x-ray structure of d(C-G-C-G-A-A-T-T-C-G-C-G)2. The average sugar pucker phases and amplitudes are in qualitative agreement with those found in NMR studies of double-helical DNA, and a number of examples of sugar repuckering from C2' endo to C3' endo carbon conformations in the sugar ring are found. The hydrogen bond correlations as well as torsion correlations are analyzed, and some interesting long-range correlations between dihedral angles are found.

A molecular dynamics simulation of double-helical B-DNA including counterions and water.

We present the results of an atomic level molecular dynamical simulation of a 5-base-pair fragment of double-helical DNA with inclusion of water and sodium counterions and a complete description of their electrostatic interactions. The shape of the double helix is preserved throughout the simulation, and the helix repeat is calculated to be 10.0, in reasonable agreement with experimental results. The most flexible conformational angles in the structure are the glycosidic angle and the sugar pucker.