Standard molecular mechanics (MM) force fields predict a nearly linear decrease in hydration free energy with each successive addition of a methyl group to ammonia or acetamide, whereas a nonadditive relationship is observed experimentally. In contrast, the non-additive hydration behavior is reproduced directly using a quantum mechanics (QM)/MM-based free-energy perturbation (FEP) method wherein the solute partial atomic charges are updated at every window. Decomposing the free energies into electrostatic and van der Waals contributions and comparing the results with the corresponding free energies obtained using a conventional FEP method and a QM/MM method wherein the charges are not updated suggests that inaccuracies in the electrostatic free energies are the primary reason for the inability of the conventional FEP method to predict the experimental findings. The QM/MM-based FEP method was subsequently used to evaluate inhibitors of the diabetes drug target fructose-1,6-bisphosphatase adenosine 5'-monophosphate and 6-methylamino purine riboside 5'-monophosphate. The predicted relative binding free energy was consistent with the experimental findings, whereas the relative binding free energy predicted using the conventional FEP method differed from the experimental finding by an amount consistent with the overestimated relative solvation free energies calculated for alkylamines. Accordingly, the QM/MM-based FEP method offers potential advantages over conventional FEP methods, including greater accuracy and reduced user input. Moreover, since drug candidates often contain either functionality that is inadequately treated by MM (e.g., simple alkylamines and alkylamides) or new molecular scaffolds that require time-consuming development of MM parameters, these advantages could enable future automation of FEP calculations as well as greatly increase the use and impact of FEP calculations in drug discovery.

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. © 2006 Wiley Periodicals, Inc. J Comput Chem 28: 491–494, 2007

The structures of the mirror image (+)and (—)-m/n.s-a/ir/-adductsof 7,8-dihydroxj-9JO-epox) -7,8,9, lÃ1⁄4-tetrahydroben/o(a)p>rene lo guaninc A/2have been of great interest because the high biological activity of 7,8dih) droxy-9,10-epo\y-7,8,9,10-tetrahydrobenzo(a)py rene in mammalian mutagenesis and tumorigcnesis has been attributed to the predominant (+)-fca/ii-anfi-adduct. \Ve have carried out new potential energy mini mization studies, involving wide-scale «informational searches on small modified DNA subunits, followed by energy-minimized build-up tech niques, to generate atomic resolution views of these adducts. These energy-minimi/ed 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 (+)adduci 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 pyrcnyl long axis is directed more toward the 3' end of the modified strand, and with much less helix axis bend. Structures with intercalation character com puted for these adducts are less preferred. The favored computed struc tures agree with spectroscopic data on the (+)and (-)-trans-antiadducts, whereas recent experimental evidence suggests that m-adducts assume intercalation-type structures. Perhaps the «informationaldistinc tions elucidated for the (+)and (—)-rra/iÃ--anÃ-i-adducts play a role in their differential tumorigenic properties in mammalian systems.

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.

Solvation free energy is an important molecular characteristic useful in drug discovery because it represents the desolvation cost of a ligand binding to a receptor. Most of the recent developments in the estimation of solvation free energy require the use of molecular mechanics and dynamics calculations. Group contribution methods have been rarely used in the past for calculating salvation free energy because automated prediction methods have not been developed in this regard. As an aid to combinatorial library design, we explored rapid and accurate means of computing salvation free energies from the covalent structures of organic molecules and compared the results on a test set with the GB/SA solvation model.. Two independent additive-constitutive QSPR methods have been developed for the computation of solvation free energy. The first is a QSPR model (HLOGS) derived using a technique that uses the counts of distinct/similar fragments and substructures for each molecule as variables in a PLS regression. The second method (ALOGS) uses an extensive atom classification scheme developed earlier for the calculation of Log P. A database of 265 molecules with experimentally determined salvation free energies is used to derive the HLOGS (r = 0.97; rms = 0.58) and ALOGS (r = 0.98; rms = 0.38) models, which were then tested on 27 molecules not present in the training set. A detailed comparison of the HLOGS, ALOGS, GB/SA (with AMBER* and OPLSA* potentials) on the test set showed that the HLOGS and ALOGS models give better results than the GB/SA model. Among the three methods tested, the ALOGS method gives the best result on the test set (r = 0.96; rms = 0.86), though the parametrization for this method is incomplete as many atom types are undetermined due to their absence in the current training set. The HLOGS method appears to handle intramolecular interactions better than the ALOGS method.

Three-dimensional structures of Dendrotoxin (DtX), Toxin-I (DpI), and Toxin-K (DpK) were determined using molecular mechanics and molecular dynamics techniques. The overall molecular conformation and protein folding of the three dendrotoxins are very similar to the published crystal structures of bovine pancreatic trypsin inhibitor (BPTI) and alpha-DtX. Major secondary structural regions of the dendrotoxins are stable without much fluctuation during the dynamics simulation; the regions corresponding to the turns and bends (rich in lysines and arginines) exhibit more fluctuations. The conformational angles and the C alpha...C alpha' distances of the three disulfides (in each of the dendrotoxins) are different from each other. Comparative model building studies, involving the dendrotoxins and the proteinases, reveal that the key interactions (observed in BPTI-trypsin complex) needed for anti-protease activity are absent due to structural differences between the dendrotoxins and BPTI at the anti-protease loop; this explains the inability of the dendrotoxins to inhibit proteinases. The model also suggests that the solvent-exposed beta-turn region, rich in lysines (residues 26-28), might bind directly to the extracellular anionic sites of the receptors (K+ channels) by ionic interactions. The strikingly homologous cysteine distribution (Cys-x-x-x-Cys) in DtX, DpI, and DpK, at the C-terminus, induces the occurrence of a characteristic conformational motif, consisting of an alpha-helix (in an amphiphilic environment) stabilized by two disulfides, one involving a cysteine at the beta-strand, and the other at the N-terminus. This amphiphilic secondary structural element seems to provide the rigid frame work needed for exposing the proposed active site region of the dendrotoxins to the anionic sites of the K+ channel receptors.

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.

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.

Methylphosphonate (MP)-substituted antisense DNA oligomers have significant Potential for genome targeted therapy. The physicochemical basis of the difference in R- versus S-MP diastereomers incorporated into DNA oligomers with regard to complementary DNA target hybridization has been poorly understood. State of the art advanced molecular computational methods and supercomputer technology were applied to identify key physicochemical determinants involved in stabilizing and destabilizing target hybridization. MP-oligomer:DNA target hybridization is more stable with R-MP substitution and is due to favorable hydrophobic interactions between the equatorial projecting methyl group and water molecules. S-MP destabilizes the double strand helix formation by less favorable local hydrophobic interactions with waters, which result in DNA helix unwinding, and by promoting local changes from C2' endo to C3' endo in the 5' furanose ring. In the oligomer studied, R-MP thymine is more stable than R-MP uracil, resulting from hydrophobic interaction to locally stabilize the helix due to the presence of the pyrimidine C5 methyl group. The presence of the C5 methyl group in thymine also appears to influence helical winding and the local conformation of the furanose ring. These numerical simulations are supported by experimental data, enabling us to propose several mechanisms which underlie some of the experimental observations. These results support further development of diastereomerically pure MP analogues as diagnostic/therapeutic agents involving sequence-specific DNA interactions.

Significant advancements in the treatment of cancer have developed during the last four decades as a result of dedicated clinical and experimental efforts. Many of the therapeutic gains are related to discovery and development of more effective medicinal agents, technologic advancements, and an enhanced understanding of the fundamental chemical and biologic interactions involving the pathogenesis and pathophysiology of these heterogenous diseases. It is now possible to cure several types of malignancy and to achieve significant palliation in a variety of other tumors. Unfortunately, many of the more common types of neoplasms (e.g. carcinomas of the lung, breast, GI tract, and melanoma) are refractory to therapy with currently available agents. A most significant obstacle to address in the coming years is multiple drug resistance in which the cytotoxic action of pharmacologic agents is rendered ineffectual by a transmembrane pump in tumor cells (1–3). Early in this past decade, infection with the human immunodeficiency virus (HIV) has presented a major health problem because of the lethal nature of the disease, significant latent interval between infection and disease manifestation, fluctuating/ evolving epidemiologic patterns, and the lack of effective therapy. An important complication of pharmacologic agents that must be considered in the development of new therapeutic agents for the AIDS and neoplastic disorders are the immediate and long term clinical toxicities that frequently affect the quality of life (4). Because of the immediate and life-threatening nature of these diseases the untoward effects of these agents have been monitored and managed expectantly, since some toxicities can be lethal to the patient. Our primary research is directed towards generating new classes of pharmacologic agents possessing highly specific cytotoxicity for neoplastic cells, and cells that have incorporated the HIV genome. The major goals of developing such agents will be to cure these diseases which are currently refractory to therapy with minimal patient toxicity.

Diels-Alder and nitrile oxide intramolecular cycloadditions were studied using several methods. The structures found using all methods are similar when the forming bonds lengths are constrained, but the stereochemical predictions are quite different. The experimental stereochemical differences found for the parent Diels-Alder reactions forming 6-5 and 6-6 systems are rationalized. When the addends are linked by three methylene groups (formation of a five-membered ring), the strain in the transition structure (TS) causes the addends to twist about the forming bonds, resulting in a skewed TS as compared to the intermolecular TS. However, when the addends are linked by four methylene groups (formation of a six-membered ring), there is little strain in the TS, and the addends do not twist.

The free energy perturbation method has been employed to determine the binding free energy contributions of different groups of two classes of HIV-1 proteinase inhibitors: (1) a hydroxyethylene isostere inhibitor, Ala-Ala-Phe[CH(OH)-CH2]Gly-Val-Val-OMe (reported by Dreyer et al. 1), and a reduced peptide inhibitor, MVT-101 (reported by Miller et al. 2). For the first inhibitor, the configuration of the central hydroxyl group is changed from S to R in two steps. In the first step, the hydroxyl group in the S configuration was mutated to a hydrogen, and in the second step, the hydroxyl group of the (R)-OH analogue of the inhibitor was mutated to a hydrogen. In this way the binding contributions of the hydroxyl group in the two diastereomers are determined separately in addition to obtaining the effect of changing the hydroxyl group configuration from S to R. The calculated free energy difference between the binding of the two diastereomers is 3.37 +/- 0.64 kcal/mol, which is close the experimental value of 2.6 kcal/mol. The calculations on the substitution of Gly by Nle at the P'1 position of the same inhibitor predict an enhancement in the binding by about 1.7 kcal/mol. Similar calculations on the substitution of Nle with Met in MVT-101 inhibitor predict a decrease in binding by about 0.7 kcal/mol. The details of these results will be discussed and compared with the results of a similar study on pepstatin-rhizopus pepsin complex.

Substitution of 5-methylcytosine for cytosine in DNA oligomers increases the stability of the Hoogsteen-paired third strand hybridizing to complementary double strand DNA under physiological pH conditions. The physicochemical mechanism(s) underlying the increased stability of DNA triplex formation which accompanies 5-methylcytosine substitution in the third DNA strand is poorly understood. To address these objectives, we performed ab initio quantum mechanical and statistical mechanical studies on the equilibrium geometries and solution proton affinities of 5-methylcytosine and cytosine and incorporated these into large-scale numerical simulation models of a DNA triple helical system d[CT]10-d[GA]10-d[5mC+T]10. The purpose of these molecular simulations was to develop accurate models of cytosine and 5-methylcytosine nucleosides and to incorporate these structural and proton affinity models in thermodynamic and structural studies of DNA oligomer hybridization. We find that the models correctly.represent the net proton affinities of cytosine and 5-methylcytosine in water and that the calculated hybridization free energy difference (DELTA-DELTA-G(hybrid) = 13.5 kcal/mol) between 5-methylcytosine and cytosine-substituted triple helices is qualitatively accurate. Using experimental entropy data (DELTA-S-degrees), we predict that 5-methylcytosine substitution into DNA oligomers stabilizes the net transition enthalpy (DELTA-H-degrees) of the triple helix by 73.2 kcal/mol, and a net DELTA-H-degrees(MEC) 7.3 kcal/mol per base triplet relative to the cytosine-substituted oligomer in the DNA triple helix. We observe no major conformational differences in 5-methylcytosine versus cytosine-substituted triple helices during molecular dynamics. These studies provide a greater understanding of key physicochemical mechanism(s) and properties of DNA oligomers incorporated into DNA triple helices containing 5-methylcytosine versus cytosine involved in modulating the stability of hybridization.

The ab initio quantum mechanics, molecular mechanics, and free energy perturbation methods have been applied to study the energetics of the active site of Rhizopus pepsin and its interactions with several inhibitors derived from pepstatin. The studies on the Asp diad in the active site of the enzyme show that the energetics of the diad are very sensitive to small changes in the relative orientations of the diad and hence the energetic equivalence of the two charge states of the diad (arising due to the protonation of either of the two aspartates) can be easily attained by small changes in the atomic positions of the diad. Further, the studies point out that the proton possibly shuttles between the two inner oxygens of the diad. The barrier for the proton shuttle could be as low as 1.0 kcal/mol when the inner oxygen distance is around 2.5 angstrom and it increases with the increase in this distance. Although the present studies show that the configurations of the Asp diad distorted from planarity are lower in energy than the coplanar configuration found in the crystal structure, the latter configuration is found to be crucial for optimal inhibitor binding. This is also borne out in the calculated binding free energy differences between pepstatin and its derivatives. The calculated values obtained with the lower energy configuration of the Asp diad were found to be lower than those obtained with the Asp diad configuration found in the crystal structure, and the latter values were closer to the experimental results. For the mutation of the central statine residue of pepstatin to dehydroxystatine, the calculated free energy difference of 5.17 kcal/mol is in good agreement with the experimental value. This shows that the contribution of about 5 kcal/mol to binding from the hydroxyl group of the central statine residue is mainly due to the strong interaction of this group with the negatively charged Asp diad. The results of the other mutations on pepstatin are also in support of this view.

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.

A free energy perturbation study of solvation in hydrazine and carbon tetrachloride is carried out to examine the process of solvation of different solutes in these two solvents. For this purpose, models of liquid hydrazine and liquid carbon tetrachloride were generated by molecular dynamics simulations. The structure of liquid hydrazine obtained from the molecular dynamics simulation is discussed in detail. Differences in the free energies of solvation of different solutes belonging to different classes of ions and molecules have been determined. The solutes studied include closed shell ions, tetraalkylammonium ions, normal alkanes, and tetraalkylmethane molecules. The calculated differences in free energy of solvation, DELTA-G, between two different solutes in these solvents compare well with the experimental values. Detailed analysis of the solvation behavior in these solvents and their comparison with the behavior of solvation in water suggest that solvation behavior in hydrazine resembles that in water for many solutes, whereas the behavior in carbon tetrachloride is different. The results of this study support the view that the special phenomenon observed in the hydration of apolar solutes is a result of the structural peculiarity of liquid water.

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.