A point-charge force field for molecular mechanics simulations of proteins based on condensed-phase quantum mechanical calculations

Molecular mechanics models have been applied extensively to study the dynamics of proteins and nucleic acids. Here we report the development of a third-generation point-charge all-atom force field for proteins. Following the earlier approach of Cornell et al., the charge set was obtained by fitting to the electrostatic potentials of dipeptides calculated using B3LYP/cc-pVTZ//HF/6-31G** quantum mechanical methods. The main-chain torsion parameters were obtained by fitting to the energy profiles of Ace-Ala-Nme and Ace-Gly-Nme di-peptides calculated using MP2/cc-pVTZ//HF/6-31G** quantum mechanical methods. All other parameters were taken from the existing AMBER data base. The major departure from previous force fields is that all quantum mechanical calculations were done in the condensed phase with continuum solvent models and an effective dielectric constant of epsilon = 4. We anticipate that this force field parameter set will address certain critical short comings of previous force fields in condensed-phase simulations of proteins. Initial tests on peptides demonstrated a high-degree of similarity between the calculated and the statistically measured Ramanchandran maps for both Ace-Gly-Nme and Ace-Ala-Nme di-peptides. Some highlights of our results include (1) well-preserved balance between the extended and helical region distributions, and (2) favorable type-II poly-proline helical region in agreement with recent experiments. Backward compatibility between the new and Cornell et al. charge sets, as judged by overall agreement between dipole moments, allows a smooth transition to the new force field in the area of ligand-binding calculations. Test simulations on a large set of proteins are also discussed.     

Application of molecular dynamics simulations in molecular property prediction II: diffusion coefficient

In this work, we have evaluated how well the general assisted model building with energy refinement (AMBER) force field performs in studying the dynamic properties of liquids. Diffusion coefficients (D) have been predicted for 17 solvents, five organic compounds in aqueous solutions, four proteins in aqueous solutions, and nine organic compounds in nonaqueous solutions. An efficient sampling strategy has been proposed and tested in the calculation of the diffusion coefficients of solutes in solutions. There are two major findings of this study. First of all, the diffusion coefficients of organic solutes in aqueous solution can be well predicted: the average unsigned errors and the root mean square errors are 0.137 and 0.171 × 10(-5) cm(-2) s(-1), respectively. Second, although the absolute values of D cannot be predicted, good correlations have been achieved for eight organic solvents with experimental data (R(2) = 0.784), four proteins in aqueous solutions (R(2) = 0.996), and nine organic compounds in nonaqueous solutions (R(2) = 0.834). The temperature dependent behaviors of three solvents, namely, TIP3P water, dimethyl sulfoxide, and cyclohexane have been studied. The major molecular dynamics (MD) settings, such as the sizes of simulation boxes and with/without wrapping the coordinates of MD snapshots into the primary simulation boxes have been explored. We have concluded that our sampling strategy that averaging the mean square displacement collected in multiple short-MD simulations is efficient in predicting diffusion coefficients of solutes at infinite dilution.     

A new set of molecular mechanics parameters for hydroxyproline and its use in molecular dynamics simulations of collagen-like peptides

Recently, the importance of proline ring pucker conformations in collagen has been suggested in the context of hydroxylation of prolines. The previous molecular mechanics parameters for hydroxyproline, however, do not reproduce the correct pucker preference. We have developed a new set of parameters that reproduces the correct pucker preference. Our molecular dynamics simulations of proline and hydroxyproline monomers as well as collagen-like peptides, using the new parameters, support the theory that the role of hydroxylation in collagen is to stabilize the triple helix by adjusting to the right pucker conformation (and thus the right phi angle) in the Y position.     

An extensible interface for QM/MM molecular dynamics simulations with AMBER

We present an extensible interface between the AMBER molecular dynamics (MD) software package and electronic structure software packages for quantum mechanical (QM) and mixed QM and classical molecular mechanical (MM) MD simulations within both mechanical and electronic embedding schemes. With this interface, ab initio wave function theory and density functional theory methods, as available in the supported electronic structure software packages, become available for QM/MM MD simulations with AMBER. The interface has been written in a modular fashion that allows straight forward extensions to support additional QM software packages and can easily be ported to other MD software. Data exchange between the MD and QM software is implemented by means of files and system calls or the message passing interface standard. Based on extensive tests, default settings for the supported QM packages are provided such that energy is conserved for typical QM/MM MD simulations in the microcanonical ensemble. Results for the free energy of binding of calcium ions to aspartate in aqueous solution comparing semiempirical and density functional Hamiltonians are shown to demonstrate features of this interface. © 2013 Wiley Periodicals, Inc.     

Importance of explicit smeared lone‐pairs in anisotropic polarizable molecular mechanics. Torture track angular tests for exchange‐repulsion and charge transfer contributions

A correct representation of the short‐range contributions such as exchange‐repulsion (E _rep) and charge‐transfer (E _ct) is essential for the soundness of separable, anisotropic polarizable molecular mechanics potentials. Within the context of the SIBFA procedure, this is aimed at by explicit representations of lone pairs in their expressions. It is necessary to account for their anisotropic behaviors upon performing not only in‐plane, but also out‐of‐plane, variations of a probe molecule or cation interacting with a target molecule or molecular fragment. Thus, E _rep and E _ct have to reproduce satisfactorily the corresponding anisotropies of their quantum chemical (QC) counterparts. A significant improvement of the out‐of‐plane dependencies was enabled when the sp² and sp localized lone‐pairs are, even though to a limited extent, delocalized on both sides of the plane, above and below the atom bearer but at the closely similar angles as the in‐plane lone pair. We report calibration and validation tests on a series of monoligated complexes of a probe Zn(II) cation with several biochemically relevant ligands. Validations are then performed on several polyligated Zn(II) complexes found in the recognition sites of Zn‐metalloproteins. Such calibrations and validations are extended to representative monoligated and polyligated complexes of Mg(II) and Ca(II). It is emphasized that the calibration of all three cations was for each ΔE contribution done on a small training set bearing on a limited number of representative N , O , and S monoligated complexes. Owing to the separable nature of ΔE , a secure transferability is enabled to a diversity of polyligated complexes. For these the relative errors with respect to the target ΔE (QC) values are generally < 3%. Overall, the article proposes a full set of benchmarks that could be useful for force field developers. © 2017 Wiley Periodicals, Inc.     

Ergodicity and model quality in template‐restrained canonical and temperature/Hamiltonian replica exchange coarse‐grained molecular dynamics simulations of proteins

Molecular simulations restrained to single or multiple templates are commonly used in protein‐structure modeling. However, the restraints introduce additional barriers, thus impairing the ergodicity of simulations, which can affect the quality of the resulting models. In this work, the effect of restraint types and simulation schemes on ergodicity and model quality was investigated by performing template‐restrained canonical molecular dynamics (MD), multiplexed replica‐exchange molecular dynamics, and Hamiltonian replica exchange molecular dynamics (HREMD) simulations with the coarse‐grained UNRES force field on nine selected proteins, with pseudo‐harmonic log‐Gaussian (unbounded) or Lorentzian (bounded) restraint functions. The best ergodicity was exhibited by HREMD. It has been found that non‐ergodicity does not affect model quality if good templates are used to generate restraints. However, when poor‐quality restraints not covering the entire protein are used, the improved ergodicity of HREMD can lead to significantly improved protein models. © 2017 Wiley Periodicals, Inc.     

Introduction of steered molecular dynamics into UNRES coarse‐grained simulations package

In this article, an implementation of steered molecular dynamics (SMD) in coarse‐grain UNited RESidue (UNRES) simulations package is presented. Two variants of SMD have been implemented: with a constant force and a constant velocity. The huge advantage of SMD implementation in the UNRES force field is that it allows to pull with the speed significantly lower than the accessible pulling speed in simulations with all‐atom representation of a system, with respect to a reasonable computational time. Therefore, obtaining pulling speed closer to those which appear in the atomic force spectroscopy is possible. The newly implemented method has been tested for behavior in a microcanonical run to verify the influence of introduction of artificial constrains on keeping total energy of the system. Moreover, as time dependent artificial force was introduced, the thermostat behavior was tested. The new method was also tested via unfolding of the Fn3 domain of human contactin 1 protein and the I27 titin domain. Obtained results were compared with Gø‐like force field, all‐atom force field, and experimental results. © 2017 Wiley Periodicals, Inc.     

Atomistic insight into chondroitin‐6‐sulfate glycosaminoglycan chain through quantum mechanics calculations and molecular dynamics simulation

Chondroitin‐6‐sulfate (C6S) is a glycosaminoglycan (GAG) constituent in the extracellular matrix, which participates actively in crucial biological processes, as well as in various pathological conditions, such as atherosclerosis and cancer. Molecular interactions involving the C6S chain are therefore of considerable interest. A computational model for atomistic simulation was built. This work describes the design and validation of a force field for a C6S dodecasaccharide chain. The results of an extensive molecular dynamics simulation performed with the new force field provide a novel insight into the structure and dynamics of the C6S chain. The intramolecular H‐bonds in the disaccharide linkage region are suggested to play a major role in determining the chain structural dynamics. Moreover, the unravelling of an additional H‐bond involving the sulfate groups in C6S is interesting as changes in sulfation have been claimed to be an important factor in several diseases. The force field will prove useful for future studies of crucial interactions between C6S and various nanoassemblies. It can also be used as a basis for modeling of other GAGs. © 2010 Wiley Periodicals, Inc. J Comput Chem, 2010     

Substrate‐Guided Front‐Face Reaction Revealed by Combined Structural Snapshots and Metadynamics for the Polypeptide N‐Acetylgalactosaminyltransferase 2

The retaining glycosyltransferase GalNAc‐T2 is a member of a large family of human polypeptide GalNAc‐transferases that is responsible for the post‐translational modification of many cell‐surface proteins. By the use of combined structural and computational approaches, we provide the first set of structural snapshots of the enzyme during the catalytic cycle and combine these with quantum‐mechanics/molecular‐mechanics (QM/MM) metadynamics to unravel the catalytic mechanism of this retaining enzyme at the atomic‐electronic level of detail. Our study provides a detailed structural rationale for an ordered bi–bi kinetic mechanism and reveals critical aspects of substrate recognition, which dictate the specificity for acceptor Thr versus Ser residues and enforce a front‐face S_Ni‐type reaction in which the substrate N‐acetyl sugar substituent coordinates efficient glycosyl transfer.     

Strike a balance: optimization of backbone torsion parameters of AMBER polarizable force field for simulations of proteins and peptides

Based on the AMBER polarizable model (ff02), we have re-optimized the parameters related to the main-chain (Phi, Psi) torsion angles by fitting to the Boltzmann-weighted average quantum mechanical (QM) energies of the important regions (i.e., beta, P(II), alpha(R), and alpha(L) regions). Following the naming convention of the AMBER force field series, this release will be called ff02pol.rl The force field has been assessed both by energetic comparison against the QM data and by the replica exchange molecular dynamics simulations of short alanine peptides in water. For Ace-Ala-Nme, the simulated populations in the beta, P(II) and alpha(R) regions were approximately 30, 43, and 26%, respectively. For Ace-(Ala)(7)-Nme, the populations in these three regions were approximately 24, 49, and 26%. Both were in qualitative agreement with the NMR and CD experimental conclusions. In comparison with the previous force field, ff02pol.rl demonstrated good balance among these three important regions. The optimized torsion parameters, together with those in ff02, allow us to carry out simulations on proteins and peptides with the consideration of polarization.     

Nucleic acid reactivity: Challenges for next‐generation semiempirical quantum models

Semiempirical quantum models are routinely used to study mechanisms of RNA catalysis and phosphoryl transfer reactions using combined quantum mechanical (QM)/molecular mechanical methods. Herein, we provide a broad assessment of the performance of existing semiempirical quantum models to describe nucleic acid structure and reactivity to quantify their limitations and guide the development of next‐generation quantum models with improved accuracy. Neglect of diatomic differential overlap and self‐consistent density‐functional tight‐binding semiempirical models are evaluated against high‐level QM benchmark calculations for seven biologically important datasets. The datasets include: proton affinities, polarizabilities, nucleobase dimer interactions, dimethyl phosphate anion, nucleoside sugar and glycosidic torsion conformations, and RNA phosphoryl transfer model reactions. As an additional baseline, comparisons are made with several commonly used density‐functional models, including M062X and B3LYP (in some cases with dispersion corrections). The results show that, among the semiempirical models examined, the AM1/d‐PhoT model is the most robust at predicting proton affinities. AM1/d‐PhoT and DFTB3‐3ob/OPhyd reproduce the MP2 potential energy surfaces of 6 associative RNA phosphoryl transfer model reactions reasonably well. Further, a recently developed linear‐scaling “modified divide‐and‐conquer” model exhibits the most accurate results for binding energies of both hydrogen bonded and stacked nucleobase dimers. The semiempirical models considered here are shown to underestimate the isotropic polarizabilities of neutral molecules by approximately 30%. The semiempirical models also fail to adequately describe torsion profiles for the dimethyl phosphate anion, the nucleoside sugar ring puckers, and the rotations about the nucleoside glycosidic bond. The modeling of pentavalent phosphorus, particularly with thio substitutions often used experimentally as mechanistic probes, was problematic for all of the models considered. Analysis of the strengths and weakness of the models suggests that the creation of robust next‐generation models should emphasize the improvement of relative conformational energies and barriers, and nonbonded interactions. © 2015 Wiley Periodicals, Inc.     

Parameterization of Ca+2-protein interactions for molecular dynamics simulations

Molecular dynamics simulations of Ca+2 ions near protein were performed with three force fields: GROMOS96, OPLS-AA, and CHARMM22. The simulations reveal major, force-field dependent, inconsistencies in the interaction between the Ca+2 ions with the protein. The variations are attributed to the nonbonded parameterizations of the Ca+2-carboxylates interactions. The simulations results were compared to experimental data, using the Ca+2-HCOO- equilibrium as a model. The OPLS-AA force field grossly overestimates the binding affinity of the Ca+2 ions to the carboxylate whereas the GROMOS96 and CHARMM22 force fields underestimate the stability of the complex. Optimization of the Lennard-Jones parameters for the Ca+2-carboxylate interactions were carried out, yielding new parameters which reproduce experimental data.     

Ab initio quantum chemistry and molecular dynamics simulations studies of LiPF6/poly(ethylene oxide) interactions

Ab initio and molecular mechanics studies of LiPF₆ and the interaction of the salt with the poly(ethylene oxide) (PEO) oligomer dimethylether have been performed. Optimized geometries and energies of Li⁺/PF₆^? complexes obtained from quantum chemistry revealed a preference for C_3V symmetry structures for Li⁺–P separations under 2.8 Å, C_2V symmetry for Li⁺–P in the range of 2.8–3.3 Å and C_4V symmetry for Li⁺–P separations larger than 3.3 Å. Electron correlation effects were found to make an insignificant contribution to binding in the Li⁺/PF₆^? complex. By contrast, analogous studies of PF₆^?/PF₆^? and PF₆^?/dimethyl ether complexes revealed important contributions of electron correlation to the complex interaction energy. A molecular mechanics force field for simulations of PEO/LiPF₆ melts was parameterized to reproduce the geometries and energies of Li⁺/PF₆^?, PF₆^?/PF₆^?, PF₆^?/dimethylether complexes. Molecular dynamics simulations of PEO/LiPF₆ melts were performed to validate this quantum chemistry‐based force field. Accurate reproduction of the increase in solution density with addition of salt was found while the electrical conductivity of PEO/LiPF₆ solutions was found to be within an order of magnitude of the experimental values. © 2001 John Wiley & Sons, Inc. J Comput Chem 22: 641–654, 2001     

Revised CHARMM force field parameters for iron‐containing cofactors of photosystem II

Photosystem II is a complex protein–cofactor machinery that splits water molecules into molecular oxygen, protons, and electrons. All‐atom molecular dynamics simulations have the potential to contribute to our general understanding of how photosystem II works. To perform reliable all‐atom simulations, we need accurate force field parameters for the cofactor molecules. We present here CHARMM bonded and non‐bonded parameters for the iron‐containing cofactors of photosystem II that include a six‐coordinated heme moiety coordinated by two histidine groups, and a non‐heme iron complex coordinated by bicarbonate and four histidines. The force field parameters presented here give water interaction energies and geometries in good agreement with the quantum mechanical target data. © 2017 Wiley Periodicals, Inc.     

A box-counting-based algorithm for computing Shannon entropy in molecular dynamics simulations

A box-counting-based algorithm (SEBC) has been developed for the numerical computation of the Shannon entropy from samples of continuous functions. Its performance was tested by applying it to several samples of known continuous distribution functions. The results obtained with SEBC reproduced those obtained by analytical or numerical integration. SEBC was also employed for computing the Shannon entropies of the steric energy, Sh(E(S)), of several amino acids from their in vacuo NVE molecular dynamics simulations using the AMBER-4 force field. The results obtained correlate linearly with the experimental standard thermodynamic entropies of these compounds. This work points to the possibility of introducing straightforward and reliable calculations of thermodynamic entropies from empirical linear relationships with Sh(E(S)) obtained from MD simulations.     

CHARMM fluctuating charge force field for proteins: II protein/solvent properties from molecular dynamics simulations using a nonadditive electrostatic model

A fluctuating charge (FQ) force field is applied to molecular dynamics simulations for six small proteins in explicit polarizable solvent represented by the TIP4P-FQ potential. The proteins include 1FSV, 1ENH, 1PGB, 1VII, 1H8K, and 1CRN, representing both helical and beta-sheet secondary structural elements. Constant pressure and temperature (NPT) molecular dynamics simulations are performed on time scales of several nanoseconds, the longest simulations yet reported using explicitly polarizable all-atom empirical potentials (for both solvent and protein) in the condensed phase. In terms of structure, the FQ force field allows deviations from native structure up to 2.5 A (with a range of 1.0 to 2.5 A). This is commensurate to the performance of the CHARMM22 nonpolarizable model and other currently existing polarizable models. Importantly, secondary structural elements maintain native structure in general to within 1 A (both helix and beta-strands), again in good agreement with the nonpolarizable case. In qualitative agreement with QM/MM ab initio dynamics on crambin (Liu et al. Proteins 2001, 44, 484), there is a sequence dependence of average condensed phase atomic charge for all proteins, a dependence one would anticipate considering the differing chemical environments around individual atoms; this is a subtle quantum mechanical feature captured in the FQ model but absent in current state-of-the-art nonpolarizable models. Furthermore, there is a mutual polarization of solvent and protein in the condensed phase. Solvent dipole moment distributions within the first and second solvation shells around the protein display a shift towards higher dipole moments (increases on the order of 0.2-0.3 Debye) relative to the bulk; protein polarization is manifested via the enhanced condensed phase charges of typical polar atoms such as backbone carbonyl oxygens, amide nitrogens, and amide hydrogens. Finally, to enlarge the sample set of proteins, gas-phase minimizations and 1 ps constant temperature simulations are performed on various-sized proteins to compare to earlier work by Kaminsky et al. (J Comp Chem 2002, 23, 1515). The present work establishes the feasibility of applying a fully polarizable force field for protein simulations and demonstrates the approach employed in extending the CHARMM force field to include these effects.     

Force-field development and molecular dynamics simulations of ferrocene-peptide conjugates as a scaffold for hydrogenase mimics

The increasing importance of hydrogenase enzymes in the new energy research field has led us to examine the structure and dynamics of potential hydrogenase mimics, based on a ferrocene-peptide scaffold, using molecular dynamics (MD) simulations. To enable this MD study, a molecular mechanics force field for ferrocene-bearing peptides was developed and implemented in the CHARMM simulation package, thus extending the usefulness of the package into peptide-bioorganometallic chemistry. Using the automated frequency-matching method (AFMM), optimized intramolecular force-field parameters were generated through quantum chemical reference normal modes. The partial charges for ferrocene were derived by fitting point charges to quantum-chemically computed electrostatic potentials. The force field was tested against experimental X-ray crystal structures of dipeptide derivatives of ferrocene-1,1'-dicarboxylic acid. The calculations reproduce accurately the molecular geometries, including the characteristic C2-symmetrical intramolecular hydrogen-bonding pattern, that were stable over 0.1 micros MD simulations. The crystal packing properties of ferrocene-1-(D)alanine-(D)proline-1'-(D)alanine-(D)proline were also accurately reproduced. The lattice parameters of this crystal were conserved during a 0.1 micros MD simulation and match the experimental values almost exactly. Simulations of the peptides in dichloromethane are also in good agreement with experimental NMR and circular dichroism (CD) data in solution. The developed force field was used to perform MD simulations on novel, as yet unsynthesized peptide fragments that surround the active site of [Ni-Fe] hydrogenase. The results of this simulation lead us to propose an improved design for synthetic peptide-based hydrogenase models. The presented MD simulation results of metallocenes thereby provide a convincing validation of our proposal to use ferrocene-peptides as minimal enzyme mimics.     

Explicit proton transfer in classical molecular dynamics simulations

We present Hydrogen Dynamics (HYDYN), a method that allows explicit proton transfer in classical force field molecular dynamics simulations at thermodynamic equilibrium. HYDYN reproduces the characteristic properties of the excess proton in water, from the special pair dance, to the continuous fluctuation between the limiting Eigen and Zundel complexes, and the water reorientation beyond the first solvation layer. Advantages of HYDYN with respect to existing methods are computational efficiency, microscopic reversibility, and easy parameterization for any force field. © 2014 Wiley Periodicals, Inc.     

The adaptive buffered force QM/MM method in the CP2K and AMBER software packages

The implementation and validation of the adaptive buffered force (AdBF) quantum‐mechanics/molecular‐mechanics (QM/MM) method in two popular packages, CP2K and AMBER are presented. The implementations build on the existing QM/MM functionality in each code, extending it to allow for redefinition of the QM and MM regions during the simulation and reducing QM‐MM interface errors by discarding forces near the boundary according to the buffered force‐mixing approach. New adaptive thermostats, needed by force‐mixing methods, are also implemented. Different variants of the method are benchmarked by simulating the structure of bulk water, water autoprotolysis in the presence of zinc and dimethyl‐phosphate hydrolysis using various semiempirical Hamiltonians and density functional theory as the QM model. It is shown that with suitable parameters, based on force convergence tests, the AdBF QM/MM scheme can provide an accurate approximation of the structure in the dynamical QM region matching the corresponding fully QM simulations, as well as reproducing the correct energetics in all cases. Adaptive unbuffered force‐mixing and adaptive conventional QM/MM methods also provide reasonable results for some systems, but are more likely to suffer from instabilities and inaccuracies. © 2015 The Authors. Journal of Computational Chemistry Published by Wiley Periodicals, Inc.