Standard Molecular Mechanics Force Field Research Articles

Ligand-Binding Affinity Estimates Supported by Quantum-Mechanical Methods.

One of the largest challenges of computational chemistry is calculation of accurate free energies for the binding of a small molecule to a biological macromolecule, which has immense implications in drug development. It is well-known that standard molecular-mechanics force fields used in most such calculations have a limited accuracy. Therefore, there has been a great interest in improving the estimates using quantum-mechanical (QM) methods. We review here approaches involving explicit QM energies to calculate binding affinities, with an emphasis on the methods, rather than on specific applications. Many different QM methods have been employed, ranging from semiempirical QM calculations, via density-functional theory, to strict coupled-cluster calculations. Dispersion and other empirical corrections are mandatory for the approximate methods, as well as large basis sets for the stricter methods. QM has been used for the ligand, for a few crucial groups around the ligand, for all the closest atoms (200-1000 atoms), or for the full receptor-ligand complex, but it is likely that with a proper embedding it might be enough to include all groups within ∼6 Å of the ligand. Approaches involving minimized structures, simulations of the end states of the binding reaction, or full free-energy simulations have been tested.

General Model for Treating Short-Range Electrostatic Penetration in a Molecular Mechanics Force Field.

Classical molecular mechanics force fields typically model interatomic electrostatic interactions with point charges or multipole expansions, which can fail for atoms in close contact due to the lack of a description of penetration effects between their electron clouds. These short-range penetration effects can be significant and are essential for accurate modeling of intermolecular interactions. In this work we report parametrization of an empirical charge–charge function previously reported (PiquemalJ.-P.; J. Phys. Chem. A2003, 107, 1035326313624) to correct for the missing penetration term in standard molecular mechanics force fields. For this purpose, we have developed a database (S101×7) of 101 unique molecular dimers, each at 7 different intermolecular distances. Electrostatic, induction/polarization, repulsion, and dispersion energies, as well as the total interaction energy for each complex in the database are calculated using the SAPT2+ method (ParkerT. M.; J. Chem. Phys.2014, 140, 09410624606352). This empirical penetration model significantly improves agreement between point multipole and quantum mechanical electrostatic energies across the set of dimers and distances, while using only a limited set of parameters for each chemical element. Given the simplicity and effectiveness of the model, we expect the electrostatic penetration correction will become a standard component of future molecular mechanics force fields.

Electrostatics of proteins in dielectric solvent continua. I. An accurate and efficient reaction field description

We present a reaction field (RF) method which accurately solves the Poisson equation for proteins embedded in dielectric solvent continua at a computational effort comparable to that of an electrostatics calculation with polarizable molecular mechanics (MM) force fields. The method combines an approach originally suggested by Egwolf and Tavan [J. Chem. Phys. 118, 2039 (2003)] with concepts generalizing the Born solution [Z. Phys. 1, 45 (1920)] for a solvated ion. First, we derive an exact representation according to which the sources of the RF potential and energy are inducible atomic anti-polarization densities and atomic shielding charge distributions. Modeling these atomic densities by Gaussians leads to an approximate representation. Here, the strengths of the Gaussian shielding charge distributions are directly given in terms of the static partial charges as defined, e.g., by standard MM force fields for the various atom types, whereas the strengths of the Gaussian anti-polarization densities are calculated by a self-consistency iteration. The atomic volumes are also described by Gaussians. To account for covalently overlapping atoms, their effective volumes are calculated by another self-consistency procedure, which guarantees that the dielectric function ε(r) is close to one everywhere inside the protein. The Gaussian widths σ(i) of the atoms i are parameters of the RF approximation. The remarkable accuracy of the method is demonstrated by comparison with Kirkwood's analytical solution for a spherical protein [J. Chem. Phys. 2, 351 (1934)] and with computationally expensive grid-based numerical solutions for simple model systems in dielectric continua including a di-peptide (Ac-Ala-NHMe) as modeled by a standard MM force field. The latter example shows how weakly the RF conformational free energy landscape depends on the parameters σ(i). A summarizing discussion highlights the achievements of the new theory and of its approximate solution particularly by comparison with so-called generalized Born methods. A follow-up paper describes how the method enables Hamiltonian, efficient, and accurate MM molecular dynamics simulations of proteins in dielectric solvent continua.

Testing physical models of passive membrane permeation.

The biophysical basis of passive membrane permeability is well-understood, but most methods for predicting membrane permeability in the context of drug design are based on statistical relationships that indirectly capture the key physical aspects. Here, we investigate molecular mechanics-based models of passive membrane permeability and evaluate their performance against different types of experimental data, including parallel artificial membrane permeability assays (PAMPA), cell-based assays, in vivo measurements, and other in silico predictions. The experimental data sets we use in these tests are diverse, including peptidomimetics, congeneric series, and diverse FDA approved drugs. The physical models are not specifically trained for any of these data sets; rather, input parameters are based on standard molecular mechanics force fields, such as partial charges, and an implicit solvent model. A systematic approach is taken to analyze the contribution from each component in the physics-based permeability model. A primary factor in determining rates of passive membrane permeation is the conformation-dependent free energy of desolvating the molecule, and this measure alone provides good agreement with experimental permeability measurements in many cases. Other factors that improve agreement with experimental data include deionization and estimates of entropy losses of the ligand and the membrane, which lead to size-dependence of the permeation rate.

Case Studies of ONIOM(DFT:DFTB) and ONIOM(DFT:DFTB:MM) for Enzymes and Enzyme Mimics

The replacement of standard molecular mechanics force fields by inexpensive molecular orbital (QM′) methods in multiscale models has many advantages, e.g., a more straightforward description of mutual polarization and charge transfer between layers. The ONIOM(QM:QM′) scheme with mechanical embedding can combine any two methods without prior parametrization or significant coding effort. In this scheme, the environmental effect is evaluated fully at the QM′ level, and the accuracy therefore depends on how well the low-level QM′ method describes the changes in electron density of the reacting region. To examine the applicability of the QM:QM′ approach, we perform case studies with density-functional tight-binding (DFTB) as the low-level QM′ method in two-layer ONIOM(B3LYP/6-31G(d):DFTB) models. The investigated systems include simple amino acid models, one nonheme iron enzyme mimic, and the enzymatic reactions of Zn-β-lactamase and trypsin. For the last example, we also illustrate the use of a three-layer ONIOM(B3LYP/6-31G(d):DFTB:Amber96) model. The ONIOM extension, compared to the QM calculation for the small model system, improves the relative energies, but high accuracy (deviations below 1 kcal/mol) is not achieved even with relatively large QM models. Polarization effects are fairly well described using DFTB, but in some cases QM and QM′ methods converge to different electronic states. We discuss when the QM:QM′ approach is appropriate and the possibilities of estimating the quality of the ONIOM extension without having to make explicit benchmarks of the entire system.

Aten—An application for the creation, editing, and visualization of coordinates for glasses, liquids, crystals, and molecules

Aten is a tool for Linux, Mac, and Windows platforms to ease the creation, editing, and visualization of coordinates for use in, for example, molecular dynamics simulations. The code handles gas-phase molecules in addition to crystals, surfaces, and liquids, providing standard tools to edit "by the atom" along with specific "by the box" methods suited to periodic systems, including full crystallographic spacegroup packing definitions. Visualization of systems encompasses the standard drawing styles and may be mixed with an arbitrary number of other basic objects, such as arrows and geometric objects allowing creation of scenes involving vector fields, coarse-grained particles, etc. Standard molecular mechanics forcefields can be read in, edited, applied to molecules using the built-in chemical typing language, and subsequently used to calculate energies and forces, minimize energies with respect to coordinates, decompose energies into contribution by molecule type, etc. Monte Carlo methods are available with which to generate random configurations of N-component systems, or solvate around molecules or in specific regions within existing configurations. Supported file formats are governed by user-defined "filters," which provide flexibility of input/output as well as the ability to define custom or extended formats of traditional files. Forcefield descriptions for loaded systems can be output and formatted for use as input to common codes. All features can be accessed by a comprehensive GUI and scripting language based on the C syntax.

Accuracy of typical approximations in classical models of intermolecular polarization

One of the largest limitations of standard molecular-mechanics force fields is the neglect of intermolecular polarization. Several attempts to cure this problem have been made, but the results have not always been fully satisfactory. In this paper, we present a quantitative study of the fundamental approximations that underlie polarization models for classical force fields. The induced charge density of a large set of molecular dimers is compared to supermolecular calculations for a hierarchy of simplified models. We study the effect of the Pauli principle, the local inhomogeneity of the electric field, the intramolecular coupling of the polarization response, and the fact that the induced density is a continuous function. We show that standard point-polarizability models work rather well, despite their lack of all these effects, because (1) there is a systematic error cancellation between the neglect of effects of the Pauli principle and the locally inhomogeneous electric field, and (2) the lack of intramolecular coupling and the use of a dipole expansion of the induced density have only minor effects on the polarization. However, the cancellation in (1) is not perfect, and therefore polarizable force-fields could be improved if both effects are explicitly treated.

FlexOracle: predicting flexible hinges by identification of stable domains.

BackgroundProtein motions play an essential role in catalysis and protein-ligand interactions, but are difficult to observe directly. A substantial fraction of protein motions involve hinge bending. For these proteins, the accurate identification of flexible hinges connecting rigid domains would provide significant insight into motion. Programs such as GNM and FIRST have made global flexibility predictions available at low computational cost, but are not designed specifically for finding hinge points.ResultsHere we present the novel FlexOracle hinge prediction approach based on the ideas that energetic interactions are stronger within structural domains than between them, and that fragments generated by cleaving the protein at the hinge site are independently stable. We implement this as a tool within the Database of Macromolecular Motions, MolMovDB.org. For a given structure, we generate pairs of fragments based on scanning all possible cleavage points on the protein chain, compute the energy of the fragments compared with the undivided protein, and predict hinges where this quantity is minimal. We present three specific implementations of this approach. In the first, we consider only pairs of fragments generated by cutting at a single location on the protein chain and then use a standard molecular mechanics force field to calculate the enthalpies of the two fragments. In the second, we generate fragments in the same way but instead compute their free energies using a knowledge based force field. In the third, we generate fragment pairs by cutting at two points on the protein chain and then calculate their free energies.ConclusionQuantitative results demonstrate our method's ability to predict known hinges from the Database of Macromolecular Motions.

Combined valence bond-molecular mechanics potential-energy surface and direct dynamics study of rate constants and kinetic isotope effects for the H+C2H6 reaction

This article presents a multifaceted study of the reaction H+C(2)H(6)-->H(2)+C(2)H(5) and three of its deuterium-substituted isotopologs. First we present high-level electronic structure calculations by the W1, G3SX, MCG3-MPWB, CBS-APNO, and MC-QCISD/3 methods that lead to a best estimate of the barrier height of 11.8+/-0.5 kcal/mol. Then we obtain a specific reaction parameter for the MPW density functional in order that it reproduces the best estimate of the barrier height; this yields the MPW54 functional. The MPW54 functional, as well as the MPW60 functional that was previously parametrized for the H+CH(4) reaction, is used with canonical variational theory with small-curvature tunneling to calculate the rate constants for all four ethane reactions from 200 to 2000 K. The final MPW54 calculations are based on curvilinear-coordinate generalized-normal-mode analysis along the reaction path, and they include scaled frequencies and an anharmonic C-C bond torsion. They agree with experiment within 31% for 467-826 K except for a 38% deviation at 748 K; the results for the isotopologs are predictions since these rate constants have never been measured. The kinetic isotope effects (KIEs) are analyzed to reveal the contributions from subsets of vibrational partition functions and from tunneling, which conspire to yield a nonmonotonic temperature dependence for one of the KIEs. The stationary points and reaction-path potential of the MPW54 potential-energy surface are then used to parametrize a new kind of analytical potential-energy surface that combines a semiempirical valence bond formalism for the reactive part of the molecule with a standard molecular mechanics force field for the rest; this may be considered to be either an extension of molecular mechanics to treat a reactive potential-energy surface or a new kind of combined quantum-mechanical/molecular mechanical (QM/MM) method in which the QM part is semiempirical valence bond theory; that is, the new potential-energy surface is a combined valence bond molecular mechanics (CVBMM) surface. Rate constants calculated with the CVBMM surface agree with the MPW54 rate constants within 12% for 534-2000 K and within 23% for 200-491 K. The full CVBMM potential-energy surface is now available for use in variety of dynamics calculations, and it provides a prototype for developing CVBMM potential-energy surfaces for other reactions.

Analysis and Application of Potential Energy Smoothing and Search Methods for Global Optimization

Global energy optimization of a molecular system is difficult due to the well-known “multiple minimum” problem. The rugged potential energy surface (PES) characteristic of multidimensional systems can be transformed reversibly using potential smoothing to generate a new surface that is easier to search for favorable configurations. The diffusion equation method (DEM) is one example of a potential smoothing algorithm. Potential smoothing as implemented in DEM is intuitively appealing and has certain appropriate statistical mechanical properties, but often fails to identify the global minimum even for relatively small problems. In the present paper, extensions to DEM capable of correcting its empirical behavior are systematically investigated. Two types of local search (LS) procedures are applied during the reversing schedule from the smooth deformed PES to the undeformed surface. Changes needed to generate smoothable versions of standard molecular mechanics force fields such as AMBER/OPLS and MM2 are also described. The resulting methods are applied in an attempt to determine the global energy minimum for a variety of systems in different coordinate representations. The problems studied include argon clusters, cycloheptadecane, capped polyalanine, and the docking of α-helices. Depending on the specific problem, potential smoothing and search (PSS) is performed in Cartesian, torsional, or rigid body space. For example, PSS finds a very low energy structure for cycloheptadecane with much greater efficiency than a search restricted to the undeformed potential surface. It is shown that potential smoothing is characterized by three salient features. As the level of smoothing is increased, unique minima merge into a common basin, crossings can occur in the relative energies of a pair of minima, and the spatial locations of minima are shifted due to the averaging effects of smoothing. Local search procedures improve the ability of smoothing methods to locate global minima because they facilitate the post facto correction of errors due to energy crossings that may have occurred at higher levels of smoothing. PSS methods should serve as useful tools for global energy optimization on a variety of difficult problems of practical interest.

A comparison of the geometry and of the energy results obtained by application of different molecular mechanics force fields to methyl α-lactoside and the C-analogue of lactose

The applicability of a number of standard molecular mechanics force fields, namely AMBER/Homans, CVFF, CFF91, ESFF, and MM3 ∗, and of the AMI semiempirical Hamiltonian, to the molecular modelling of methyl α-lactoside is described. A protocol based on the calculation of relaxed energy maps for a significant number of initial structures has been employed. Different dielectric constants, as well as a continuum model solvent have been used. The geometric characteristics (distances and torsion and valence angles) of the low energy minima have been gathered and compared to experimental solution NMR and solid state X-ray data. Concerning the energies, the relative steric energy values provided by the different programs have been translated into ensemble average distributions of conformers which, in turn, have been employed to generate NMR parameters, i.e. NOE, NOESY, and ROESY intensities and homo- and hetero-nuclear coupling constants. The results indicate that MM3 ∗ and CVFF provide conformer distributions which account for the solution NMR observables in a satisfactory way. AMBER/Homans provide a reasonable agreement, while conformational analysis of oligosaccharides by AM1, ESFF, and CFF91 is not advised. The best fit between the experimental X-ray distances and angles and the modeling studies is found when CVFF is used. Results from the different force fields (and AM1) for C-lactose were compared.

Interpretation of AFM images: the graphite surface with a diamond tip

Constant force AFM images of the graphite surface using a diamond tip are calculated from a standard molecular mechanics force field. It is shown that there is no need for an asymmetric tip apex nor for a graphite flake under the apex to reproduce the experimental images. A good orientation of the tip apex and a good choice of the force applied is enough to evolute from a trigonal to a honeycomb AFM apparent surface structure for graphite.

Peptide conformational analysis using the TRIPOS force field

AbstractTo explain observed electron transfer rates between the terminal aromatic moieties in the arylsubstituted alanine dimer 1, a conformational analysis of the dipeptide was performed using a standard molecular mechanics force field. The resulting low‐energy conformers could conveniently be grouped into two families with an average energy difference of ca. 2 kcal mol−1 and populations of 64% and 30%, respectively. These results correlate well with the 78:22 ratio of the two experimentally distinguishable decay processes for the radical anion of 1. © 1992 John Wiley & Sons, Inc.

Computational method for the design of enzymes with altered substrate specificity

A combination of enzyme kinetics and X-ray crystallographic analysis of site-specific mutants has been used to probe the determinants of substrate specificity for the enzyme α-lytic protease. We now present a generalized model for understanding the effects of mutagenesis on enzyme substrate specificity. This algorithm uses a library of side-chain rotamers to sample conformation space within the binding site for the enzyme—substrate complex. The free energy of each conformation is evaluated with a standard molecular mechanics force field, modified to include a solvation energy term. This rapid energy calculation based on coarse conformation sampling quite accurately predicts the relative catalytic efficiency of over 40 different α-lytic protease—substrate combinations. Unlike other computational approaches, with this method it is feasible to evaluate all possible mutations within the binding site. Using this algorithm, we have successfully designed a protease that is both highly active and selective for a non-natural substrate. These encouraging results indicate that it is possible to design altered enzymes solely on the basis of empirical energy calculations.