Evaluation of DNA Force Fields in Implicit Solvation

DNA structural deformations and dynamics are crucial to its interactions in the cell. Theoretical simulations are essential tools to explore the structure, dynamics, and thermodynamics of biomolecules in a systematic way. Molecular mechanics force fields for DNA have benefited from constant improvements during the last decades. Several studies have evaluated and compared available force fields when the solvent is modeled by explicit molecules. On the other hand, few systematic studies have assessed the quality of duplex DNA models when implicit solvation is employed. The interest of an implicit modeling of the solvent consists in the important gain in the simulation performance and conformational sampling speed. In this study, respective influences of the force field and the implicit solvation model choice on DNA simulation quality are evaluated. To this end, extensive implicit solvent duplex DNA simulations are performed, attempting to reach both conformational and sequence diversity convergence. Structural parameters are extracted from simulations and statistically compared to available experimental and explicit solvation simulation data. Our results quantitatively expose the respective strengths and weaknesses of the different DNA force fields and implicit solvation models studied. This work can lead to the suggestion of improvements to current DNA theoretical models.     

Assessing the performance of the g_mmpbsa tools to simulate the inhibition of oseltamivir to influenza virus neuraminidase by molecular mechanics Poisson–Boltzmann surface area methods

The molecular mechanics Poisson–Boltzmann surface area (MM-PBSA) method for GROMACS (g_mmpbsa) is an open-source tool that is capable of reading the trajectories generated by GROMACS and calculating the binding free energy using the MM-PBSA method. However, there are multiple force fields available for users to choose from in the GROMACS software, and there are also different solvent water models to combine with the chosen force fields. These different combinations of parameters may significantly impact the results of g_mmpbsa calculation. Unfortunately, the exact combination of force field and solvent water that can well calculate the free energy of the receptor–ligand binding in GROMACS has not been explored yet. To resolve the above issues, this study mainly explored the molecular dynamics (MD) simulations by GROMACS with the six commonly used force fields and three solvent water models, in combination with g_mmpbsa, to calculate the binding free energies of the influenza virus neuraminidase and its mutants with inhibitor oseltamivir carboxylate and compared the present results with previous published results of Amber software from ours and other researchers. Finally, we provided an optimized calculation model, as well as suggestions that may serve as advice and guidance for future computer-aided designs of drug molecules.     

力场/溶剂水模型对EGFRvⅢ抗原-MR1(scFv)抗体复合物的MM-PBSA计算精度的影响与实验验证

以表皮生长因子Ⅲ型突变体(EGFRvⅢ)抗原多肽与其抗体(MR1)及其人源化突变体的复合物结构为出发点,采用分子动力学中的6种常用力场及3种常用溶剂水模型,分别对上述抗原-抗体复合物进行100ns的分子动力学模拟与分子力学和连续介质模型计算自由能(MM-PBSA),并在实验上利用等温滴定量热(ITC)仪测定了抗原和抗体相互作用的热力学参数.通过在结构变化、能量变化及野生型与突变体比较等几个方面进行综合分析,给出了最佳的计算模型.对不同力场及水模型计算精度等相关问题进行了探讨.     

MATCH: an atom-typing toolset for molecular mechanics force fields

We introduce a toolset of program libraries collectively titled multipurpose atom-typer for CHARMM (MATCH) for the automated assignment of atom types and force field parameters for molecular mechanics simulation of organic molecules. The toolset includes utilities for the conversion of multiple chemical structure file formats into a molecular graph. A general chemical pattern-matching engine using this graph has been implemented whereby assignment of molecular mechanics atom types, charges, and force field parameters are achieved by comparison against a customizable list of chemical fragments. While initially designed to complement the CHARMM simulation package and force fields by generating the necessary input topology and atom-type data files, MATCH can be expanded to any force field and program, and has core functionality that makes it extendable to other applications such as fragment-based property prediction. In this work, we demonstrate the accurate construction of atomic parameters of molecules within each force field included in CHARMM36 through exhaustive cross validation studies illustrating that bond charge increment rules derived from one force field can be transferred to another. In addition, using leave-one-out substitution it is shown that it is also possible to substitute missing intra and intermolecular parameters with ones included in a force field to complete the parameterization of novel molecules. Finally, to demonstrate the robustness of MATCH and the coverage of chemical space offered by the recent CHARMM general force field (Vanommeslaeghe, et al., J Comput Chem 2010, 31, 671), one million molecules from the PubChem database of small molecules are typed, parameterized, and minimized.     

Mixed ab initio QM/MM modeling using frozen orbitals and tests with alanine dipeptide and tetrapeptide

Methodology is discussed for mixed ab initio quantum mechanics/molecular mechanics modeling of systems where the quantum mechanics (QM) and molecular mechanics (MM) regions are within the same molecule. The ab initio QM calculations are at the restricted Hartree–Fock level using the pseudospectral method of the Jaguar program while the MM part is treated with the OPLS force fields implemented in the IMPACT program. The interface between the QM and MM regions, in particular, is elaborated upon, as it is dealt with by “breaking” bonds at the boundaries and using Boys-localized orbitals found from model molecules in place of the bonds. These orbitals are kept frozen during QM calculations. Results from tests of the method to find relative conformational energies and geometries of alanine dipeptides and alanine tetrapeptides are presented along with comparisons to pure QM and pure MM calculations. ©1999 John Wiley & Sons, Inc. J Comput Chem 20: 1468–1494, 1999     

AGBNP: an analytic implicit solvent model suitable for molecular dynamics simulations and high-resolution modeling

We have developed an implicit solvent effective potential (AGBNP) that is suitable for molecular dynamics simulations and high-resolution modeling. It is based on a novel implementation of the pairwise descreening Generalized Born model for the electrostatic component and a new nonpolar hydration free energy estimator. The nonpolar term consists of an estimator for the solute-solvent van der Waals dispersion energy designed to mimic the continuum solvent solute-solvent van der Waals interaction energy, in addition to a surface area term corresponding to the work of cavity formation. AGBNP makes use of a new parameter-free algorithm to calculate the scaling coefficients used in the pairwise descreening scheme to take into account atomic overlaps. The same algorithm is also used to calculate atomic surface areas. We show that excellent agreement is achieved for the GB self-energies and surface areas in comparison to accurate, but much more expensive, numerical evaluations. The parameter-free approach used in AGBNP and the sensitivity of the AGBNP model with respect to large and small conformational changes makes the model suitable for high-resolution modeling of protein loops and receptor sites as well as high-resolution prediction of the structure and thermodynamics of protein-ligand complexes. We present illustrative results for these kinds of benchmarks. The model is fully analytical with first derivatives and is computationally efficient. It has been incorporated into the IMPACT molecular simulation program.     

Helix formation in a pentapeptide: experiment and force-field dependent dynamics

We used a combined approach of experiment and simulation to determine the helical population and folding pathway of a small helix forming blocked pentapeptide, Ac-(Ala)(5)-NH(2). Experimental structural characterization of this blocked peptide was carried out with far UV circular dichroism spectroscopy, FTIR, and NMR measurements. These measurements confirm the presence of the α-helical state in a buffer solution. Direct molecular dynamics and replica-exchange simulations of the pentapeptide were performed using several popular force fields with explicit solvent. The simulations yielded statistically reliable estimates of helix populations, melting curves, folding, and nucleation times. The distributions of conformer populations are used to measure folding cooperativity. Finally, a statistical analysis of the sample of helix-coil transition paths was performed. The details of the calculated helix populations, folding kinetics and pathways vary with the employed force field. Interestingly, the helix populations, folding, and unfolding times obtained from most of the studied force fields are in qualitative agreement with each other and with available experimental data, with the deviations corresponding to several kcal/mol in energy at 300 K. Most of the force fields also predict qualitatively similar transition paths, with unfolding initiated at the C-terminus. Accuracy of potential energy parameters, rather than conformational sampling may be the limiting factor in current molecular simulations.     

Analyzing the robustness of the MM/PBSA free energy calculation method: application to DNA conformational transitions

The ability to predict and characterize free energy differences associated with conformational equilibria or the binding of biomolecules is vital to understanding the molecular basis of many important biological functions. As biological studies focus on larger molecular complexes and properties of the genome, proteome, and interactome, the development and characterization of efficient methods for calculating free energy becomes increasingly essential. The aim of this study is to examine the robustness of the end-point free energy method termed the molecular mechanics Poisson-Boltzmann solvent accessible surface area (MM/PBSA) method. Specifically, applications of MM/PBSA to the conformational equilibria of nucleic acid (NA) systems are explored. This is achieved by comparing A to B form DNA conformational free energy differences calculated using MM/PBSA with corresponding free energy differences determined with a more rigorous and time-consuming umbrella sampling algorithm. In addition, the robustness of NA MM/PBSA calculations is also evaluated in terms of the sensitivity towards the choice of force field and the choice of solvent model used during conformational sampling. MM/PBSA calculations of the free energy difference between A-form and B-form DNA are shown to be in very close agreement with the PMF result determined using an umbrella sampling approach. Further, it is found that the MM/PBSA conformational free energy differences were also in agreement using either the CHARMM or AMBER force field. The influence of ionic strength on conformational stability was particularly insensitive to the choice of force field. Finally, it is also shown that the use of a generalized Born implicit solvent during conformational sampling results in free energy estimates that deviate slightly from those obtained using explicitly solvated MD simulations in these NA systems.     

Comparison of Secondary Structure Formation Using 10 Different Force Fields in Microsecond Molecular Dynamics Simulations

We have compared molecular dynamics (MD) simulations of a β-hairpin forming peptide derived from the protein Nrf2 with 10 biomolecular force fields using trajectories of at least 1 μs. The total simulation time was 37.2 μs. Previous studies have shown that different force fields, water models, simulation methods, and parameters can affect simulation outcomes. The MD simulations were done in explicit solvent with a 16-mer Nrf2 β-hairpin forming peptide using Amber ff99SB-ILDN, Amber ff99SB*-ILDN, Amber ff99SB, Amber ff99SB*, Amber ff03, Amber ff03*, GROMOS96 43a1p, GROMOS96 53a6, CHARMM27, and OPLS-AA/L force fields. The effects of charge-groups, terminal capping, and phosphorylation on the peptide folding were also examined. Despite using identical starting structures and simulation parameters, we observed clear differences among the various force fields and even between replicates using the same force field. Our simulations show that the uncapped peptide folds into a native-like β-hairpin structure at 310 K when Amber ff99SB-ILDN, Amber ff99SB*-ILDN, Amber ff99SB, Amber ff99SB*, Amber ff03, Amber ff03*, GROMOS96 43a1p, or GROMOS96 53a6 were used. The CHARMM27 simulations were able to form native hairpins in some of the elevated temperature simulations, while the OPLS-AA/L simulations did not yield native hairpin structures at any temperatures tested. Simulations that used charge-groups or peptide capping groups were not largely different from their uncapped counterparts with single atom charge-groups. On the other hand, phosphorylation of the threonine residue located at the β-turn significantly affected the hairpin formation. To our knowledge, this is the first study comparing such a large set of force fields with respect to β-hairpin folding. Such a comprehensive comparison will offer useful guidance to others conducting similar types of simulations.     

乙醇醛的分子动态结构

以量子化学计算作为起点, 为最简单的糖类分子——乙醇醛开发了两套分子力学力场参数: 基于肽类的力场和基于醛类的力场. 分子动力学模拟结果表明, 所开发的类醛力场参数能够较好地描述乙醇醛分子在水中的结构以及水分子在其周围的分布. 通过瞬时简正模式分析, 得到了3N-6个模式的瞬时振动频率和振动跃迁偶极矩等振动光谱参数的统计分布及其相关性. 结合量子化学计算和分子动力学模拟对生物分子体系的多元振动光谱参数进行预测和评估, 为从化学键水平出发模拟宽带飞秒二维红外相关光谱提供了一个新方法.     

Harmonic IR spectra from empirical force fields and ab initio dipole-moment derivatives

Infrared spectra simulations require ab initio techniques to get reliable intensities. On the other hand, recent force fields can provide accurate molecular geometries and frequencies. Therefore, it is suggested that these new force fields could be used to simulate infrared spectra, dipole-moment surfaces being described at high levels of theory in order to get satisfactory intensities. As pointed out, for a system with N atoms, the cost of such a simulation would be reduced N-fold with respect to standard quantum approaches. Preliminary calculations based on this scheme are reported here. Encouraging results are obtained since no significant lost of accuracy is noted on going from the ab initio to the molecular mechanics potential energy surface. © 1998 John Wiley & Sons, Inc. Int J Quant Chem 69: 705–711, 1998     

Advance of computational technology for simulating solvent extraction.

Due to recent significant enhancement of computer performance as well as computational techniques, molecular modeling and molecular simulations using computational chemistry can be achieved at the level of practical applications. Even in solvent extraction, the application of computational chemistry to simulations of extraction processes and the molecular design of high-performance extracting agents have gradually been increasing during the last decade. With combining the quantitative structure-property relationship between the molecule properties calculated by the computational chemistry methods and the thermodynamic properties obtained from experiments, researchers can precisely predict the next-generation of extracting agents and novel extraction processes. In this review, the concept of computational chemistry, such as molecular mechanics, molecular orbitals and molecular dynamics calculations, frequently used in the filed of solvent extraction, are outlined. Our systematic research on the solvent-extraction process utilizing MM, MO and MD calculations is also presented.     

A test on peptide stability of AMBER force fields with implicit solvation

We used replica exchange molecular dynamics (REMD) simulations to evaluate four different AMBER force fields and three different implicit solvent models. Our aim was to determine if these physics-based models captured the correct secondary structures of two alpha-helical and two beta-peptides: the 14-mer EK helix of Baldwin and co-workers, the C-terminal helix of ribonuclease, the 16-mer C-terminal hairpin of protein G, and the trpzip2 miniprotein. The different models gave different results, but generally we found that AMBER ff96 plus the implicit solvent model of Onufriev, Bashford, and Case gave reasonable structures, and is fairly well-balanced between helix and sheet. We also observed differences in the strength of ion pairing in the solvent models, we but found that the native secondary structures were retained even when salt bridges were prevented in the conformational sampling. Overall, this work indicates that some of these all-atom physics-based force fields may be good starting points for protein folding and protein structure prediction.     

The inclusion of electrostatic hydration energies in molecular mechanics calculations

Summary The problem of including solvent effects in molecular mechanics calculations is discussed. It is argued that the neglect of charge-solvent (solvation) interactions can introduce significant errors. The finite difference Poisson-Boltzmann (FDPB) method for calculating electrostatic interactions is summarized and is used as a basis for introducing a new pairwise energy term which accounts for charge-solvent interactions. This term acts between all pairs of atoms usually considered in molecular mechanics calculations and can be easily incorporated into existing force fields. As an example, a parameterization is developed for the CHARMm force field and the results compared to the predictions of the FDPB method. An approach to the realistic incorporation of solvent screening into force fields is also outlined.     

Nonuniform charge scaling (NUCS): a practical approximation of solvent electrostatic screening in proteins

In molecular mechanics calculations, electrostatic interactions between chemical groups are usually represented by a Coulomb potential between the partial atomic charges of the groups. In aqueous solution these interactions are modified by the polarizable solvent. Although the electrostatic effects of the polarized solvent on the protein are well described by the Poisson--Boltzmann equation, its numerical solution is computationally expensive for large molecules such as proteins. The procedure of nonuniform charge scaling (NUCS) is a pragmatic approach to implicit solvation that approximates the solvent screening effect by individually scaling the partial charges on the explicit atoms of the macromolecule so as to reproduce electrostatic interaction energies obtained from an initial Poisson--Boltzmann analysis. Once the screening factors have been determined for a protein the scaled charges can be easily used in any molecular mechanics program that implements a Coulomb term. The approach is particularly suitable for minimization-based simulations, such as normal mode analysis, certain conformational reaction path or ligand binding techniques for which bulk solvent cannot be included explicitly, and for combined quantum mechanical/molecular mechanical calculations when the interface to more elaborate continuum solvent models is lacking. The method is illustrated using reaction path calculations of the Tyr 35 ring flip in the bovine pancreatic trypsin inhibitor.     

Yammp: Development of a molecular mechanics program using the modular programming method

Molecular mechanics is a fast developing discipline with new methods and potential fields appearing every year. A versatile molecular mechanics program supports many methods and potential fields that make it unavoidably large. There are problems writing and maintaining large programs with traditional methods because of data and other dependencies. Modular programming provides a solution. A program is developed as a collection of highly independent modules containing only related data structures and procedures. These entities are isolated in the module and access to them is provided through a well-defined and controlled interface. The high degree of independence circumscribes programming errors. Most of all, it reduces the cost of revising the program as only a small part of the program needs to be read and understood for each revision. We implemented a molecular mechanics program, yammp, using the modular programming method. © 1993 John Wiley & Sons, Inc.     

A New Method for Molecular Mechanics

A promising new method for optimizing molecular structures is described. In place of the terms involving bond angles and torsion angles, used in the force fields of conventional molecular mechanics, two-body central forces between atoms are used exclusively, resulting in a considerable computational advantage. The program STRFIT, using this method has been tested by comparing geometries obtained with those found using the popular molecular mechanics program MM2 (Allinger) for a variety of cyclic and acyclic molecules. For unstrained molecules, the difference in steric energy between geometries refined by STRFIT and MM2 is only a few tenths of a kilocalorie and up to about a kilocalorie for strained molecules. Geometry optimization with STRFIT, to a structure that is slightly higher in energy than the structure arrived at by MM2 starting from the same initial starting geometry, is three to eight times faster. A complete new package of programs for conveniently and rapidly doing molecular mechanics calculations is described. It includes a convenient algorithm for the input of approximate molecular structures, a rapid structure-optimizing module using a pure Central force-field approach, and a drawing program designed for use with a dot-matrix printer or a laser printer.     

Development of glycyl radical parameters for the OPLS-AA/L force field

On the basis of quantum chemical calculations C(alpha)-glycyl radical parameters have been developed for the OPLS-AA/L force field. The molecular mechanics hypersurface was fitted to the calculated quantum chemical surface by minimizing their molecular mechanics parameter dependent sum-of-squares deviations. To do this, a computer program in which the molecular mechanics energy derivatives with respect to the parameters were calculated analytically was developed, implementing the general method of Lifson and Warshel (J Chem Phys 1968, 49, 5116) for force field parameter optimization. This program, in principle, can determine the optimal parameter set in one calculation if enough representative value points on the quantum chemical potential energy surface are available and there is no linear dependency between the parameters. Some of the parameters in quantum calculations, including several new torsion types around a bond as well as angle parameters at a new central atom type, are not completely separable. Consequently, some restrictions and/or presumptions were necessary during parameter optimization. The relative OPLS-AA energies reproduced those calculated quantum chemically almost perfectly.     

Molecular mechanics force-field parameterization procedures

A set of procedures and guidelines are presented for the estimation of bond length, bond angle, and torsional potential constants for molecular mechanics force fields. The force field constants are ultimately derived by “subtracting” nonbonded molecular mechanics energies from corresponding molecular orbital energies using a model compound containing the chemical structure to be parameterized. Case study examples of bond length, bond angle, and torsional rotation force field parameterizations are presented. A general discussion of molecular mechanics force field parameterization strategy is included for reference and completeness. Finally, a curve-fitting program to generate force field parameters from raw data is given in Appendix I.