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.     

New out-of-plane angle and bond angle internal coordinates and related potential energy functions for molecular mechanics and dynamics simulations

With currently used definitions of out-of-plane angle and bond angle internal coordinates, Cartesian derivatives have singularities, at ±π/2 in the former case and π in the latter. If either of these occur during molecular mechanics or dynamics simulations, the forces are not well defined. To avoid such difficulties, we provide new out-of-plane and bond angle coordinates and associated potential energy functions that inherently avoid these singularities. The application of these coordinates is illustrated by ab initio calculations on ammonia, water, and carbon dioxide. ©1999 John Wiley & Sons, Inc. J Comput Chem 20: 1067–1084, 1999     

Parameters for molecular dynamics simulations of iron‐sulfur proteins

Iron‐sulfur proteins involved in electron transfer reactions have finely tuned redox potentials, which allow them to be highly efficient and specific. Factors such as metal center solvent exposure, interaction with charged residues, or hydrogen bonds between the ligand residues and amide backbone groups have all been pointed out to cause such specific redox potentials. Here, we derived parameters compatible with the AMBER force field for the metal centers of iron‐sulfur proteins and applied them in the molecular dynamics simulations of three iron‐sulfur proteins. We used density‐functional theory (DFT) calculations and Seminario's method for the parameterization. Parameter validation was obtained by matching structures and normal frequencies at the quantum mechanics and molecular mechanics levels of theory. Having guaranteed a correct representation of the protein coordination spheres, the amide H‐bonds and the water exposure to the ligands were analyzed. Our results for the pattern of interactions with the metal centers are consistent to those obtained by nuclear magnetic resonance spectroscopy (NMR) experiments and DFT calculations, allowing the application of molecular dynamics to the study of those proteins. © 2013 Wiley Periodicals, Inc.     

Minimum sequence requirements for selective RNA-ligand binding: a molecular mechanics algorithm using molecular dynamics and free-energy techniques

In vitro evolution techniques allow RNA molecules with unique functions to be developed. However, these techniques do not necessarily identify the simplest RNA structures for performing their functions. Determining the simplest RNA that binds to a particular ligand is currently limited to experimental protocols. Here, we introduce a molecular-mechanics based algorithm employing molecular dynamics simulations and free-energy methods to predict the minimum sequence requirements for selective ligand binding to RNA. The algorithm involves iteratively deleting nucleotides from an experimentally determined structure of an RNA-ligand complex, performing energy minimizations and molecular dynamics on each truncated structure, and assessing which truncations do not prohibit RNA binding to the ligand. The algorithm allows prediction of the effects of sequence modifications on RNA structural stability and ligand-binding energy. We have implemented the algorithm in the AMBER suite of programs, but it could be implemented in any molecular mechanics force field parameterized for nucleic acids. Test cases are presented to show the utility and accuracy of the methodology.     

A quantum mechanics/molecular mechanics study of the protein-ligand interaction for inhibitors of HIV-1 integrase

Human immunodeficiency virus type-1 integrase (HIV-1 IN) is an essential enzyme for effective viral replication. Diketo acids such as L-731,988 and S-1360 are potent and selective inhibitors of HIV-1 IN. In this study, we used molecular dynamics simulations, within the hybrid quantum mechanics/molecular mechanics (QM/MM) approach, to determine the protein-ligand interaction energy between HIV-1 IN and L-731,988 and 10 of its derivatives and analogues. This hybrid methodology has the advantage that it includes quantum effects such as ligand polarisation upon binding, which can be very important when highly polarisable groups are embedded in anisotropic environments, as for example in metal-containing active sites. Furthermore, an energy decomposition analysis was performed to determine the contributions of individual residues to the enzyme-inhibitor interactions on averaged structures obtained from rather extensive conformational sampling. Analysis of the results reveals first that there is a correlation between protein-ligand interaction energy and experimental strand transfer into human chromosomes and secondly that the Asn-155, Lys-156 and Lys-159 residues and the Mg(2+) ion are crucial to anti-HIV IN activity. These results may explain the available experimental data.     

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.     

Comparison of properties of Aib-rich peptides in crystal and solution: a molecular dynamics study.

In order to study the differences of the structural properties of Aib-rich peptides in solution and in the crystalline state, molecular dynamics (MD) simulations of the Aib-containing peptide II (pBrBz-(Aib)5-Leu-(Aib)2-OMe) were performed in the crystalline state, starting from two different conformers obtained experimentally by X-ray diffraction. The structural properties as derived from X-ray crystallography (e.g., torsional angles and hydrogen bonds) are well-reproduced in both constant-volume and constant-pressure simulations, although the force-field parameters used result in a too-high density of the crystals. Through comparison with the results from previous MD and nuclear magnetic resonance (NMR) studies of the very similar peptide I (Z-(Aib)s-Leu-(Aib)2-OMe) in dimethylsulfoxide (DMSO) solution, it is found that, in the crystal simulation, the conformational distribution of peptide II is much narrower than that in the solution simulation of peptide. I. This leads to a significant difference in 3 [symbol: see text] (HN, HC alpha) coupling constant values, in agreement with experimental data, whereas the NOE intensities or proton-proton distance bounds appear insensitive to the difference in conformational distribution. For small peptides the differences between their conformational distribution in the crystalline form and in solution may be much larger than for proteins, a fact which should be kept in mind when interpreting molecular properties in the solution state by using X-ray crystallographic data.     

Paramfit: Automated optimization of force field parameters for molecular dynamics simulations

The generation of bond, angle, and torsion parameters for classical molecular dynamics force fields typically requires fitting parameters such that classical properties such as energies and gradients match precalculated quantum data for structures that scan the value of interest. We present a program, Paramfit, distributed as part of the AmberTools software package that automates and extends this fitting process, allowing for simplified parameter generation for applications ranging from single molecules to entire force fields. Paramfit implements a novel combination of a genetic and simplex algorithm to find the optimal set of parameters that replicate either quantum energy or force data. The program allows for the derivation of multiple parameters simultaneously using significantly fewer quantum calculations than previous methods, and can also fit parameters across multiple molecules with applications to force field development. Paramfit has been applied successfully to systems with a sparse number of structures, and has already proven crucial in the development of the Assisted Model Building with Energy Refinement Lipid14 force field. © 2014 Wiley Periodicals, Inc.     

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.     

MkVsites: A tool for creating GROMACS virtual sites parameters to increase performance in all-atom molecular dynamics simulations

The absolute performance of any all-atom molecular dynamics simulation is typically limited by the length of the individual timesteps taken when integrating the equations of motion. In the GROMACS simulation software, it has for a long time been possible to use so-called virtual sites to increase the length of the timestep, resulting in a large gain of simulation efficiency. Up until now, support for this approach has in practice been limited to the standard 20 amino acids however, shrinking the applicability domain of virtual sites. MkVsites is a set of python tools which provides a convenient way to obtain all parameters necessary to use virtual sites for virtually any molecules in a simulation. Required as input to MkVsites is the molecular topology of the molecule(s) in question, along with a specification of where to find the parent force field. As such, MkVsites can be a very valuable tool suite for anyone who is routinely using GROMACS for the simulation of molecular systems.     

A computationally inexpensive modification of the point dipole electrostatic polarization model for molecular simulations

We present an approximation, which allows reduction of computational resources needed to explicitly incorporate electrostatic polarization into molecular simulations utilizing empirical force fields. The proposed method is employed to compute three-body energies of molecular complexes with dipolar electrostatic probes, gas-phase dimerization energies, and pure liquid properties for five systems that are important in biophysical and organic simulations-water, methanol, methylamine, methanethiol, and acetamide. In all the cases, the three-body energies agreed with high level ab initio data within 0.07 kcal/mol, dimerization energies-within 0.43 kcal/mol (except for the special case of the CH(3)SH), and computed heats of vaporization and densities differed from the experimental results by less than 2%. Moreover, because the presented method allows a significant reduction in computational cost, we were able to carry out the liquid-state calculations with Monte Carlo technique. Comparison with the full-scale point dipole method showed that the computational time was reduced by 3.5 to more than 20 times, depending on the system in hand and on the desired level of the full-scale model accuracy, while the difference in energetic results between the full-scale and the presented approximate model was not great in the most cases. Comparison with the nonpolarizable OPLS-AA force field for all the substances involved and with the polarizable POL3 and q90 models for water and methanol, respectively, demonstrates that the presented technique allows reduction of computational cost with no sacrifice of accuracy. We hope that the proposed method will be of benefit to research employing molecular modeling technique in the biophysical and physical organic chemistry areas.     

Extending the treatment of backbone energetics in protein force fields: limitations of gas-phase quantum mechanics in reproducing protein conformational distributions in molecular dynamics simulations

Computational studies of proteins based on empirical force fields represent a powerful tool to obtain structure-function relationships at an atomic level, and are central in current efforts to solve the protein folding problem. The results from studies applying these tools are, however, dependent on the quality of the force fields used. In particular, accurate treatment of the peptide backbone is crucial to achieve representative conformational distributions in simulation studies. To improve the treatment of the peptide backbone, quantum mechanical (QM) and molecular mechanical (MM) calculations were undertaken on the alanine, glycine, and proline dipeptides, and the results from these calculations were combined with molecular dynamics (MD) simulations of proteins in crystal and aqueous environments. QM potential energy maps of the alanine and glycine dipeptides at the LMP2/cc-pVxZ//MP2/6-31G* levels, where x = D, T, and Q, were determined, and are compared to available QM studies on these molecules. The LMP2/cc-pVQZ//MP2/6-31G* energy surfaces for all three dipeptides were then used to improve the MM treatment of the dipeptides. These improvements included additional parameter optimization via Monte Carlo simulated annealing and extension of the potential energy function to contain peptide backbone phi, psi dihedral crossterms or a phi, psi grid-based energy correction term. Simultaneously, MD simulations of up to seven proteins in their crystalline environments were used to validate the force field enhancements. Comparison with QM and crystallographic data showed that an additional optimization of the phi, psi dihedral parameters along with the grid-based energy correction were required to yield significant improvements over the CHARMM22 force field. However, systematic deviations in the treatment of phi and psi in the helical and sheet regions were evident. Accordingly, empirical adjustments were made to the grid-based energy correction for alanine and glycine to account for these systematic differences. These adjustments lead to greater deviations from QM data for the two dipeptides but also yielded improved agreement with experimental crystallographic data. These improvements enhance the quality of the CHARMM force field in treating proteins. This extension of the potential energy function is anticipated to facilitate improved treatment of biological macromolecules via MM approaches in general.     

Hybrid particle‐field molecular dynamics simulations: Parallelization and benchmarks

The parallel implementation of a recently developed hybrid scheme for molecular dynamics (MD) simulations (Milano and Kawakatsu, J Chem Phys 2009, 130, 214106) where self‐consistent field theory (SCF) and particle models are combined is described. Because of the peculiar formulation of the hybrid method, considering single particles interacting with density fields, the most computationally expensive part of the hybrid particle‐field MD simulation can be efficiently parallelized using a straightforward particle decomposition algorithm. Benchmarks of simulations, including comparisons of serial MD and MD‐SCF program profiles, serial MD‐SCF and parallel MD‐SCF program profiles, and parallel benchmarks compared with efficient MD program GROMACS 4.5.4 are tested and reported. The results of benchmarks indicate that the proposed parallelization scheme is very efficient and opens the way to molecular simulations of large scale systems with reasonable computational costs. © 2012 Wiley Periodicals, Inc.     

Coordination number of zinc ions in the phosphotriesterase active site by molecular dynamics and quantum mechanics

We have run several molecular dynamics (MD) simulations on zinc-containing phosphotriesterase (PTE) with two bound substrates, sarin and paraoxon, and with the substrate analog diethyl 4-methylbenzylphosphonate. A standard nonbonded model was employed to treat the zinc ions with the commonly used charge of +2. In all the trajectories, we observed a tightly bound water (TBW) molecule in the active site that was coordinated to the less buried zinc ion. The phosphoryl oxygen of the substrate/inhibitor was found to be coordinated to the same zinc ion so that, considering all ligands, the less buried zinc was hexa-coordinated. The hexa-coordination of this zinc ion was not seen in the deposited X-ray pdb files for PTE. Several additional MD simulations were then performed using different charges (+1, +1.5) on the zinc ions, along with ab initio and density functional theory (DFT) calculations, to evaluate the following possibilities: the crystal diffraction data were not correctly interpreted; the hexa-coordinated zinc ion in PTE is only present in solution and not in the crystal; and the hexa-coordinated zinc ion in PTE is an artifact of the force field used. A charge of +1.5 leads to a coordination number (CN) of 5 on both zinc ions, which is consistent with the results from ab initio and DFT calculations and with the latest high resolution X-ray crystal structure. The commonly used charge of +2 produces a CN of 6 on the less buried zinc. The CN on the more buried zinc ion is 5 when the substrate/inhibitor is present in the simulation, and increases to 6 when the substrate/inhibitor is removed prior to the simulation. The results of both of the MD and quantum mechanical calculations lead to the conclusion that the zinc ions in the PTE active site are both penta-coordinated, and that the MD simulations performed with the charge of +2 overestimate the CN of the zinc ions in the PTE active site. The overall protein structures in the simulations remain unaffected by the change in zinc charge from +2 to +1.5. The results also suggest that the charge +1.5 is the most appropriate for the molecular dynamics simulations on zinc-containing PTE when a nonbonded model is used and no global thermodynamic conclusion is sought. We also show that the standard nonbonded model is not able to properly treat the CN and energy at the same time. A preliminary, promising charge-transfer model is discussed with the use of the zinc charge of +1.5.     

Estimating the temperature dependence of peptide folding entropies and free enthalpies from total energies in molecular dynamics simulations

The temperature dependence of thermodynamic quantities, such as heat capacity, entropy and free enthalpy, may be obtained by using well-known equations that relate these quantities to the enthalpy of the molecular system of interest at a range of temperatures. In turn, the enthalpy of a molecular system can be estimated from molecular dynamics simulations of an appropriate model. To demonstrate this, we have investigated the temperature dependence of the enthalpy, heat capacity, entropy and free enthalpy of a system that consists of a beta-heptapeptide in methanol and have used the statistical mechanics relationships to describe the thermodynamics of the folding/unfolding equilibrium of the peptide. The results illustrate the power of current molecular simulation force fields and techniques in establishing the link between thermodynamic quantities and conformational distributions.     

MDAnalysis: A toolkit for the analysis of molecular dynamics simulations

MDAnalysis is an object‐oriented library for structural and temporal analysis of molecular dynamics (MD) simulation trajectories and individual protein structures. It is written in the Python language with some performance‐critical code in C. It uses the powerful NumPy package to expose trajectory data as fast and efficient NumPy arrays. It has been tested on systems of millions of particles. Many common file formats of simulation packages including CHARMM, Gromacs, Amber, and NAMD and the Protein Data Bank format can be read and written. Atoms can be selected with a syntax similar to CHARMM's powerful selection commands. MDAnalysis enables both novice and experienced programmers to rapidly write their own analytical tools and access data stored in trajectories in an easily accessible manner that facilitates interactive explorative analysis. MDAnalysis has been tested on and works for most Unix‐based platforms such as Linux and Mac OS X. It is freely available under the GNU General Public License from http://mdanalysis.googlecode.com . © 2011 Wiley Periodicals, Inc. J Comput Chem 2011     

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.     

Conformational behavior of antineoplastic peptides Dolastatin 10 and Dolastatin 15 from Monte Carlo and molecular dynamics simulations

In this work, we report the results of a combined Monte Carlo (MC) and molecular dynamics (MD) approach applied to the conformational study of natural potent mitosis inhibitor Dolastatin 15. A GRID and CLogP analysis has been performed with the aim to provide some hints to the interpretation of the different behavior of Dolastatin 15 and Dolastatin 10, which has been examined in a former study. © 2006 Wiley Periodicals, Inc. Int J Quantum Chem, 2007