A force field for monosaccharides and (1 → 4) linked polysaccharides

A force field for monosaccharides that can be extended to (1 → 4) linked polysaccharides has been developed for the AMBER potential function. The resulting force field is consistent with the existing AMBER force field for proteins and nucleic acids. Modifications to the standard AMBER OH force constant and to the Lennard-Jones parameters were made. Furthermore, a 10–12 nonbonded term was included between the hydroxyl hydrogen of the saccharide and the water oxygen (TIP3P, SPC/E, etc.) to reproduce better the water–saccharide intermolecular distances. STO-3G electrostatic potential (ESP) charges were used to represent the electrostatic interactions between the saccharide and its surrounding environment. To obtain charges for polysaccharides, a scheme was developed to piece together saccharide residues through 1 → 4 connections while still retaining a net neutral charge on the molecule as a whole. Free energy perturbation (FEP) simulations of D -glucose and D -mannose in water were performed to test the resulting force field. The FEP simulations demonstrate that AMBER overestimates intramolecular interaction energies, suggesting that further improvements are needed in this part of the force field. To test further the reliability of the parameters, a molecular dynamics (MD) simulation of α-D -glucose in water was also performed. The MD simulation was able to produce structural and conformational results that are in accord with experimental evidence and previous theoretical results. Finally, a relaxed conformational map of β-maltose was assembled and it was found that the present force field is consistent with available theoretical and experimental results. © 1994 by John Wiley & Sons, Inc.     

On achieving experimental accuracy from molecular dynamics simulations of flexible molecules: aqueous glycerol

The rotational isomeric states (RIS) of glycerol at infinite dilution have been characterized in the aqueous phase via a 1 micros conventional molecular dynamics (MD) simulation, a 40 ns enhanced sampling replica exchange molecular dynamics (REMD) simulation, and a reevaluation of the experimental NMR data. The MD and REMD simulations employed the GLYCAM06/AMBER force field with explicit treatment of solvation. The shorter time scale of the REMD sampling method gave rise to RIS and theoretical scalar 3J(HH) coupling constants that were comparable to those from the much longer traditional MD simulation. The 3J(HH) coupling constants computed from the MD methods were in excellent agreement with those observed experimentally. Despite the agreement between the computed and the experimental J-values, there were variations between the rotamer populations computed directly from the MD data and those derived from the experimental NMR data. The experimentally derived populations were determined utilizing limiting J-values from an analysis of NMR data from substituted ethane molecules and may not be completely appropriate for application in more complex molecules, such as glycerol. Here, new limiting J-values have been derived via a combined MD and quantum mechanical approach and were used to decompose the experimental 3J(HH) coupling constants into population distributions for the glycerol RIS.     

Nonreactive molecular dynamics force field for crystalline hexahydro-1,3,5-trinitro-1,3,5 triazine

An empirical nonreactive force field has been developed for molecular dynamics (MD)/Monte Carlo simulation of the formation, diffusion, and agglomeration of point defects in the crystal lattice of the alpha modification of hexahydro-1,3,5-trinitro-1,3,5 triazine (RDX) using flexible molecules. Bond stretching and angle bending are represented by Morse and harmonic functions, and torsion by a truncated cosine series. Nonbonded interactions, both inter- and intramolecular, are described by Buckingham potentials separately parametrized. Intermolecular electrostatic interactions are treated via a Coulomb term coupled with a smooth 15.0 A cutoff radius. Parameters were taken in part from earlier published works and were determined partly by fitting to known molecular and crystal properties of RDX. In MD simulations at constant pressure and temperature, the model was able to stabilize and maintain the correct crystal structure, symmetry, and molecular conformation of alpha-RDX. Vibrational frequencies, lattice binding energy and dimensions, coefficients of thermal expansion, and several unusually short intermolecular distances are all reproduced in satisfactory agreement with experimental data.     

Communication: Quantum polarized fluctuating charge model: a practical method to include ligand polarizability in biomolecular simulations

We present a simple and practical method to include ligand electronic polarization in molecular dynamics (MD) simulation of biomolecular systems. The method involves periodically spawning quantum mechanical (QM) electrostatic potential (ESP) calculations on an extra set of computer processors using molecular coordinate snapshots from a running parallel MD simulation. The QM ESPs are evaluated for the small-molecule ligand in the presence of the electric field induced by the protein, solvent, and ion charges within the MD snapshot. Partial charges on ligand atom centers are fit through the multi-conformer restrained electrostatic potential (RESP) fit method on several successive ESPs. The RESP method was selected since it produces charges consistent with the AMBER/GAFF force-field used in the simulations. The updated charges are introduced back into the running simulation when the next snapshot is saved. The result is a simulation whose ligand partial charges continuously respond in real-time to the short-term mean electrostatic field of the evolving environment without incurring additional wall-clock time. We show that (1) by incorporating the cost of polarization back into the potential energy of the MD simulation, the algorithm conserves energy when run in the microcanonical ensemble and (2) the mean solvation free energies for 15 neutral amino acid side chains calculated with the quantum polarized fluctuating charge method and thermodynamic integration agree better with experiment relative to the Amber fixed charge force-field.     

Effective force fields for condensed phase systems from ab initio molecular dynamics simulation: a new method for force-matching

A novel least-squares fitting approach is presented to obtain classical force fields from trajectory and force databases produced by ab initio (e.g., Car-Parrinello) molecular dynamics (MD) simulations. The method was applied to derive effective nonpolarizable three-site force fields for liquid water at ambient conditions from Car-Parrinello MD simulations in the Becke-Lee-Yang-Parr approximation to the electronic density functional theory. The force-matching procedure includes a fit of short-ranged nonbonded forces, bonded forces, and atomic partial charges. The various parameterizations of the water force field differ by an enforced smooth cut-off applied to the short-ranged interaction term. These were obtained by fitting to the trajectory and force data produced by Car-Parrinello MD simulations of systems of 32 and 64 H(2)O molecules. The new water force fields were developed assuming both flexible or rigid molecular geometry. The simulated structural and self-diffusion properties of liquid water using the fitted force fields are in close agreement with those observed in the underlying Car-Parrinello MD simulations. The resulting empirical models compare to experiment much better than many conventional simple point charge (SPC) models. The fitted potential is also shown to combine well with more sophisticated intramolecular potentials. Importantly, the computational cost of the new models is comparable to that for SPC-like potentials.     

Molecular dynamics study of peptides in implicit water: ab initio folding of beta-hairpin, beta-sheet, and beta beta alpha-motif

In this communication, we have demonstrated that molecular dynamics simulations using a GB implicit solvation model with the all-atom based force field (CHARMM19) can describe the spontaneous folding of small peptides in aqueous solution. The native structures of peptides with various structural motifs (beta-hairpin, beta-sheet, and betabetaalpha-moiety) were successfully predicted within reasonable time scales by MD simulations at moderately elevated temperatures. It is expected that the present simulations provide further insight into mechanism/pathways of the peptide folding.     

Nature of minor-groove binders-DNA complexes in the gas phase

The structure of noncovalent complexes of DNA duplex with minor groove binders (mG-binders) has been analyzed by state of the art molecular dynamics (MD) simulations. More than 3.3 micros of MD trajectories (including 4 x 0.5 micros trajectories) were collected for the Dickerson's dodecamer bound to DAPI, Hoechst 33258, and Netropsin. Comparison of these trajectories with control simulations in water allowed us to determine that the extreme dehydration and partial neutralization occurring during electrospray experiments does not produce the disruption of the DNA:mG-binder complexes or the dissociation of the two strands of the duplex. Irrespective of the drug and the simulation conditions the mG-binders remains bound to the DNA near the preferential binding position in aqueous conditions. Large distortions appear in the two DNA strands, which maintain however a memory of the original DNA duplex structure in water, and a general helical-like conformation.     

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.     

Molecular dynamics study of the beta amyloid peptide of Alzheimer's disease and its divalent copper complexes

The Abeta1-42 monomer structure was assessed with a 790 ns molecular dynamics (MD) simulation, and the results were compared with the NMR experiment on Abeta10-35 and Abeta1-40. Previous theoretical work in a model of the His13-His14 region of Abeta defined the possible Cu(II) binding geometries at this site (Raffa et al. J. Biol. Inorg. Chem. 2005, 10, 887-902). MD simulations totalling almost 2 micros were also carried out on Cu(II)/Abeta1-42 systems, using the ab initio structures as templates for the copper binding site. This work finds that the copper-free Abeta1-42 system may stabilize after approximately 350 ns into a collapsed coil conformation, and we find good agreement with some, but not all, of the structural features determined experimentally for the Abeta10-35 and Abeta1-40 peptides. The results of the Cu(II)/Abeta1-42 systems are compared to the Cu(II)-free Abeta1-42 simulation.     

Convergence of Free Energy Profile of Coumarin in Lipid Bilayer

Atomistic molecular dynamics (MD) simulations of druglike molecules embedded in lipid bilayers are of considerable interest as models for drug penetration and positioning in biological membranes. Here we analyze partitioning of coumarin in dioleoylphosphatidylcholine (DOPC) bilayer, based on both multiple, unbiased 3 μs MD simulations (total length) and free energy profiles along the bilayer normal calculated by biased MD simulations (～7 μs in total). The convergences in time of free energy profiles calculated by both umbrella sampling and z-constraint techniques are thoroughly analyzed. Two sets of starting structures are also considered, one from unbiased MD simulation and the other from "pulling" coumarin along the bilayer normal. The structures obtained by pulling simulation contain water defects on the lipid bilayer surface, while those acquired from unbiased simulation have no membrane defects. The free energy profiles converge more rapidly when starting frames from unbiased simulations are used. In addition, z-constraint simulation leads to more rapid convergence than umbrella sampling, due to quicker relaxation of membrane defects. Furthermore, we show that the choice of RESP, PRODRG, or Mulliken charges considerably affects the resulting free energy profile of our model drug along the bilayer normal. We recommend using z-constraint biased MD simulations based on starting geometries acquired from unbiased MD simulations for efficient calculation of convergent free energy profiles of druglike molecules along bilayer normals. The calculation of free energy profile should start with an unbiased simulation, though the polar molecules might need a slow pulling afterward. Results obtained with the recommended simulation protocol agree well with available experimental data for two coumarin derivatives.     

Hierarchical multiscale simulation of electrokinetic transport in silica nanochannels at the point of zero charge

Effects of nanoscale confinement and partial charges that stem from quantum calculations are investigated in silica slit channels filled with 1 M KCl at the point of zero charge by using a hierarchical multiscale simulation methodology. Partial charges of both bulk and surface atoms from ab initio quantum calculations that take into account bond polarization and electronegativity are used in molecular dynamics (MD) simulations to obtain ion and water concentration profiles for channel widths of 1.1, 2.1, 2.75, and 4.1 nm. The interfacial electron density profiles of simulations matched well with that of recent X-ray reflectivity experiments. By simulating corresponding channels with no partial charges, it was observed that the partial charges affect the concentration profiles and transport properties such as diffusion coefficients and mobilities up to a distance of about 3 sigma(O)(-)(O) from the surface. Both in uncharged and partially charged cases, oscillations in concentration profiles of K(+) and Cl(-) ions give rise to an electro-osmotic flow in the presence of an external electric field, indicating the presence of an electric double layer at net zero surface charge, contrary to the expectations from classical continuum theory. I-V curves in a channel-bath system using ionic mobilities from MD simulations were significantly different for channels with and without partial charges for channel widths less than 4.1 nm.     

Nitrile and thiocyanate IR probes: molecular dynamics simulation studies

Nitrile- and thiocyanate-derivatized amino acids have been found to be useful IR probes for investigating their local electrostatic environments in proteins. To shed light on the CN stretch frequency shift and spectral lineshape change induced by interactions with hydrogen-bonding solvent molecules, we carried out both classical and quantum mechanical/molecular mechanical (QM/MM) molecular dynamics (MD) simulations for MeCN and MeSCN in water. These QM/MM and conventional force field MD simulation results were found to be inconsistent with the experimental results as well as with the high-level ab initio calculation results of MeCN-water and MeSCN-water potential energies. Thus, a new set of atomic partial charges of MeCN and MeSCN is obtained. By using the MD simulation trajectories and the electrostatic potential model recently developed, the CN and SCN stretching mode frequency trajectories were obtained and used to simulate the IR spectra. The C[Triple Bond]N frequency blueshifts of MeCN and MeSCN in water are estimated to be 9.0 and 1.9 cm(-1), respectively, in comparison with those of gas phase values. These values are found to be in reasonable agreement with the experimentally measured IR spectra of MeCN, MeSCN, beta-cyano-L-alanine, and cyanylated cysteine in water and other polar solvents.     

Ab initio protein structure prediction with force field parameters derived from water-phase quantum chemical calculation

Molecular dynamics (MD) simulations are extensively used in the study of the structures and functions of proteins. Ab initio protein structure prediction is one of the most important subjects in computational biology, and many trials have been performed using MD simulation so far. Since the results of MD simulations largely depend on the force field, reliable force field parameters are indispensable for the success of MD simulation. In this work, we have modified atom charges in a standard force field on the basis of water-phase quantum chemical calculations. The modified force field turned out appropriate for ab initio protein structure prediction by the MD simulation with the generalized Born method. Detailed analysis was performed in terms of the conformational stability of amino acid residues, the stability of secondary structure of proteins, and the accuracy for prediction of protein tertiary structure, comparing the modified force field with a standard one. The energy balance between alpha-helix and beta-sheet structures was significantly improved by the modification of charge parameters.     

An extensible and systematic force field,ESFF, for molecular modeling of organic,inorganic, and organometallic systems

ESFF is a rule-based force field designed for modeling organic, inorganic, and organometallic systems. To cover this broad range of molecular systems, ESFF was developed in an extensible and systematic manner. Several unique features were introduced including pseudoangle and a dot product function representing torsion energy terms. The partial atomic charges that are topology-dependent are determined from ab initio (DFT) calculated electronegativity and hardness for valence orbitals. The van der Waals parameters are charge-dependent, and correlated with the ionization potential for atoms in various valence states. To obtain a set of well-defined and physically meaningful parameters, ESFF employs semiempirical rules to translate atomic-based parameters to parameters typically associated with a covalent valence force field. The atomic parameters depend not only on atom type, but also on internal type, thus resulting in a more accurate force field. This article presents the theory and the method used to develop the force field. The force field has been applied to molecular simulations of a wide variety of systems including nucleic acids, peptides, hydrocarbons, porphyrins, transition metal complexes, zeolites, and organometallic compounds. Agreement with the experimental results indicates that ESFF is a valuable tool in molecular simulations for understanding and predicting both crystal and gas phase molecular structures.     

The binding domain of the HMGB1 inhibitor carbenoxolone: Theory and experiment

We present a combined computational and experimental study of the interaction of the Box A of the HMGB1 protein and carbenoxolone, an inhibitor of its pro-inflammatory activity. The computational approach consists of classical molecular dynamics (MD) simulations based on the GROMOS force field with quantum-refined (QRFF) atomic charges for the ligand. Experimental data consist of fluorescence intensities, chemical shift displacements, saturation transfer differences and intermolecular Nuclear Overhauser Enhancement signals. Good agreement is found between observations and the conformation of the ligand–protein complex resulting from QRFF–MD. In contrast, simple docking procedures and MD based on the unrefined force field provide models inconsistent with experiment. The ligand–protein binding is dominated by non-directional interactions.     

A molecular mechanics force field for lignin

A CHARMM molecular mechanics force field for lignin is derived. Parameterization is based on reproducing quantum mechanical data of model compounds. Partial atomic charges are derived using the RESP electrostatic potential fitting method supplemented by the examination of methoxybenzene:water interactions. Dihedral parameters are optimized by fitting to critical rotational potentials and bonded parameters are obtained by optimizing vibrational frequencies and normal modes. Finally, the force field is validated by performing a molecular dynamics simulation of a crystal of a lignin fragment molecule and comparing simulation-derived structural features with experimental results. Together with the existing force field for polysaccharides, this lignin force field will enable full simulations of lignocellulose.     

New all-atom force field for molecular dynamics simulation of an AlPO(4)-34 molecular sieve

A force field of the triclinic framework of AlPO(4)-34, important in methanol-hydrocarbon conversion reactions, was developed using an empirical potential function. Molecular dynamics simulation of an AlPO(4)-34 triclinic framework segment of 1216 atoms, containing the template molecules isopropylamine and water, was performed with explicit consideration of atomic charges. The average RMS difference between instantaneous positions of the framework atoms during 1 ns simulation and their positions in the structure determined from single crystal X-ray diffraction was calculated, and the average structure of the flexible framework was determined. The computed Debye-Waller factors and simulated FTIR spectra are in good agreement with the experimental data. The new force field permits detailed molecular dynamics simulations of flexible, charged aluminophosphate molecular sieves which should lead to a better understanding of the catalytic processes and the crucial role played by templating molecules.     

Car-Parrinello molecular dynamics simulation of Fe 3+ (aq)

The optimized geometry and energetic properties of Fe(D2O)n 3+ clusters, with n = 4 and 6, have been studied with density-functional theory calculations and the BLYP functional, and the hydration of a single Fe 3+ ion in a periodic box with 32 water molecules at room temperature has been studied with Car-Parrinello molecular dynamics and the same functional. We have compared the results from the CPMD simulation with classical MD simulations, using a flexible SPC-based water model and the same number of water molecules, to evaluate the relative strengths and weaknesses of the two MD methods. The classical MD simulations and the CPMD simulations both give Fe-water distances in good agreement with experiment, but for the intramolecular vibrations, the classical MD yields considerably better absolute frequencies and ion-induced frequency shifts. On the other hand, the CPMD method performs considerably better than the classical MD in describing the intramolecular geometry of the water molecule in the first hydration shell and the average first shell...second shell hydrogen-bond distance. Differences between the two methods are also found with respect to the second-shell water orientations. The effect of the small box size (32 vs 512 water molecules) was evaluated by comparing results from classical simulations using different box sizes; non-negligible effects are found for the ion-water distance and the tilt angles of the water molecules in the second hydration shell and for the O-D stretching vibrational frequencies of the water molecules in the first hydration shell.