Accurate interaction energies of hydrogen-bonded nucleic acid base pairs

Hydrogen-bonded nucleic acids base pairs substantially contribute to the structure and stability of nucleic acids. The study presents reference ab initio structures and interaction energies of selected base pairs with binding energies ranging from -5 to -47 kcal/mol. The molecular structures are obtained using the RI-MP2 (resolution of identity MP2) method with extended cc-pVTZ basis set of atomic orbitals. The RI-MP2 method provides results essentially identical with the standard MP2 method. The interaction energies are calculated using the Complete Basis Set (CBS) extrapolation at the RI-MP2 level. For some base pairs, Coupled-Cluster corrections with inclusion of noniterative triple contributions (CCSD(T)) are given. The calculations are compared with selected medium quality methods. The PW91 DFT functional with the 6-31G basis set matches well the RI-MP2/CBS absolute interaction energies and reproduces the relative values of base pairing energies with a maximum relative error of 2.6 kcal/mol when applied with Becke3LYP-optimized geometries. The Becke3LYP DFT functional underestimates the interaction energies by few kcal/mol with relative error of 2.2 kcal/mol. Very good performance of nonpolarizable Cornell et al. force field is confirmed and this indirectly supports the view that H-bonded base pairs are primarily stabilized by electrostatic interactions.     

Continuum polarizable force field within the Poisson-Boltzmann framework

We have developed and tested a complete set of nonbonded parameters for a continuum polarizable force field. Our analysis shows that the new continuum polarizable model is consistent with B3LYP/cc-pVTZ in modeling electronic response upon variation of dielectric environment. Comparison with experiment also shows that the new continuum polarizable model is reasonable, with accuracy similar to that of B3LYP/cc-pVTZ in reproduction of dipole moments of selected organic molecules in the gas phase. We have further tested the validity to interchange the Amber van der Waals parameters between the explicit and continuum polarizable force fields with a series of dimers. It can be found that the continuum polarizable model agrees well with MP2/cc-pVTZ, with deviations in dimer binding energies less than 0.9 kcal/mol in the aqueous dielectric environment. Finally, we have optimized atomic cavity radii with respect to experimental solvation free energies of 177 training molecules. To validate the optimized cavity radii, we have tested these parameters against 176 test molecules. It is found that the optimized Poisson-Boltzmann atomic cavity radii transfer well from the training set to the test set, with an overall root-mean-square deviation of 1.30 kcal/mol, an unsigned average error of 1.07 kcal/mol, and a correlation coefficient of 92% for all 353 molecules in both the training and test sets. Given the development documented here, the next natural step is the construction of a full protein/nucleic acid force field within the new continuum polarization framework.     

Predicting hydration free energies using all-atom molecular dynamics simulations and multiple starting conformations

Molecular dynamics simulations in explicit solvent were applied to predict the hydration free energies for 23 small organic molecules in blind SAMPL2 test. We found good agreement with experimental results, with an RMS error of 2.82 kcal/mol over the whole set and 1.86 kcal/mol over all the molecules except several hydroxyl-rich compounds where we find evidence for a systematic error in the force field. We tested two different solvent models, TIP3P and TIP4P-Ew, and obtained very similar hydration free energies for these two models; the RMS difference was 0.64 kcal/mol. We found that preferred conformation of the carboxylic acids in water differs from that in vacuum. Surprisingly, this conformational change is not adequately sampled on simulation timescales, so we apply an umbrella sampling technique to include free energies associated with the conformational change. Overall, the results of this test reveal that the force field parameters for some groups of molecules (such as hydroxyl-rich compounds) still need to be improved, but for most compounds, accuracy was consistent with that seen in our previous tests.     

Surveying implicit solvent models for estimating small molecule absolute hydration free energies

Implicit solvent models are powerful tools in accounting for the aqueous environment at a fraction of the computational expense of explicit solvent representations. Here, we compare the ability of common implicit solvent models (TC, OBC, OBC2, GBMV, GBMV2, GBSW, GBSW/MS, GBSW/MS2 and FACTS) to reproduce experimental absolute hydration free energies for a series of 499 small neutral molecules that are modeled using AMBER/GAFF parameters and AM1-BCC charges. Given optimized surface tension coefficients for scaling the surface area term in the nonpolar contribution, most implicit solvent models demonstrate reasonable agreement with extensive explicit solvent simulations (average difference 1.0-1.7 kcal/mol and R(2)=0.81-0.91) and with experimental hydration free energies (average unsigned errors=1.1-1.4 kcal/mol and R(2)=0.66-0.81). Chemical classes of compounds are identified that need further optimization of their ligand force field parameters and others that require improvement in the physical parameters of the implicit solvent models themselves. More sophisticated nonpolar models are also likely necessary to more effectively represent the underlying physics of solvation and take the quality of hydration free energies estimated from implicit solvent models to the next level.     

Combining the polarizable Drude force field with a continuum electrostatic Poisson–Boltzmann implicit solvation model

In this work, we have combined the polarizable force field based on the classical Drude oscillator with a continuum Poisson–Boltzmann/solvent‐accessible surface area (PB/SASA) model. In practice, the positions of the Drude particles experiencing the solvent reaction field arising from the fixed charges and induced polarization of the solute must be optimized in a self‐consistent manner. Here, we parameterized the model to reproduce experimental solvation free energies of a set of small molecules. The model reproduces well‐experimental solvation free energies of 70 molecules, yielding a root mean square difference of 0.8 kcal/mol versus 2.5 kcal/mol for the CHARMM36 additive force field. The polarization work associated with the solute transfer from the gas‐phase to the polar solvent, a term neglected in the framework of additive force fields, was found to make a large contribution to the total solvation free energy, comparable to the polar solute–solvent solvation contribution. The Drude PB/SASA also reproduces well the electronic polarization from the explicit solvent simulations of a small protein, BPTI. Model validation was based on comparisons with the experimental relative binding free energies of 371 single alanine mutations. With the Drude PB/SASA model the root mean square deviation between the predicted and experimental relative binding free energies is 3.35 kcal/mol, lower than 5.11 kcal/mol computed with the CHARMM36 additive force field. Overall, the results indicate that the main limitation of the Drude PB/SASA model is the inability of the SASA term to accurately capture non‐polar solvation effects. © 2018 Wiley Periodicals, Inc.     

Crystal structure prediction of flexible molecules using parallel genetic algorithms with a standard force field

This article describes the application of our distributed computing framework for crystal structure prediction (CSP) the modified genetic algorithms for crystal and cluster prediction (MGAC), to predict the crystal structure of flexible molecules using the general Amber force field (GAFF) and the CHARMM program. The MGAC distributed computing framework includes a series of tightly integrated computer programs for generating the molecule's force field, sampling crystal structures using a distributed parallel genetic algorithm and local energy minimization of the structures followed by the classifying, sorting, and archiving of the most relevant structures. Our results indicate that the method can consistently find the experimentally known crystal structures of flexible molecules, but the number of missing structures and poor ranking observed in some crystals show the need for further improvement of the potential. © 2009 Wiley Periodicals, Inc. J Comput Chem, 2009     

The tetracycline: Mg2+ complex: a molecular mechanics force field

Tetracycline (Tc) is an important antibiotic, which binds specifically to the ribosome and several proteins, in the form of a Tc-:Mg2+ complex. To model Tc:protein and Tc:RNA interactions, we have developed a molecular mechanics force field model of Tc, which is consistent with the CHARMM force field for proteins and nucleic acids. We used structures from the Cambridge Crystallographic Data Base to identify the main Tc conformations that are likely to be present in solution and in biomolecular complexes. A conformational search was also done, using the MM3 force field to perform simulated annealing of Tc. Several resulting, low-energy structures were optimized with an ab initio model and used in developing the new Tc force field. Atomic charges and Lennard-Jones parameters were derived from a supermolecule ab initio approach. We considered the ab initio energies and geometries of a probe water molecule interacting with Tc at 36 different positions. We considered both a neutral and a zwitterionic Tc form, with and without bound Mg2+. The final rms deviation between the ab initio and force field energies, averaged over all forms, was just 0.35 kcal/mol. The model also reproduces the ab initio geometry and flexibility of Tc. As further tests, we did simulations of a Tc crystal, of Tc:Mg2+ and Tc:Ca2+ complexes in aqueous solution, and of a solvated complex between Tc:Mg2+ and the Tet repressor protein (TetR). With slight, ad hoc adjustments, the model can reproduce the experimental, relative, Tc binding affinities of Mg2+ and Ca2+. It performs well for the structure and fluctuations of the Tc:Mg2+:TetR complex. The model should therefore be suitable to investigate the interactions of Tc with proteins and RNA. It provides a starting point to parameterize other compounds in the large Tc family.     

OPLS all-atom force field for carbohydrates

The OPLS all-atom (AA) force field for organic and biomolecular systems has been expanded to include carbohydrates. Starting with reported nonbonded parameters of alcohols, ethers, and diols, torsional parameters were fit to reproduce results from ab initio calculations on the hexopyranoses, α,β-d -glucopyranose, α,β-d -mannopyranose, α,β-d -galactopyranose, methyl α,β-d -glucopyranoside, and methyl α,β-d -mannopyranoside. In all, geometry optimizations were carried out for 144 conformers at the restricted Hartree–Fock (RHF)/6–31G* level. For the conformers with a relative energy within 3 kcal/mol of the global minima, the effects of electron correlation and basis-set extension were considered by performing single-point calculations with density functional theory at the B3LYP/6–311+G** level. The torsional parameters for the OPLS-AA force field were parameterized to reproduce the energies and structures of these 44 conformers. The resultant force field reproduces the ab initio calculated energies with an average unsigned error of 0.41 kcal/mol. The α/β ratios as well as the relative energies between the isomeric hexopyranoses are in good accord with the ab initio results. The predictive abilities of the force field were also tested against RHF/6–31G* results for d -allopyranose with excellent success; a surprising discovery is that the lowest energy conformer of d -allopyranose is a β anomer. © 1997 John Wiley & Sons, Inc. J Comput Chem 18 : 1955–1970, 1997     

Evaluation of solvation free energies for small molecules with the AMOEBA polarizable force field

The effects of electronic polarization in biomolecular interactions will differ depending on the local dielectric constant of the environment, such as in solvent, DNA, proteins, and membranes. Here the performance of the AMOEBA polarizable force field is evaluated under nonaqueous conditions by calculating the solvation free energies of small molecules in four common organic solvents. Results are compared with experimental data and equivalent simulations performed with the GAFF pairwise‐additive force field. Although AMOEBA results give mean errors close to “chemical accuracy,” GAFF performs surprisingly well, with statistically significantly more accurate results than AMOEBA in some solvents. However, for both models, free energies calculated in chloroform show worst agreement to experiment and individual solutes are consistently poor performers, suggesting non‐potential‐specific errors also contribute to inaccuracy. Scope for the improvement of both potentials remains limited by the lack of high quality experimental data across multiple solvents, particularly those of high dielectric constant. © 2016 The Authors. Journal of Computational Chemistry Published by Wiley Periodicals, Inc.     

Fast,efficient generation of high-quality atomic charges. AM1-BCC model: II. Parameterization and validation

We present the first global parameterization and validation of a novel charge model, called AM1-BCC, which quickly and efficiently generates high-quality atomic charges for computer simulations of organic molecules in polar media. The goal of the charge model is to produce atomic charges that emulate the HF/6-31G* electrostatic potential (ESP) of a molecule. Underlying electronic structure features, including formal charge and electron delocalization, are first captured by AM1 population charges; simple additive bond charge corrections (BCCs) are then applied to these AM1 atomic charges to produce the AM1-BCC charges. The parameterization of BCCs was carried out by fitting to the HF/6-31G* ESP of a training set of >2700 molecules. Most organic functional groups and their combinations were sampled, as well as an extensive variety of cyclic and fused bicyclic heteroaryl systems. The resulting BCC parameters allow the AM1-BCC charging scheme to handle virtually all types of organic compounds listed in The Merck Index and the NCI Database. Validation of the model was done through comparisons of hydrogen-bonded dimer energies and relative free energies of solvation using AM1-BCC charges in conjunction with the 1994 Cornell et al. forcefield for AMBER.(13) Homo- and hetero-dimer hydrogen-bond energies of a diverse set of organic molecules were reproduced to within 0.95 kcal/mol RMS deviation from the ab initio values, and for DNA dimers the energies were within 0.9 kcal/mol RMS deviation from ab initio values. The calculated relative free energies of solvation for a diverse set of monofunctional isosteres were reproduced to within 0.69 kcal/mol of experiment. In all these validation tests, AMBER with the AM1-BCC charge model maintained a correlation coefficient above 0.96. Thus, the parameters presented here for use with the AM1-BCC method present a fast, accurate, and robust alternative to HF/6-31G* ESP-fit charges for general use with the AMBER force field in computer simulations involving organic small molecules.     

Polarizable Atomic Multipole-based Molecular Mechanics for Organic Molecules

An empirical potential based on permanent atomic multipoles and atomic induced dipoles is reported for alkanes, alcohols, amines, sulfides, aldehydes, carboxylic acids, amides, aromatics and other small organic molecules. Permanent atomic multipole moments through quadrupole moments have been derived from gas phase ab initio molecular orbital calculations. The van der Waals parameters are obtained by fitting to gas phase homodimer QM energies and structures, as well as experimental densities and heats of vaporization of neat liquids. As a validation, the hydrogen bonding energies and structures of gas phase heterodimers with water are evaluated using the resulting potential. For 32 homo- and heterodimers, the association energy agrees with ab initio results to within 0.4 kcal/mol. The RMS deviation of hydrogen bond distance from QM optimized geometry is less than 0.06 ?. In addition, liquid self-diffusion and static dielectric constants computed from molecular dynamics simulation are consistent with experimental values. The force field is also used to compute the solvation free energy of 27 compounds not included in the parameterization process, with a RMS error of 0.69 kcal/mol. The results obtained in this study suggest the AMOEBA force field performs well across different environments and phases. The key algorithms involved in the electrostatic model and a protocol for developing parameters are detailed to facilitate extension to additional molecular systems.     

Alcohols,ethers, carbohydrates,and related compounds. IV. Carbohydrates

Ab initio calculations [B3LYP/6-311++G(2d,2p)] have been carried out on 84 conformations of 12 different sugars (hexoses), in both pyranose and furanose forms, with the idea of generating a data base for carbohydrate structural energies that may be used for developing the predictive value of molecular mechanics calculations for carbohydrates. The average value for the apparent gas phase anomeric effect for a series of 31 pairs of pyranose conformations was found to be 1.83 kcal/mol (vs. 2.67 kcal/mol with a smaller basis set used in earlier calculations). In developing MM4 to reproduce these data, it was necessary first to have good energies for simple alcohols and ethers, together with an adequate treatment of hydrogen bonding, and then to include the anomeric effect, and the ethylene glycol type system, as was previously recognized. It was also found that the so-called delta-2 effect, long recognized in carbohydrates, must be explicitly included, in order to obtain acceptable results. When a force field that included all of these items as developed from the small molecules based on the MM4 hydrocarbon force field was applied without any parameter adjustment to the set of hexopyranose and furanose conformations mentioned earlier, the E(beta) - E(alpha) was found to have an average value of 1.88 kcal/mol, versus 1.74 for the quantum calculations. The signed average and RMS deviations of the MM4 from the QM results were +0.15 and 0.87 kcal/mol.     

Polarizable model potential function for nucleic acid bases

A polarizable model potential (PMP) function for adenine (A), cytosine (C), guanine (G), thymine (T), and uracil (U) is developed on the basis of ab initio molecular orbital calculations at the MP2/6-31+G* level. The PMP function consists of Coulomb, van der Waals, and polarization terms. The permanent atomic charges of the Coulomb term are determined by using electrostatic potential (ESP) optimization. The multicenter polarizabilities of the polarization term are determined by using polarized one-electron potential (POP) optimization in which the electron density changes induced by a test charge are target. Isotropic and anisotropic polarizabilities are adopted as the multicenter polarizabilities. In the PMP calculations using the optimized parameters, the interaction energies of Watson-Crick type A-T and C-G base pairs were -15.6 and -29.4 kcal/mol, respectively. The interaction energy of Hoogsteen type A-T base pair was -17.8 kcal/mol. These results reproduce well the quantum chemistry calculations at the MP2/6-311++G(3df,2pd) level within the differences of 0.6 kcal/mol. The stacking energies of A-T and C-G were -9.7 and -10.9 kcal/mol. These reproduce well the calculation results at the MP2/6-311++G (2d,2p) level within the differences of 1.3 kcal/mol. The potential energy surfaces of the system in which a sodium ion or a chloride ion is adjacent to the nucleic acid base are calculated. The interaction energies of the PMP function reproduced well the calculation results at the MP2/6-31+G* or MP2/6-311++G(2d,2p) level. The reason why the PMP function reproduces well the high-level quantum mechanical interaction energies is addressed from the viewpoint of each energy terms.     

Blind prediction of solvation free energies from the SAMPL4 challenge

Here, we give an overview of the small molecule hydration portion of the SAMPL4 challenge, which focused on predicting hydration free energies for a series of 47 small molecules. These gas-to-water transfer free energies have in the past proven a valuable test of a variety of computational methods and force fields. Here, in contrast to some previous SAMPL challenges, we find a relatively wide range of methods perform quite well on this test set, with RMS errors in the 1.2 kcal/mol range for several of the best performing methods. Top-performers included a quantum mechanical approach with continuum solvent models and functional group corrections, alchemical molecular dynamics simulations with a classical all-atom force field, and a single-conformation Poisson–Boltzmann approach. While 1.2 kcal/mol is still a significant error, experimental hydration free energies covered a range of nearly 20 kcal/mol, so methods typically showed substantial predictive power. Here, a substantial new focus was on evaluation of error estimates, as predicting when a computational prediction is reliable versus unreliable has considerable practical value. We found, however, that in many cases errors are substantially underestimated, and that typically little effort has been invested in estimating likely error. We believe this is an important area for further research.     

Internal coordinate modeling of DNA: Force field comparisons

The Jumna internal coordinate program for modeling nucleic acids was extended to include the force field developed for the Amber program. This forms a bridge between internal and Cartesian coordinate modeling approaches. Using the extensive conformational mapping and substate search facilities available within Jumna, we rigorously compared the behavior of the different force fields and also of different continuum solvent models. The results, which help to explain trends seen in earlier minimization and molecular dynamics simulations, point to the superiority of the latest Amber parameterization (Parm94) and to a surprising degree of agreement with the Flex force field originally developed for Jumna. © 1997 John Wiley & Sons, Inc. J Comput Chem 18 :1043–1055, 1997     

A sobering assessment of small‐molecule force field methods for low energy conformer predictions

We have carried out a large scale computational investigation to assess the utility of common small‐molecule force fields for computational screening of low energy conformers of typical organic molecules. Using statistical analyses on the energies and relative rankings of up to 250 diverse conformers of 700 different molecular structures, we find that energies from widely used classical force fields (MMFF94, UFF, and GAFF) show unconditionally poor energy and rank correlation with semiempirical (PM7) and Kohn–Sham density functional theory (DFT) energies calculated at PM7 and DFT optimized geometries. In contrast, semiempirical PM7 calculations show significantly better correlation with DFT calculations and generally better geometries. With these results, we make recommendations to more reliably carry out conformer screening.     

Study of peptide conformation in terms of the ABEEM/MM method

The ABEEM/MM model (atom-bond electronegativity equalization method fused into molecular mechanics) is applied to study of the polypeptide conformations. The Lennard-Jones and torsional parameters were optimized to be consistent with the ABEEM/MM fluctuating charge electrostatic potential. The hydrogen bond was specially treated with an electrostatic fitting function. Molecular dipole moments, dimerization energies, and hydrogen bond lengths of complexes are reasonably achieved by our model, compared to ab initio results. The ABEEM/MM fluctuating charge model reproduces both the peptide conformational energies and structures with satisfactory accuracy with low computer cost. The transferability is tested by applying the parameters of our model to the tetrapeptide of alanine and another four dipeptides. The overall RMS deviations in conformational energies and key dihedral angles for four di- or tetrapeptide, is 0.39 kcal/mol and 7.7 degrees . The current results agree well with those by the accurate ab initio method, and are comparable to those from the best existing force fields. The results make us believe that our fluctuating charge model can obtain more promising results in protein and macromolecular modeling with good accuracy but less computer cost.     

Parameterization of OPLS-AA force field for the conformational analysis of macrocyclic polyketides

The parameters for the OPLS-AA potential energy function have been extended to include some functional groups that are present in macrocyclic polyketides. Existing OPLS-AA torsional parameters for alkanes, alcohols, ethers, hemiacetals, esters, and ketoamides were improved based on MP2/aug-cc-pVTZ and MP2/aug-cc-pVDZ calculations. Nonbonded parameters for the sp(3) carbon and oxygen atoms were refined using Monte Carlo simulations of bulk liquids. The resulting force field predicts conformer energies and torsional barriers of alkanes, alcohols, ethers, and hemiacetals with an overall RMS deviation of 0.40 kcal/mol as compared to reference data. Densities of 19 bulk liquids are predicted with an average error of 1.1%, and heats of vaporization are reproduced within 2.4% of experimental values. The force field was used to perform conformational analysis of smaller analogs of the macrocyclic polyketide drug FK506. Structures that adopted low-energy conformations similar to that of bound FK506 were identified. The results show that a linker of four ketide units constitutes the shortest effector domain that allows binding of the ketide drugs to FKBP proteins. It is proposed that the exact chemical makeup of the effector domain has little influence on the conformational preference of tetraketides.     

Ab initio force field for simulations of proteins and nucleic acids

We present a new self-consistent set of ab initio analytical pair potential to predict specific nonbonded interactions of protein with nucleic acid, of protein with protein, and of nucleic acid with nucleic acid. The purpose of this study is to represent the interaction between biological molecules with an accuracy equivalent to the ab initio molecular orbital calculations, which are used as reference data to obtain the pair potentials. Atoms in nucleic acids and proteins are classified according to their chemical environments. An “effective charge,” a modification of a charge obtained from the Mulliken population analysis, is introduced and used to represent the electrostatic energy. More than 30,000 SCF interaction energies have been calculated to provide the reference data for the fitting procedure that we have adopted in the parameterization of the potentials. The standard deviation is 1.61 kcal/mol for interaction energies spanning the range from about ?220 kcal/mol to +20 kcal/mol. Molecular dynamics simulations, using the above new set of force field, have been performed successfully for the systems where adequate treatments of specific interactions are required: The stability of α-helix of C-peptide and the interaction of spermine with oligonucleotide are examined as preliminary examples.     

Hydrogen-bonded nucleic acid base pairs containing unusual base tautomers: complete basis set calculations at the MP2 and CCSD(T) levels

The total interaction energies of altogether 15 hydrogen-bonded nucleic acid base pairs containing unusual base tautomers were calculated. The geometry properties of all selected adenine-thymine and guanine-cytosine hydrogen-bonded base pairs enable their incorporation into DNA. Unusual base pairing patterns were compared with Watson-Crick H-bonded structures of the adenine-thymine and guanine-cytosine pairs. The complete basis set (CBS) limit of the MP2 interaction energy and the CCSD(T) correction term, determined as the difference between the CCSD(T) and MP2 interaction energies, was evaluated. Extrapolation to the MP2 CBS limit was done using the aug-cc-pVDZ and aug-cc-pVTZ results, and the CCSD(T) correction term was determined with the 6-31G*(0.25) basis set. Final interaction energies were corrected while taking into account both tautomeric penalization determined at the CBS level and solvation/desolvation free energies. The situation for the adenine-thymine pairs is straightforward, and tautomeric pairs are significantly less stable than the Watson-Crick pair consisting of the canonical forms. In the case of the guanine-cytosine pair, the Watson-Crick structure made by canonical forms is again the most stable. The other two structures are, however, energetically rather similar (by 5 and 6 kcal/mol), which provides a very small but non-negligible chance of detecting these structures in the DNA double helix (1:5000). Due to the fact that DNA bases and base pairs incorporated into DNA are solvated less favorably than in isolated systems, this probability represents the very upper limit. The results clearly show how precisely the canonical building blocks of DNA molecules were chosen and how well their stability is maintained.