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.

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.

Ab initio quantum mechanical and X-ray crystallographic studies of gemcitabine and 2'-deoxy cytosine.

Gemcitabine 2',2'-difluoro 2'-deoxy cytosine (GEM) is a novel nucleoside which has demonstrated broad preclinical anti-cancer activity and appears promising in early stage human clinical trials. One purpose of this study was to characterize the energetically favored conformational modes of GEM by means of ab initio quantum mechanical studies with comparison to a novel X-ray crystallographic structure, and to determine the performance of ab initio quantum mechanical theory by comparison with X-ray structural data for GEM and 2'-deoxy cytosine (CYT). Another objective of this study was to attempt to determine key structural and electronic atomic interactions relating to the 2',2'-difluoro substitution in GEM by the application of ab initio quantum mechanical methods. To our knowledge, these are the first reported ab initio quantum mechanical geometry optimizations of nucleosides using large (e.g. 6-31G*) slit valence function basis sets. The development of accurate physicochemical models on a small scale enables us to extend our studies of GEM to more complex studies including DNA incorporation, deamination, ribonucleotide reductase inhibition, and triphosphorylation.

Structural changes by sulfoxidation of phenothiazine drugs.

The side-chain conformations of psychoactive phenothiazine drugs in crystals are different from those of biologically inactive ring sulfoxide metabolites. This study examines the potential energies, molecular conformations and electrostatic potentials in chlorpromazine, levomepromazine (methotrimeprazine), their sulfoxide metabolites and methoxypromazine. The purpose of the study was to examine the significance of the different crystal conformations of active and inactive phenothiazine derivatives, and to determine why phenothiazine drugs lose most of their biological activity by sulfoxidation. Quantum mechanics and molecular mechanics calculations demonstrated that conformations with the side chain folded over the ring structure had lowest potential energy in vacuo, both in the drugs and in the sulfoxide metabolites. In the sulfoxides, side chain conformations corresponding to the crystal structure of chlorpromazine sulfoxide were characterized by stronger negative electrostatic potentials around the ring system than in the parent drugs. This may weaken the electrostatic interaction of sulfoxide metabolites with negatively charged domains in dopamine receptors, and cause the sulfoxides to be virtually inactive in dopamine receptor binding and related pharmacological tests.

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.

Can oligonucleoside methylphosphonates form a stable triplet with a double DNA helix?

Molecular dynamics simulations of oligonucleoside methylphosphonates (MP) and naturally occurring oligodeoxynucleotides (ODN) as the third strand to double-stranded (ds) DNA targets predict that the third strand with MP backbone is more favorable to triple helix formation than the native phosphodiester ODN third strand. In contrast to experimental data, the calculated DNA conformations in both fully solvated triple helical system are found to be in hybrid A and B conformations. The calculations predict that a third strand with MP backbone is readily accommodated in the major groove of the ds DNA target, and adopts a different conformation from a helical triplet formed with native ODN as the third strand through the Hoogsteen base pairing scheme.

500-picosecond molecular dynamics in water of the Man alpha 1----2Man alpha glycosidic linkage present in Asn-linked oligomannose-type structures on glycoproteins.

Molecular dynamics simulations of the Man alpha 1----2Man alpha glycosidic linkage found in the N-linked glycans of glycoproteins were performed in vacuo and in the presence of water. In the latter case significant dampening of the molecular fluctuations was found when compared to the in vacuo simulation. A 500-ps dynamics simulation in water showed only occasional short-lived deviations from the minimum-energy conformation, more consistent with carbohydrate "breathing" than flexibility. These studies add further evidence that oligosaccharides can maintain "fixed" geometries with relatively long lifetimes and are in agreement with experimental NMR-derived parameters for the same linkage in oligomannose structures.

Identification of local determinants of DNA interstrand crosslink formation by cyclophosphamide metabolites.

The anti-cancer drug cyclophosphamide produces cytotoxic effects by DNA interstrand crosslink formation. Current knowledge of local physicochemical factors in DNA influencing crosslink formation in DNA has been based on physical models of Watson-Crick DNA. Specific interactions between DNA and cyclophosphamide metabolites determining formation and stabilization of interstrand crosslinks have been identified by advanced molecular computational methods. These results predict that the more favorable DNA sequence for interstrand crosslinking is 5'-GC-3' rather than the previously proposed 5'-CG-3' and that thymine methyl groups adjacent to guanine will inhibit interstrand crosslinking. In addition the conformational and energetic consequences of formamido-pyrimidine adducts and local water interactions with the drug-DNA complex were identified. These simulations of interstrand crosslinking provide a fundamental basis for understanding the mechanism of action of nitrogen mustard-type alkylating agents, and a rationale for the prospective design of more effective anti-tumor agents.

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.

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.

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.

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.

Free energy perturbation calculations on binding and catalysis after mutating Asn 155 in subtilisin.

Site-directed mutagenesis is a very powerful approach to altering the biological functions of proteins, the structural stability of proteins and the interactions of proteins with other molecules. Several experimental studies in recent years have been directed at estimating the changes in catalytic properties, (rates of binding and catalysis) in site-directed mutants of enzymes compared to the native enzymes. Simulation approaches to the study of complex molecules have also become more powerful, in no small measure owing to the increase in computer power. These simulations have often allowed results of experiments to be rationalized and understood mechanistically. A new approach called the free-energy pertubation method, which uses statistical mechanics and molecular dynamics can often be used for quantitative calculation of free energy differences. We have applied such a technique to calculate the differential free energy of binding and free energy of activation for catalysis of a tripeptide substrate by native subtilisin and a subtilisin mutant (Asn 155----Ala 155). Our studies lead to a calculated difference in free energy of binding which is relatively small, but a calculated change in free energy of catalysis which is substantial. These energies are very close to those determined experimentally (J. A. Wells and D. A. Estell, personal communication), which were not known to us until the simulations were completed. This demonstrates the predictive power and utility of theoretical simulation methods in studies of the effects of site-specific mutagenesis on both enzyme binding and catalysis.

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.

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.

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.

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 mechanical studies of d(CGTACG)2: complex of triostin A with the middle A - T base pairs in either Hoogsteen or Watson-Crick pairing.

Computer graphics model building and molecular mechanical calculations have been carried out on d(CGTACG)2 and its bis-intercalation complexes with triostin A and an N-Me-Ala analogue of triostin A. Two conformations of the DNA have been considered both for the uncomplexed and for complexed nucleic acid: in one the central A - T base pairs are Watson-Crick base paired; in the other they are Hoogsteen base paired. The calculations offer a clear explanation why Hoogsteen base pairing is not favorable in isolated A + T-rich DNA and also suggest reasons why the bis-intercalation of triostin A might help stabilize the neighboring A - T base pairs into a Hoogsteen form. To our knowledge, this is the first study to use molecular mechanical and dynamical methods to investigate the bis-intercalator-DNA complex.