Molecular modeling of organic and biomolecular systems using BOSS and MCPRO

A remote arene-binding site on prostate specific membrane antigen revealed by antibody-recruiting small molecules

Prostate specific membrane antigen (PSMA) is a membrane-bound glutamate carboxypeptidase overexpressed in many forms of prostate cancer. Our laboratory has recently disclosed a class of small molecules, called ARM-Ps (antibody-recruiting molecule targeting prostate cancer) that are capable of enhancing antibody-mediated immune recognition of prostate cancer cells. Interestingly, during the course of these studies, we found ARM-Ps to exhibit extraordinarily high potencies toward PSMA, compared to previously reported inhibitors. Here, we report in-depth biochemical, crystallographic, and computational investigations which elucidate the origin of the observed affinity enhancement. These studies reveal a previously unreported arene-binding site on PSMA, which we believe participates in an aromatic stacking interaction with ARMs. Although this site is composed of only a few amino acid residues, it drastically enhances small molecule binding affinity. These results provide critical insights into the design of PSMA-targeted small molecules for prostate cancer diagnosis and treatment; more broadly, the presence of similar arene-binding sites throughout the proteome could prove widely enabling in the optimization of small molecule-protein interactions.

E/Z Energetics for Molecular Modeling and Design

Thermochemical data have been obtained from G3B3 quantum mechanical calculations for 18 prototypical organic molecules, which exhibit E/Z conformational equilibria. The results are fundamentally important for molecular design including evaluation of structures from protein-ligand docking. For the 18 E/Z pairs, relative energies, enthalpies, free energies, and dipole moments are reported; the E - Z free-energy differences at 298 K range from +8.2 kcal/mol for 1,3-dimethyl carbamate to -6.4 kcal/mol for acetone oxime. A combination of steric and electronic effects can rationalize the variations. Free energies of hydration were also estimated using the GB/SA continuum solvent model. These results indicate that differential hydration is unlikely to qualitatively change the preferred direction of the E/Z equilibria, though further study with free-energy methods using explicit solvent is desirable.

Computer tools in the discovery of HIV-1 integrase inhibitors

Computer-aided drug design (CADD) methodologies have made great advances and contributed significantly to the discovery and/or optimization of many clinically used drugs in recent years. CADD tools have likewise been applied to the discovery of inhibitors of HIV-1 integrase, a difficult and worthwhile target for the development of efficient anti-HIV drugs. This article reviews the application of CADD tools, including pharmacophore search, quantitative structure-activity relationships, model building of integrase complexed with viral DNA and quantum-chemical studies in the discovery of HIV-1 integrase inhibitors. Different structurally diverse integrase inhibitors have been identified by, or with significant help from, various CADD tools.

Exploring solvent effects upon the Menshutkin reaction using a polarizable force field

The energetics of the Menshutkin reaction between triethylamine and ethyl iodide have been computed using B3LYP and MP2 with the LANL2DZ, LANL2DZd, SVP, MIDI!, 6-311G(d,p), and aug-cc-PVTZ basis sets. Small- and large-core energy-consistent relativistic pseudopotentials were employed. Solvent effect corrections were computed from QM/MM Monte Carlo simulations utilizing free-energy perturbation theory, PDDG/PM3, and both a nonpolarizable OPLS and polarizable OPLS-AAP force field. The B3LYP/MIDI! theory level provided the best DeltaG(++) values with a mean absolute error (MAE) of 4.9 kcal/mol from experiment in cyclohexane, CCl(4), THF, DMSO, acetonitrile, water, and methanol. However, the relative rates in cyclohexane, and to a certain extent CCl(4), were determined to be greatly underestimated when using the nonpolarizable OPLS force field. An overall reduction in the MAE to 3.1 kcal/mol using B3LYP/MIDI!/OPLS-AAP demonstrated the need for a fully polarizable force field when computing solvent effects for highly dipolar transition structures in low-dielectric media. The MAEs obtained with PDDG/PM3/OPLS and OPLS-AAP of 5.3 and 3.8 kcal/mol, respectively, provided comparable results to B3LYP at a fraction of the computational resources. The large rate accelerations observed in the reaction were correlated to an increased stabilization of the emerging charge separation at the transition state via favorable solute-solvent interactions.

Thorpe-Ingold acceleration of oxirane formation is mostly a solvent effect

The Thorpe-Ingold hypothesis for the gem-dimethyl effect in the cyclization reactions of 2-chloroethoxide derivatives has been investigated computationally in the gas phase and in aqueous solution. Ab initio MP2/6-311+G(d,p) and CBS-Q calculations reveal little intrinsic difference in reactivity with increasing alpha-methylation for the series of reactants 1-3. However, inclusion of continuum hydration or of explicit hydration through mixed quantum and statistical mechanics (MC/FEP) simulations does reproduce the substantial, experimentally observed rate increases with increasing alpha-methylation. Analysis of the MC/FEP results provides clear evidence that the rate increases stem primarily from increased steric hindrance to hydration of the nucleophilic oxygen atom with increasing alpha-methylation. Thus, the gem-dimethyl acceleration of oxirane formation for 1-3 is found to be predominantly a solvent effect.

Quantum mechanical/molecular mechanical modeling finds Diels-Alder reactions are accelerated less on the surface of water than in water

Quantum and molecular mechanics calculations for the Diels-Alder reactions of cyclopentadiene with 1,4-naphthoquinone, methyl vinyl ketone, and acrylonitrile have been carried out at the vacuum-water interface and in the gas phase. In conjunction with previous studies of these cycloadditions in dilute solution, a more complete picture of aqueous environmental effects emerges with implications for the origin of observed rate accelerations using heterogeneous aqueous suspensions, "on water" conditions. The pure TIP4P water slab maintains the bulk density and hydrogen-bonding properties in central water layers. The bulk region merges to vacuum over a ca. 5 A band with progressive diminution of the density and hydrogen bonding. The relative free energies of activation and transition structures for the reactions at the interface are found to be intermediate between those calculated in the gas phase and in bulk water; i.e., for the reaction with 1,4-naphthoquinone, the DeltaDeltaG(++) values relative to the gas phase are -3.6 and -7.3 kcal/mol at the interface and in bulk water, respectively. Thus, the results do not support the notion that a water surface is more effective than bulk water for catalysis of such pericyclic reactions. The trend is in qualitative agreement with expectations based on density considerations and estimates of experimental rate constants for the gas phase, a heterogeneous aqueous suspension, and a dilute aqueous solution for the reaction of cyclopentadiene with methyl vinyl ketone. Computed energy pair distributions reveal a uniform loss of 0.5-1.0 hydrogen bond for the reactants and transition states in progressing from bulk water to the vacuum-water interface. Orientational effects are apparent at the surface; e.g., the carbonyl group in the methyl vinyl ketone transition structure is preferentially oriented into the surface. Also, the transition structure for the 1,4-naphthoquinone case is buried more in the surface, and the free energy of activation for this reaction is most similar to the result in bulk water.

Eastern extension of azoles as non-nucleoside inhibitors of HIV-1 reverse transcriptase; cyano group alternatives

Design of non-nucleoside inhibitors of HIV-1 reverse transcriptase is being pursued with the assistance of free energy perturbation (FEP) calculations to predict relative free energies of binding. Extension of azole-containing inhibitors into an 'eastern' channel between Phe227 and Pro236 has led to the discovery of potent and structurally novel derivatives.

Discovery of novel fibroblast growth factor receptor 1 kinase inhibitors by structure-based virtual screening

Fibroblast growth factors (FGFs) play important roles in embryonic development, angiogenesis, wound healing, and cell proliferation and differentiation. In search of inhibitors of FGFR1 kinase, 2.2 million compounds were docked into the ATP binding site of the protein. A co-crystal structure, which shows two alternative conformations for the nucleotide binding loop, is reported. Docking was performed on both conformations and, ultimately, 23 diverse compounds were purchased and assayed. Following hit validation, two compounds 10 and 16, a benzylidene derivative of pseudothiohydantoin and a thienopyrimidinone derivative, respectively, were discovered that inhibit FGFR1 kinase with IC(50) values of 23 and 50 microM. Initial optimization of 16 led to the more unsaturated 40, which has significantly enhanced potency, 1.9 microM. The core structures represent new structural motifs for FGFR1 kinase inhibitors. The study also illustrates complexities associated with the choice of protein structures for docking, possible use of multiple kinase structures to seek selectivity, and hit identification.

Role of water in the multifaceted catalytic antibody 4B2 for allylic isomerization and Kemp elimination reactions

Specificity toward a single reaction is a well-known characteristic of catalytic antibodies. However, contrary to convention, catalytic antibody 4B2 possesses the ability to efficiently catalyze two unrelated reactions: a Kemp elimination and an allylic isomerization of a beta,gamma-unsaturated ketone. To elucidate how this multifaceted antibody operates, mixed quantum and molecular mechanics calculations coupled to Monte Carlo simulations were carried out. The antibody was determined to derive its adaptability for the mechanistically different reactions through the rearrangement of water molecules in the active site into advantageous geometric orientations for enhanced electrostatic stabilization. In the case of the Kemp elimination, a general base, Glu L34, carried out the proton abstraction from the isoxazole ring of 5-nitro-benzisoxazole while water molecules delivered specific stabilization at the transition state. The role of water was found to be more pronounced in the allylic isomerization because the solvent actively participated in the stepwise mechanism. A rate-limiting abstraction of the alpha-proton from the beta,gamma-unsaturated ketone via Glu L34 led to the formation of a neutral dienol intermediate, which was rapidly reprotonated at the gamma-position via a solvent hydronium ion. Preferential channeling of H(3)O(+) in the active site ensured a stereoselective proton exchange from the alpha- to the gamma-position, in good agreement with deuterium exchange NMR and HPLC experiments. Ideas for improved water-mediated catalytic antibody designs are presented. In a technical advancement, improvements to a recent polynomial fitting and integration technique utilizing free energy perturbation theory delivered greater accuracy and speed gains.

Advances in quantum and molecular mechanical (QM/MM) simulations for organic and enzymatic reactions

Application of combined quantum and molecular mechanical (QM/MM) methods focuses on predicting activation barriers and the structures of stationary points for organic and enzymatic reactions. Characterization of the factors that stabilize transition structures in solution and in enzyme active sites provides a basis for design and optimization of catalysts. Continued technological advances allowed for expansion from prototypical cases to mechanistic studies featuring detailed enzyme and condensed-phase environments with full integration of the QM calculations and configurational sampling. This required improved algorithms featuring fast QM methods, advances in computing changes in free energies including free-energy perturbation (FEP) calculations, and enhanced configurational sampling. In particular, the present Account highlights development of the PDDG/PM3 semi-empirical QM method, computation of multi-dimensional potentials of mean force (PMF), incorporation of on-the-fly QM in Monte Carlo (MC) simulations, and a polynomial quadrature method for efficient modeling of proton-transfer reactions. The utility of this QM/MM/MC/FEP methodology is illustrated for a variety of organic reactions including substitution, decarboxylation, elimination, and pericyclic reactions. A comparison to experimental kinetic results on medium effects has verified the accuracy of the QM/MM approach in the full range of solvents from hydrocarbons to water to ionic liquids. Corresponding results from ab initio and density functional theory (DFT) methods with continuum-based treatments of solvation reveal deficiencies, particularly for protic solvents. Also summarized in this Account are three specific QM/MM applications to biomolecular systems: (1) a recent study that clarified the mechanism for the reaction of 2-pyrone derivatives catalyzed by macrophomate synthase as a tandem Michael-aldol sequence rather than a Diels-Alder reaction, (2) elucidation of the mechanism of action of fatty acid amide hydrolase (FAAH), an unusual Ser-Ser-Lys proteolytic enzyme, and (3) the construction of enzymes for Kemp elimination of 5-nitrobenzisoxazole that highlights the utility of QM/MM in the design of artificial enzymes.

In silico prediction of drug solubility: 4. Will simple potentials suffice?

In view of the extreme importance of reliable computational prediction of aqueous drug solubility, we have established a Monte Carlo simulation procedure which appears, in principle, to yield reliable solubilities even for complex drug molecules. A theory based on judicious application of linear response and mean field approximations has been found to reproduce the computationally demanding free energy determinations by simulation while at the same time offering mechanistic insight. The focus here is on the suitability of the model of both drug and solvent, i.e., the force fields. The optimized potentials for liquid simulations all atom (OPLS-AA) force field, either intact or combined with partial charges determined either by semiempirical AM1/CM1A calculations or taken from the condensed-phase optimized molecular potentials for atomistic simulation studies (COMPASS) force field has been used. The results illustrate the crucial role of the force field in determining drug solubilities. The errors in interaction energies obtained by the simple force fields tested here are still found to be too large for our purpose but if a component of this error is systematic and readily removed by empirical adjustment the results are significantly improved. In fact, consistent use of the OPLS-AA Lennard-Jones force field parameters with partial charges from the COMPASS force field will in this way produce good predictions of amorphous drug solubility within 1 day on a standard desktop PC. This is shown here by the results of extensive new simulations for a total of 47 drug molecules which were also improved by increasing the water box in the hydration simulations from 500 to 2000 water molecules.

Energetics of displacing water molecules from protein binding sites: consequences for ligand optimization

A strategy in drug design is to consider enhancing the affinity of lead molecules with structural modifications that displace water molecules from a protein binding site. Because success of the approach is uncertain, clarification of the associated energetics was sought in cases where similar structural modifications yield qualitatively different outcomes. Specifically, free-energy perturbation calculations were carried out in the context of Monte Carlo statistical mechanics simulations to investigate ligand series that feature displacement of ordered water molecules in the binding sites of scytalone dehydratase, p38-alphaMAP kinase, and EGFR kinase. The change in affinity for a ligand modification is found to correlate with the ease of displacement of the ordered water molecule. However, as in the EGFR example, the binding affinity may diminish if the free-energy increase due to the removal of the bound water molecule is not more than compensated by the additional interactions of the water-displacing moiety. For accurate computation of the effects of ligand modifications, a complete thermodynamic analysis is shown to be needed. It requires identification of the location of water molecules in the protein-ligand interface and evaluation of the free-energy changes associated with their removal and with the introduction of the ligand modification. Direct modification of the ligand in free-energy calculations is likely to trap the ordered molecule and provide misleading guidance for lead optimization.

Prediction of the water content in protein binding sites

An efficient molecular simulation methodology has been developed to determine the positioning of water molecules in the binding site of a protein or protein-ligand complex. Occupancies and absolute binding free energies of water molecules are computed using a statistical thermodynamics approach. The methodology, referred to as Just Add Water Molecules (JAWS), features "theta-water" molecules that can appear and disappear on a binding-site grid. Key approximations render the technique far more efficient than conventional free energy simulations. Testing of JAWS on five diverse examples (neuraminidase, scytalone dehydratase, major urinary protein 1, beta-lactoglobulin, and COX-2) demonstrates its accuracy in locating hydration sites in comparison to results from high-resolution crystal structures. Possible applications include aid in refinement of protein crystal structures, drug lead optimization, setup of docking calculations, and simulations of protein-ligand complexes.

Vancomycin analogs: Seeking improved binding of d-Ala-d-Ala and d-Ala-d-Lac peptides by side-chain and backbone modifications

In order to seek vancomycin analogs with improved performance against VanA and VanB resistant bacterial strains, extensive computational investigations have been performed to examine the effects of side-chain and backbone modifications. Changes in binding affinities for tripeptide cell-wall precursor mimics, Ac(2)-l-Lys-d-Ala-d-Ala (3) and Ac(2)-l-Lys-d-Ala-d-Lac (4), with vancomycin analogs were computed with Monte Carlo/free energy perturbation (MC/FEP) calculations. Replacements of the 3-hydroxyl group in residue 7 with small alkyl or alkoxy groups, which improve contacts with the methyl side chain of the ligands'd-Ala residue, are predicted to be the most promising to enhance binding for both ligands. The previously reported amine backbone modification as in 5 is shown to complement the hydrophobic modifications for binding monoacetylated tripeptides. In addition, replacement of the hydroxyl groups in residues 5 and 7 by fluorine is computed to have negligible impact on binding the tripeptides, though it may be pharmacologically advantageous.

CHARMM: the biomolecular simulation program

CHARMM (Chemistry at HARvard Molecular Mechanics) is a highly versatile and widely used molecular simulation program. It has been developed over the last three decades with a primary focus on molecules of biological interest, including proteins, peptides, lipids, nucleic acids, carbohydrates, and small molecule ligands, as they occur in solution, crystals, and membrane environments. For the study of such systems, the program provides a large suite of computational tools that include numerous conformational and path sampling methods, free energy estimators, molecular minimization, dynamics, and analysis techniques, and model-building capabilities. The CHARMM program is applicable to problems involving a much broader class of many-particle systems. Calculations with CHARMM can be performed using a number of different energy functions and models, from mixed quantum mechanical-molecular mechanical force fields, to all-atom classical potential energy functions with explicit solvent and various boundary conditions, to implicit solvent and membrane models. The program has been ported to numerous platforms in both serial and parallel architectures. This article provides an overview of the program as it exists today with an emphasis on developments since the publication of the original CHARMM article in 1983.

In Silico Improvement of beta3-peptide inhibitors of p53 x hDM2 and p53 x hDMX

There is great interest in molecules capable of inhibiting the interactions between p53 and its negative regulators hDM2 and hDMX, as these molecules have validated potential against cancers in which one or both oncoproteins are overexpressed. We reported previously that appropriately substituted beta(3)-peptides inhibit these interactions and, more recently, that minimally cationic beta(3)-peptides are sufficiently cell permeable to upregulate p53-dependent genes in live cells. These observations, coupled with the known stability of beta-peptides in a cellular environment, and the recently reported structures of hDM2 and hDMX, motivated us to exploit computational modeling to identify beta-peptides with improved potency and/or selectivity. This exercise successfully identified a new beta(3)-peptide, beta53-16, that possesses the highly desirable attribute of high affinity for both hDM2 and hDMX and identifies the 3,4-dichlorophenyl moiety as a novel determinant of hDMX affinity.

Efficient drug lead discovery and optimization

During the 1980s, advances in the abilities to perform computer simulations of chemical and biomolecular systems and to calculate free energy changes led to the expectation that such methodology would soon show great utility for guiding molecular design. Important potential applications included design of selective receptors, catalysts, and regulators of biological function including enzyme inhibitors. This time also saw the rise of high-throughput screening and combinatorial chemistry along with complementary computational methods for de novo design and virtual screening including docking. These technologies appeared poised to deliver diverse lead compounds for any biological target. As with many technological advances, realization of the expectations required significant additional effort and time. However, as summarized here, striking success has now been achieved for computer-aided drug lead generation and optimization. De novo design using both molecular growing and docking are illustrated for lead generation, and lead optimization features free energy perturbation calculations in conjunction with Monte Carlo statistical mechanics simulations for protein-inhibitor complexes in aqueous solution. The specific applications are to the discovery of non-nucleoside inhibitors of HIV reverse transcriptase (HIV-RT) and inhibitors of the binding of the proinflammatory cytokine MIF to its receptor CD74. A standard protocol is presented that includes scans for possible additions of small substituents to a molecular core, interchange of heterocycles, and focused optimization of substituents at one site. Initial leads with activities at low-micromolar concentrations have been advanced rapidly to low-nanomolar inhibitors.

Statistical thermodynamics of biomembranes

An overview of the major issues involved in the statistical thermodynamic treatment of phospholipid membranes at the atomistic level is summarized: thermodynamic ensembles, initial configuration (or the physical system being modeled), force field representation as well as the representation of long-range interactions. This is followed by a description of the various ways that the simulated ensembles can be analyzed: area of the lipid, mass density profiles, radial distribution functions (RDFs), water orientation profile, deuterium order parameter, free energy profiles and void (pore) formation; with particular focus on the results obtained from our recent molecular dynamic (MD) simulations of phospholipids interacting with dimethylsulfoxide (Me(2)SO), a commonly used cryoprotective agent (CPA).

Catalytic mechanism and performance of computationally designed enzymes for Kemp elimination

A series of enzymes for Kemp elimination of 5-nitrobenzisoxazole has been recently designed and tested. In conjunction with the design process, extensive computational analyses were carried out to evaluate the potential performance of four of the designs, as presented here. The enzyme-catalyzed reactions were modeled using mixed quantum and molecular mechanics (QM/MM) calculations in the context of Monte Carlo (MC) statistical mechanics simulations. Free-energy perturbation (FEP) calculations were used to characterize the free-energy surfaces for the catalyzed reactions as well as for reference processes in water. The simulations yielded detailed information about the catalytic mechanisms, activation barriers, and structural evolution of the active sites over the course of the reactions. The catalytic mechanism for the designed enzymes KE07, KE10(V131N), and KE15 was found to be concerted with proton transfer, generally more advanced in the transition state than breaking of the isoxazolyl N-O bond. On the basis of the free-energy results, all three enzymes were anticipated to be active. Ideas for further improvement of the enzyme designs also emerged. On the technical side, the synergy of parallel QM/MM and experimental efforts in the design of artificial enzymes is well illustrated.

Origin of the activity drop with the E50D variant of catalytic antibody 34E4 for Kemp elimination

In enzymes, multiple structural effects cooperatively lead to the high catalytic activity, while individually these effects can be small. The design of artificial enzymes requires the understanding and ability to manipulate such subtle effects. The 34E4 catalytic antibody, catalyzing Kemp elimination of 5-nitrobenzisoxazole, and its Glu50Asp (E50D) variant are the subject of the present investigation. This removal of only a methylene group yields an approximately 30-fold reduction in the rate for the catalyzed Kemp elimination. Here, the aim is to understand this difference in the catalytic performance. The mechanism of Kemp elimination catalyzed by 34E4 and the E50D mutant is elucidated using QM/MM Monte Carlo simulations and free energy perturbation theory. In both proteins, the reaction is shown to follow a single-step, concerted mechanism. In the mutant, the activation barrier rises by 2.4 kcal/mol, which corresponds to a 62-fold rate deceleration, which is in good agreement with the experimental data. The positions and functionality of the residues in the active site are monitored throughout the reaction. It is concluded that the looser contact with the base, shorter base-Asn58 contact, less favorable pi-stacking with Trp91 in the transition state of the reaction, and different solvation pattern all contribute to the reduction of the reaction rate in the E50D variant of 34E4.

Small molecule hydration free energies in explicit solvent: An extensive test of fixed-charge atomistic simulations

Using molecular dynamics free energy simulations with TIP3P explicit solvent, we compute the hydration free energies of 504 neutral small organic molecules and compare them to experiments. We find, first, good general agreement between the simulations and the experiments, with an RMS error of 1.24 kcal/mol over the whole set (i.e., about 2 kT) and a correlation coefficient of 0.89. Second, we use an automated procedure to identify systematic errors for some classes of compounds, and suggest some improvements to the force field. We find that alkyne hydration free energies are particularly poorly predicted due to problems with a Lennard-Jones well depth, and find that an alternate choice for this well depth largely rectifies the situation. Third, we study the non-polar component of hydration free energies - that is, the part that is not due to electrostatics. While we find that repulsive and attractive components of the non-polar part both scale roughly with surface area (or volume) of the solute, the total non-polar free energy does not scale with the solute surface area or volume, because it is a small difference between large components and is dominated by the deviations from the trend. While the methods used here are not new, this is a more extensive test than previous explicit solvent studies, and the size of the test set allows identification of systematic problems with force field parameters for particular classes of compounds. We believe that the computed free energies and components will be valuable to others in future development of force fields and solvation models.

Methods for Monte Carlo simulations of biomacromolecules

The state-of-the-art for Monte Carlo (MC) simulations of biomacromolecules is reviewed. Available methodologies for sampling conformational equilibria and associations of biomacromolecules in the canonical ensemble, given a continuum description of the solvent environment, are reviewed. Detailed sections are provided dealing with the choice of degrees of freedom, the efficiencies of MC algorithms and algorithmic peculiarities, as well as the optimization of simple movesets. The issue of introducing correlations into elementary MC moves, and the applicability of such methods to simulations of biomacromolecules is discussed. A brief discussion of multicanonical methods and an overview of recent simulation work highlighting the potential of MC methods are also provided. It is argued that MC simulations, while underutilized biomacromolecular simulation community, hold promise for simulations of complex systems and phenomena that span multiple length scales, especially when used in conjunction with implicit solvation models or other coarse graining strategies.

Validating CHARMM parameters and exploring charge distribution rules in structure-based drug design

Using an extensive series of TIBO compounds that are non-nucleoside inhibitors of HIV-1 reverse transcriptase, we have systematically evaluated the quality of recently developed ligand parameters that are consistent with the CHARMM22 force field. Thermodynamic integration simulations for 44 pairs of TIBO compounds achieve a high level of success with an overall average unsigned error (AUE) in the relative binding affinities of 1.3 kcal/mol; however, the accuracy is strongly dependent on the size differential between the substituents sampled as well as the class of functional group. Low errors are observed among the alkyl, allyl, aldehyde, nitrile, trifluorinated methyl, and halide TIBO derivatives and large systematic errors among thioether derivatives. We have also investigated how different charge assignment schemes for small molecules impact the quality of computed binding affinities for a subset of this series. This study demonstrates the advantage of using model compounds to derive physically meaningful charge distributions and bond-charge increments for rapidly expanding fragment libraries for drug development applications. Specifically, in the absence of a bond-charge increment for a given pair of atom types, the strategy of adopting CHELPG charges from localized regions of model compounds provides reliable results when modeling with the CHARMM force field.

Vancomycin resistance: modeling backbone variants with D-Ala-D-Ala and D-Ala-D-Lac peptides

To seek vancomycin analogs with broader antibacterial activity, effects of backbone modifications for the agylcon 2 on binding with D-Ala-D-Ala- and D-Ala-D-Lac-containing peptides were investigated by Monte Carlo/free energy perturbation (MC/FEP) calculations. The experimental trend in binding affinities for 2 with three tripeptides was well reproduced. Possible modifications of the peptide bond between residues 4 and 5 were then considered, specifically for conversion of the OCNH linkage to CH(2)NH(2)(+) (6), FCCH (7), HCCH (8), and HNCO (9). The MC/FEP results did not yield binding improvements for 7, 8, and 9, though the fluorovinyl replacement is relatively benign. The previously reported analog 6 remains as the only variant that exhibits improved affinity for the D-Ala-D-Lac sequence and acceptable affinity for the D-Ala-D-Ala sequence.

Characterizing aqueous solution conformations of a peptide backbone using Raman optical activity computations

Mounting spectroscopic evidence indicates that alanine predominantly adopts extended polyproline II (PPII) conformations in short polypeptides. Here we analyze Raman optical activity (ROA) spectra of N-acetylalanine-N'-methylamide (Ala dipeptide) in H2O and D2O using density functional theory on Monte Carlo (MC) sampled geometries to examine the propensity of Ala dipeptide to adopt compact right-handed (alpha(R)) and left-handed (alpha(L)) helical conformations. The computed ROA spectra based on MC-sampled alpha(R) and PPII peptide conformations contain all the key spectral features found in the measured spectra. However, there is no significant similarity between the measured and computed ROA spectra based on the alpha(L)- and beta-conformations sampled by the MC methods. This analysis suggests that Ala dipeptide populates the alpha(R) and PPII conformations but no substantial population of alpha(L)- or beta-structures, despite sampling alpha(L)- and beta-structures in our MC simulations. Thus, ROA spectra combined with the theoretical analysis allow us to determine the dominant populated structures. Including explicit solute-solvent interactions in the theoretical analysis is essential for the success of this approach.

CHARMM-GUI: a web-based graphical user interface for CHARMM

CHARMM is an academic research program used widely for macromolecular mechanics and dynamics with versatile analysis and manipulation tools of atomic coordinates and dynamics trajectories. CHARMM-GUI, http://www.charmm-gui.org, has been developed to provide a web-based graphical user interface to generate various input files and molecular systems to facilitate and standardize the usage of common and advanced simulation techniques in CHARMM. The web environment provides an ideal platform to build and validate a molecular model system in an interactive fashion such that, if a problem is found through visual inspection, one can go back to the previous setup and regenerate the whole system again. In this article, we describe the currently available functional modules of CHARMM-GUI Input Generator that form a basis for the advanced simulation techniques. Future directions of the CHARMM-GUI development project are also discussed briefly together with other features in the CHARMM-GUI website, such as Archive and Movie Gallery.

Optimization of azoles as anti-human immunodeficiency virus agents guided by free-energy calculations

Efficient optimization of an inactive 2-anilinyl-5-benzyloxadiazole core has been guided by free energy perturbation (FEP) calculations to provide potent non-nucleoside inhibitors of human immunodeficiency virus (HIV) reverse transcriptase (NNRTIs). An FEP "chlorine scan" was performed to identify the most promising sites for substitution of aryl hydrogens. This yielded NNRTIs 8 and 10 with activities (EC50) of 820 and 310 nM for protection of human T-cells from infection by wild-type HIV-1. FEP calculations for additional substituent modifications and change of the core heterocycle readily led to oxazoles 28 and 29, which were confirmed as highly potent anti-HIV agents with activities in the 10-20 nM range. The designed compounds were also monitored for possession of desirable pharmacological properties by use of additional computational tools. Overall, the trends predicted by the FEP calculations were well borne out by the assay results. FEP-guided lead optimization is confirmed as a valuable tool for molecular design including drug discovery; chlorine scans are particularly attractive since they are both straightforward to perform and highly informative.

The implementation of a fast and accurate QM/MM potential method in Amber

Version 9 of the Amber simulation programs includes a new semi-empirical hybrid QM/MM functionality. This includes support for implicit solvent (generalized Born) and for periodic explicit solvent simulations using a newly developed QM/MM implementation of the particle mesh Ewald (PME) method. The code provides sufficiently accurate gradients to run constant energy QM/MM MD simulations for many nanoseconds. The link atom approach used for treating the QM/MM boundary shows improved performance, and the user interface has been rewritten to bring the format into line with classical MD simulations. Support is provided for the PM3, PDDG/PM3, PM3CARB1, AM1, MNDO, and PDDG/MNDO semi-empirical Hamiltonians as well as the self-consistent charge density functional tight binding (SCC-DFTB) method. Performance has been improved to the point where using QM/MM, for a QM system of 71 atoms within an explicitly solvated protein using periodic boundaries and PME requires less than twice the cpu time of the corresponding classical simulation.

Perspective on Free-Energy Perturbation Calculations for Chemical Equilibria

An overview is provided on the computation of free energy changes in solution using perturbation theory, overlap sampling, and related approximate methods. As a specific application, extensive results are provided for free energies of hydration of substituted benzenes using the OPLS-AA force field in explicit TIP4P water. For a similar amount of computer time, the double-wide sampling and overlap sampling methods yield very similar results in the free-energy perturbation calculations. With standard protocols, the average statistical uncertainty in computed differences in free energies of hydration is 0.1 - 0.2 kcal/mol. Application of the power-series expansion in the Peierls equation was also tested. Use of the first-order term is generally reliable, while inclusion of the slowly-convergent, second-order fluctuation term causes deterioration in the results for strongly hydrogen-bonded solutes.

Charge-dependent model for many-body polarization, exchange, and dispersion interactions in hybrid quantum mechanical/molecular mechanical calculations

This work explores a new charge-dependent energy model consisting of van der Waals and polarization interactions between the quantum mechanical (QM) and molecular mechanical (MM) regions in a combined QMMM calculation. van der Waals interactions are commonly treated using empirical Lennard-Jones potentials, whose parameters are often chosen based on the QM atom type (e.g., based on hybridization or specific covalent bonding environment). This strategy for determination of QMMM nonbonding interactions becomes tedious to parametrize and lacks robust transferability. Problems occur in the study of chemical reactions where the "atom type" is a complex function of the reaction coordinate. This is particularly problematic for reactions, where atoms or localized functional groups undergo changes in charge state and hybridization. In the present work we propose a new model for nonelectrostatic nonbonded interactions in QMMM calculations that overcomes many of these problems. The model is based on a scaled overlap model for repulsive exchange and attractive dispersion interactions that is a function of atomic charge. The model is chemically significant since it properly correlates atomic size, softness, polarizability, and dispersion terms with minimal one-body parameters that are functions of the atomic charge. Tests of the model are examined for rare-gas interactions with neutral and charged atoms in order to demonstrate improved transferability. The present work provides a new framework for modeling QMMM interactions with improved accuracy and transferability.

Improved interaction potentials for charged residues in proteins

Electrostatic interactions dominate the structure and free energy of biomolecules. To obtain accurate free energies involving charged groups from molecular simulations, OPLS-AA parameters have been reoptimized using Monte Carlo free energy perturbation. New parameters fit a self-consistent, experimental set of hydration free energies for acetate (Asp), propionate (Glu), 4-methylimidazolium (Hip), n-butylammonium (Lys), and n-propylguanidinium (Arg), all resembling charged residue side chains, including beta-carbons. It is shown that OPLS-AA free energies depend critically on the type of water model, TIP4P or TIP3P; i.e., each water model requires specific water-charged molecule interaction potentials. New models (models 1 and 3) are thus described for both water models. Uncertainties in relative free energies of charged residues are approximately 2 kcal/mol with the new parameters, due to variations in system setup (MAEs of ca. 1 kcal/mol) and noise from simulations (ca. 1 kcal/mol). The latter error of approximately 1 kcal/mol contrasts MAEs from standard OPLS-AA of up to 13 kcal/mol for the entire series of charged residues or up to 5 kcal/mol for the cationic series Lys, Arg, and Hip. The new parameters can be used directly in molecular simulations with no modification of neutral residues needed and are envisioned to be particular important in simulations where charged residues change environment.

Comparison of protein force fields for molecular dynamics simulations

In the context of molecular dynamics simulations of proteins, the term "force field" refers to the combination of a mathematical formula and associated parameters that are used to describe the energy of the protein as a function of its atomic coordinates. In this review, we describe the functional forms and parameterization protocols of the widely used biomolecular force fields Amber, CHARMM, GROMOS, and OPLS-AA. We also summarize the ability of various readily available noncommercial molecular dynamics packages to perform simulations using these force fields, as well as to use modern methods for the generation of constant-temperature, constant-pressure ensembles and to treat long-range interactions. Finally, we finish with a discussion of the ability of these force fields to support the modeling of proteins in conjunction with nucleic acids, lipids, carbohydrates, and/or small molecules.

Energetic effects for observed and unobserved HIV-1 reverse transcriptase mutations of residues L100, V106, and Y181 in the presence of nevirapine and efavirenz

The effect of mutations on amino acid residues L100, V106, and Y181 for unbound HIV-1 reverse transcriptase (RT) and RT bound to nevirapine and efavirenz was investigated using Monte Carlo/free energy perturbation calculations. Using both native and bound crystal structures of RT, mutation of the amino acid residues to both those observed and unobserved in patients was carried out. The results of the calculations revealed that the variant that survives in patients dosed with either nevirapine or efavirenz had a more positive Delta Delta G value than other variants that were not observed in patients. These data suggest that the mutation observed in patients is the most effective (the one that binds the drug most weakly) of all possible codon change mutations.

Numerical fitting of molecular properties to Hermite Gaussians

A procedure is presented to fit gridded molecular properties to auxiliary basis sets (ABSs) of Hermite Gaussians, analogous to the density fitting (DF) method (Dunlap; et al. J. Chem. Phys. 1979, 71, 4993). In this procedure, the ab initio calculated properties (density, electrostatic potential, and/or electric field) are fitted via a linear- or nonlinear-least-squares procedure to auxiliary basis sets (ABS). The calculated fitting coefficients from the numerical grids are shown to be more robust than analytic density fitting due to the neglect of the core contributions. The fitting coefficients are tested by calculating intermolecular Coulomb and exchange interactions for a set of dimers. It is shown that the numerical instabilities observed in DF are caused by the attempt of the ABS to fit the core contributions. In addition, this new approach allows us to reduce the number of functions required to obtain an accurate fit. This results in decreased computational cost, which is shown by calculating the Coulomb energy of a 4096 water box in periodic boundary conditions. Using atom centered Hermite Gaussians, this calculation is only 1 order of magnitude slower than conventional atom-centered point charges.

From docking false-positive to active anti-HIV agent

Virtual screening of the Maybridge library of ca. 70 000 compounds was performed using a similarity filter, docking, and molecular mechanics-generalized Born/surface area postprocessing to seek potential non-nucleoside inhibitors of human immunodeficiency virus-1 (HIV-1) reverse transcriptase (NNRTIs). Although known NNRTIs were retrieved well, purchase and assaying of representative, top-scoring compounds from the library failed to yield any active anti-HIV agents. However, the highest-ranked library compound, oxadiazole 1, was pursued as a potential "near-miss" with the BOMB program to seek constructive modifications. Subsequent synthesis and assaying of several polychloro-analogs did yield anti-HIV agents with EC50 values as low as 310 nM. The study demonstrates that it is possible to learn from a formally unsuccessful virtual-screening exercise and, with the aid of computational analyses, to efficiently evolve a false positive into a true active.

Polarization Effects for Hydrogen-Bonded Complexes of Substituted Phenols with Water and Chloride Ion

Variations in hydrogen-bond strengths are investigated for complexes of nine para-substituted phenols (XPhOH) with a water molecule and chloride ion. Results from ab initio HF/6-311+G(d, p) and MP2/6-311+G(d, p)//HF/6-311+G(d, p) calculations are compared with those from the OPLS-AA and OPLS/CM1A force fields. In the OPLS-AA model, the partial charges on the hydroxyl group of phenol are not affected by the choice of para substituent, while the use of CM1A charges in the OPLS/CM1A approach does provide charge redistribution. The ab initio calculations reveal a 2.0-kcal/mol range in hydrogen-bond strengths for the XPhOH⋯OH(2) complexes in the order X = NO(2) > CN > CF(3) > Cl > F > H >OH >CH(3) > NH(2). The pattern is not well-reproduced with OPLS-AA, which also compresses the variation to 0.7 kcal/mol. However, the OPLS/CM1A results are in good accord with the ab initio findings for both the ordering and range, 2.3 kcal/mol. The hydrogen bonding is, of course, weaker with XPhOH as acceptor, the order for X is largely inverted, and the range is reduced to ca. 1.0 kcal/mol. The substituent effects are found to be much greater for the chloride ion complexes with a range of 11 kcal/mol. For quantitative treatment of such strong ion-molecule interactions the need for fully polarizable force fields is demonstrated.

Conformation of alkanes in the gas phase and pure liquids

Monte Carlo (MC) statistical mechanics simulations have been carried out for the homologous alkane series of n-butane through n-dodecane in the gas phase and for the pure liquids at 298 K and 1 atm using the OPLS-AA force field. The study addresses potential cumulative deviations of computed properties and potential conformational differences between the gas phase and pure liquids, for example, from self-solvation in the gas phase. The average errors in comparison with experimental data for the computed densities and heats of vaporization are modest at 0.7% and 6.9%, respectively. Also, the invariant gas and liquid-phase results for average end-to-end distances and percentages of trans conformations for each nonterminal C-C bond assert that the conformer populations are not altered upon transfer from the gas phase to the pure liquid for the n-alkanes in this size range. Average end-to-end distances were also computed from the results of conformational searches and corroborated the MC findings. Quantitatively, the OPLS-AA result for the trans population of the C3-C4 bond in n-undecane is in close agreement with the findings from (13)C NMR experiments. Finally, previous work on determining the shortest n-alkane that does not have an all-trans global energy minimum has been extended. The smallest n-alkane with a hairpin geometry that is lower in energy than the all-trans conformer occurs for C(22)H(46) with OPLS-AA, though with a correction for GG sequences, the true turning point is likely in the C(16)-C(18) range.