Refinement of the primary hydration shell model for molecular dynamics simulations of large proteins

A realistic representation of water molecules is important in molecular dynamics simulation of proteins. However, the standard method of solvating biomolecules, that is, immersing them in a box of water with periodic boundary conditions, is computationally expensive. The primary hydration shell (PHS) method, developed more than a decade ago and implemented in CHARMM, uses only a thin shell of water around the system of interest, and so greatly reduces the computational cost of simulations. Applying the PHS method, especially to larger proteins, revealed that further optimization and a partial reworking was required and here we present several improvements to its performance. The model is applied to systems with different sizes, and both water and protein behaviors are compared with those observed in standard simulations with periodic boundary conditions and, in some cases, with experimental data. The advantages of the modified PHS method over its original implementation are clearly apparent when it is applied to simulating the 82 kDa protein Malate Synthase G. © 2009 Wiley Periodicals, Inc. J Comput Chem 2009     

Molecular dynamics simulations of fluid methane properties using ab initio intermolecular interaction potentials

Intermolecular interaction energy data for the methane dimer have been calculated at a spectroscopic accuracy and employed to construct an ab initio potential energy surface (PES) for molecular dynamics (MD) simulations of fluid methane properties. The full potential curves of the methane dimer at 12 symmetric conformations were calculated by the supermolecule counterpoise‐corrected second‐order Møller‐Plesset (MP2) perturbation theory. Single‐point coupled cluster with single and double and perturbative triple excitations [CCSD(T)] calculations were also carried out to calibrate the MP2 potentials. We employed Pople's medium size basis sets [up to 6‐311++G(3df, 3pd)] and Dunning's correlation consistent basis sets (cc‐pVXZ and aug‐cc‐pVXZ, X = D, T, Q). For each conformer, the intermolecular carbon–carbon separation was sampled in a step 0.1 Å for a range of 3–9 Å, resulting in a total of 732 configuration points calculated. The MP2 binding curves display significant anisotropy with respect to the relative orientations of the dimer. The potential curves at the complete basis set (CBS) limit were estimated using well‐established analytical extrapolation schemes. A 4‐site potential model with sites located at the hydrogen atoms was used to fit the ab initio potential data. This model stems from a hydrogen–hydrogen repulsion mechanism to explain the stability of the dimer structure. MD simulations using the ab initio PES show quantitative agreements on both the atom‐wise radial distribution functions and the self‐diffusion coefficients over a wide range of experimental conditions. © 2008 Wiley Periodicals, Inc. J Comput Chem 2009     

Extracting hydration sites around proteins from explicit water simulations

Two new methods are assessed for determining the location of hydration sites around proteins from computer simulation. Current methods extract hydration sites from peaks in the water density constructed in the protein frame. However, the dynamic nature of the water molecules, the nearby protein residues, and the protein reference frame as a whole tend to smear out the water density, making it more difficult to resolve sites. Two techniques are introduced to better resolve the water density. The first is to construct the water density from the time-averaged position of each water molecule in the protein frame while the water remains within a given distance of this averaged position. The second technique is to construct the water density from the time-averaged position of each water in the reference frame only of the nearby residues. Criteria for determining hydration sites from the water density are examined. Both techniques are found to significantly improve the detail in the water density and the number of hydration sites detected.     

聚乙烯在石墨表面结晶的分子动力学模拟

采用分子动力学方法模拟了聚乙烯在石墨(001)表面的吸附和结晶过程;直观的给出了聚乙烯链被石墨(001)面吸附并诱导形成有序的片层晶体的过程;发现结晶温度对得到的有序结构中的聚乙烯链相对石墨表面的特定取向有影响(300 K和600 K时的取向方向不同);表面覆盖率影响聚乙烯吸附层的厚度,对取向的方向无影响.     

Molecular dynamics simulation of thymine glycol-lesioned DNA reveals specific hydration at the lesion

One nanosecond molecular dynamics (MD) simulation of a thymine glycol (TG)-lesioned part of human lymphoblast AG9387 was performed to determine structural changes in DNA molecule caused by the presence of a lesion. These changes can be significant for proper recognition of lesions by a repair enzyme. Thymine glycol is the DNA oxidative lesion formed by addition of OH radicals to C5 and C6 atoms of the thymine base. This lesion is known as causing Cockayne Syndrome-inherited genetic disorder. Distribution of water molecules in a hydration shell around the DNA molecule was analyzed for its contribution to the recognition of the TG lesion by the repair enzyme. The results of MD simulation show there is a specific DNA structural configuration formed at the lesion. After 500 ps the DNA is bent in a kink at the TG site. This change dislocates the glycosyl bond at C5' to a position closer to the DNA surface, and thus its atoms are more exposed to the surrounding water shell. The increased number of water molecules that are close to the TG site indicates that the glycosyl bond may be easily contacted by the repair enzyme. In addition, the higher number of water molecules at the TG site substantiates the importance of water-mediated hydrogen bonds created between the repair enzyme and the DNA upon formation of the complex. Copyright 2001 John Wiley & Sons, Inc. J Comput Chem 22: 1723-1731, 2001     

Molecular dynamics simulations of triflic acid and triflate ion/water mixtures: a proton conducting electrolytic component in fuel cells

Triflic acid is a functional group of perflourosulfonated polymer electrolyte membranes where the sulfonate group is responsible for proton conduction. However, even at extremely low hydration, triflic acid exists as a triflate ion. In this work, we have developed a force-field for triflic acid and triflate ion by deriving force-field parameters using ab initio calculations and incorporated these parameters with the Optimized Potentials for Liquid Simulations - All Atom (OPLS-AA) force-field. We have employed classical molecular dynamics (MD) simulations with the developed force field to characterize structural and dynamical properties of triflic acid (270-450 K) and triflate ion/water mixtures (300 K). The radial distribution functions (RDFs) show the hydrophobic nature of CF(3) group and presence of strong hydrogen bonding in triflic acid and temperature has an insignificant effect. Results from our MD simulations show that the diffusion of triflic acid increases with temperature. The RDFs from triflate ion/water mixtures shows that increasing hydration causes water molecules to orient around the SO(3)(-) group of triflate ions, solvate the hydronium ions, and other water molecules. The diffusion of triflate ions, hydronium ion, and water molecules shows an increase with hydration. At λ = 1, the diffusion of triflate ion is 30 times lower than the diffusion of triflic acid due to the formation of stable triflate ion-hydronium ion complex. With increasing hydration, water molecules break the stability of triflate ion-hydronium ion complex leading to enhanced diffusion. The RDFs and diffusion coefficients of triflate ions, hydronium ions and water molecules resemble qualitatively the previous findings using per-fluorosulfonated membranes.     

A hydration study of (1-->4) and (1-->6) linked alpha-glucans by comparative 10 ns molecular dynamics simulations and 500-MHz NMR

The hydration behavior of two model disaccharides, methyl-alpha-D-maltoside (1) and methyl-alpha-D-isomaltoside (2), has been investigated by a comparative 10 ns molecular dynamics study. The detailed hydration of the two disaccharides was described using three force fields especially developed for modeling of carbohydrates in explicit solvent. To validate the theoretical results the two compounds were synthesized and subjected to 500 MHz NMR spectroscopy, including pulsed field gradient diffusion measurements (1: 4.0. 10(-6) cm(2). s(-1); 2: 4.2. 10(-6) cm(2). s(-1)). In short, the older CHARMM-based force field exhibited a more structured carbohydrate-water interaction leading to better agreement with the diffusional properties of the two compounds, whereas especially the alpha-(1-->6) linkage and the primary hydroxyl groups were inaccurately modeled. In contrast, the new generation of the CHARMM-based force field (CSFF) and the most recent version of the AMBER-based force field (GLYCAM-2000a) exhibited less structured carbohydrate-water interactions with the result that the diffusional properties of the two disaccharides were underestimated, whereas the simulations of the alpha-(1-->6) linkage and the primary hydroxyl groups were significantly improved and in excellent agreement with homo- and heteronuclear coupling constants. The difference between the two classes of force field (more structured and less structured carbohydrate-water interaction) was underlined by calculation of the isotropic hydration as calculated by radial pair distributions. At one extreme, the radial O em leader O pair distribution function yielded a peak density of 2.3 times the bulk density in the first hydration shell when using the older CHARMM force field, whereas the maximum density observed in the GLYCAM force field was calculated to be 1.0, at the other extreme.     

Molecular dynamics simulations of the hydration of poly(vinyl methyl ether):Hydrogen bonds and quasi-hydrogen bonds

Atomistic detailed hydration structures of poly(vinyl methyl ether)(PVME) have been investigated by molecular dynamics simulations under 300 K at various concentrations. Both radial distribution functions and the distance distributions between donors and acceptors in hydrogen bonds show that the hydrogen bonds between the polymer and water are shorter by 0.005 nm than those between water molecules. The Quasi-hydrogen bonds take only 7.2% of the van der Waals interaction pairs. It was found the hydrogen bonds are not evenly distributed along the polymer chain,and there still exists a significant amount(10%) of ether oxygen atoms that are not hydrogen bonded to water at a concentration as low as 3.3%. This shows that in polymer solutions close contacts occur not only between polymer chains but also between chain segments within the polymer,which leads to inefficient contacts between ether oxygen atoms and water molecules. Variation of the quasi-hydrogen bonds with the concentration is similar to that of hydrogen bonds,but the ratio of the repeat units forming quasi-hydrogen bonds to those forming hydrogen bonds approaches 0.2. A transition was found in the demixing enthalpy at around 30% measured by dynamic testing differential scanning calorimetry(DTDSC) for aqueous solutions of a mono-dispersed low molecular weight PVME,which can be related to the transition of the fractions of hydrogen bonds and quasi-hydrogen bonds at ～27%. The transition of the fractions of hydrogen bonds and quasi-hydrogen bonds at ～27% can be used to explain the demixing enthalpy transition at 30% at a molecular scale. In addition,at the concentration of 86%,each ether oxygen atom bonded with water is assigned 1.56 water molecules on average,and 'free' water molecules emerge at the concentration of around 54%.     

Molecular dynamics simulations of osmosis and reverse osmosis in solutions

Computer simulation studies using the method of molecular dynamics have been carried out to investigate osmosis and reverse osmosis in solutions separated by semi-permeable membranes. The method has been used to study the dynamic approach to equilibrium in such systems from their initial nonequilibrium state. In addition density profiles of both the solute and solvent molecules have been investigated, especially near the walls for adsorption effects. Finally the diffusion coefficients and osmotic pressure have also been measured.Our results show both osmosis and reverse osmosis, as well as a smooth transition between the two when either the solution concentration is changed, or the density (pressure) difference between the solvent and solution compartments is varied. We believe this new method can be used to improve our understanding of these two important phenomena at the molecular level.     

Liquid properties of dimethyl ether from molecular dynamics simulations using Ab Initio force fields

We have used molecular dynamic simulations to study the structural and dynamical properties of liquid dimethyl ether (DME) with a newly constructed ab initio force field in this article. The ab initio potential energy data were calculated at the second order Møller‐Plesset (MP2) perturbation theory with Dunning's correlation consistent basis sets (up to aug‐cc‐pVQZ). We considered 17 configurations of the DME dime for the orientation sampling. The calculated MP2 potential data were used to construct a 3‐site united atom force field model. The simulation results are compared with those using the empirical force field of Jorgensen and Ibrahim (Jorgensen and Ibrahim, J Am Chem Soc 1981, 103, 3976) and with available experimental measurements. We obtain quantitative agreements for the atom‐wise radial distribution functions, the self‐diffusion coefficients, and the shear viscosities over a wide range of experimental conditions. This force field thus provides a suitable starting point to predict liquid properties of DME from first principles intermolecular interactions with no empirical data input a priori. © 2012 Wiley Periodicals, Inc.     

Molecular dynamics simulations of the mechanisms of thermal conduction in methane hydrates

The thermal conductivity of methane hydrate is an important physical parameter affecting the processes of methane hydrate exploration,mining,gas hydrate storage and transportation as well as other applications.Equilibrium molecular dynamics simulations and the Green-Kubo method have been employed for systems from fully occupied to vacant occupied sI methane hydrate in order to estimate their thermal conductivity.The estimations were carried out at temperatures from 203.15 to 263.15 K and at pressures from 3 to 100 MPa.Potential models selected for water were TIP4P,TIP4P-Ew,TIP4P/2005,TIP4P-FQ and TIP4P/Ice.The effects of varying the ratio of the host and guest molecules and the external thermobaric conditions on the thermal conductivity of methane hydrate were studied.The results indicated that the thermal conductivity of methane hydrate is essentially determined by the cage framework which constitutes the hydrate lattice and the cage framework has only slightly higher thermal conductivity in the presence of the guest molecules.Inclusion of more guest molecules in the cage improves the thermal conductivity of methane hydrate.It is also revealed that the thermal conductivity of the sI hydrate shows a similar variation with temperature.Pressure also has an effect on the thermal conductivity,particularly at higher pressures.As the pressure increases,slightly higher thermal conductivities result.Changes in density have little impact on the thermal conductivity of methane hydrate.     

Molecular dynamics simulations of atomistic hydration structures of poly(vinyl methyl ether)

Molecular dynamics simulations have been performed on the aqueous solutions of poly(vinyl methyl ether) (PVME) at various concentrations. Both radial and spatial distribution functions are used to investigate the detailed hydration structures. The structures of water are found to get increasingly concentrated when polymers are introduced and the water motions are severely hindered by the polymer matrix. At low concentrations, larger populations of tt conformers in meso dyads than those at higher concentrationsare found and this phenomenon is believed to be due to the increasing in bonding of water molecule to two ether oxygens in meso dyad. At higher concentrations, the size and conformations of polymers are quite similar to those in bulk. A transition of hydrogen bond fractions between PVME and water at around the concentration of 0.3 is observed and this value is perfectly in agreement with the results of conformational analysis and Raman spectra. Second neighbor hydrogen bond statistics revealed the domination of complicated hydrogen bond networks at low concentrations, but single hydrogen bonds as well as isolated clusters composed of 2-4 water molecules are usual around each polymer repeat unit.     

Conductance of Cs+ ion in water: Molecular dynamics simulation

The continuum theory of Hubbard-Onsager predicts for the friction coefficients the following behavior: >0 and /P<0. In contrast to Hubbard-Onsager theory, experimental observations on Cs⁺ ion in water show that at low temperatures <0 and /P>0. To explain the observed behavior of Nakahara et al. proposed the passage through cavities (PTC) mechanism. We performed a molecular dynamics computer simulation to determine if the PTC mechanism is responsible for the observed behavior of . No passage through cavities was observed. Molecular dynamics computer simulations were performed on Cs⁺ ion in water at temperature of 268 K and densities of 1.00 and 1.083 g-cm^–3. Our results indicate that the observed behavior of for Cs⁺ ion is related to the difference in the reorientation times of water molecules in the solvation shell and in the bulk.     

Molecular dynamics simulations of a poly(p-phenylene) oligomer

We have studied the conformation and coefficient of thermal expansion in the poly(p-phenylene) oligomer p-sexiphenyl (C₃₆H₂₆) by molecular dynamics simulations. Studies of the backbone phenyl–phenyl torsion angle in a simulated p-sexiphenyl crystal at room temperature indicate the presence of torsional librations of approximately ±20°. Further analysis of the phenyl–phenyl backbone torsion angle in less closely packed regions of the simulated crystal (crystal ends) indicate the presence of 180° phenyl ring flips, in agreement with solid-state deuterium NMR data on poly(p-phenylene oligomers). The linear coefficient of thermal expansion was also calculated and found to be negative, in qualitative agreement with experimental data on rigid-rod compounds. © 1993 John Wiley & Sons, Inc.     

Molecular dynamics simulations of LiCl association and NaCl association in water by means of ABEEM/MM

Constrained molecular dynamics simulations have been used to investigate the LiCl and NaCl ionic association in water in terms of atom-bond electronegativity equalization method fused into molecular mechanics (ABEEM/MM). The simulations make use of the seven-site fluctuating charge and flexible ABEEM-7P water model, based on which an ion-water interaction potential has been constructed. The mean force and the potential of mean force for LiCl and NaCl in water, the charge distributions, as well as the structural and dynamical properties of contact ion pair dissociation have been investigated. The results are reasonable and informative. For LiCl ion pair in water, the solvent-separated ion pair configurations are more stable than contact ion pair configurations. The calculated PMF for NaCl in water indicates that contact ion pair and solvent-separated ion pair configurations are of comparable stability.     

Quantum mechanical/molecular mechanical simulations of the Tl(III) ion in water

Classical molecular dynamics (MD) and combined quantum mechanical/molecular mechanical (QM/MM) MD simulations have been performed to investigate the structural and dynamical properties of the Tl(III) ion in water. A six-coordinate hydration structure with a maximum probability of the Tl-O distance at 2.21 A was observed, which is in good agreement with X-ray data. The librational and vibrational spectra of water molecules in the first hydration shell are blue-shifted compared with those of pure liquid water, and the Tl-O stretching force constant was evaluated as 148 Nm(-1). Both structural and dynamical properties show a distortion of the first solvation shell structure. The second shell ligands' mean residence time was determined as 12.8 ps. The Tl(III) ion can be classified as "structure forming" ion; the calculated hydration energy of -986 +/- 9 kcal mol agrees well with the experimental value of -986 kcal mol.     

Rapid approximate calculation of water binding free energies in the whole hydration domain of (bio)macromolecules

The evaluation of water binding free energies around solute molecules is important for the thermodynamic characterization of hydration or association processes. Here, a rapid approximate method to estimate water binding free energies around (bio)macromolecules from a single molecular dynamics simulation is presented. The basic idea is that endpoint free‐energy calculation methods are applied and the endpoint quantities are monitored on a three‐dimensional grid around the solute. Thus, a gridded map of water binding free energies around the solute is obtained, that is, from a single short simulation, a map of favorable and unfavorable water binding sites can be constructed. Among the employed free‐energy calculation methods, approaches involving endpoint information pertaining to actual thermodynamic integration calculations or endpoint information as exploited in the linear interaction energy method were examined. The accuracy of the approximate approaches was evaluated on the hydration of a cage‐like molecule representing either a nonpolar, polar, or charged water binding site and on α‐ and β‐cyclodextrin molecules. Among the tested approaches, the linear interaction energy method is considered the most viable approach. Applying the linear interaction energy method on the grid around the solute, a semi‐quantitative thermodynamic characterization of hydration around the whole solute is obtained. Disadvantages are the approximate nature of the method and a limited flexibility of the solute. © 2016 Wiley Periodicals, Inc.     

Molecular dynamics of electrosprayed water nanodroplets containing sodium bis(2‐ethylhexyl)sulfosuccinate

The behavior of aqueous solutions of sodium bis(2‐ethylhexyl)sulfosuccinate (AOTNa) subject to electrospray ionization (ESI) has been investigated by molecular dynamics (MD) simulations at three temperatures (350, 500 and 800 K). We consider several types of water nanodroplets containing AOTNa molecules and composed of a fixed number of water molecules (1000), N⁰_AOT AOT^? anions (N⁰_AOT = 0, 5, 10) and N⁰_Na sodium ions (N⁰_Na = 0, 5, 10, 15, 20): in a short time scale (less than 1 ns), the AOTNa molecules, initially forming direct micelles in the interior of the water nanodroplets, are observed in all cases to diffuse nearby the nanodroplet surface, so that the hydrophilic heads and sodium ions become surrounded by water molecules, whereas the alkyl chains lay at the droplet surface. Meanwhile, evaporation of water molecules and of solvated sodium ions occurs, leading to a decrease of the droplet size and charge. At 350 K, no ejection of neutral or charged surfactant molecules is observed, whereas at 500 K, some fragmentation occurs, and at 800 K, this event becomes more frequent. The interplay of all these processes, which depend on the values of temperature, N⁰_AOT and N⁰_Na eventually leads to anhydrous charged surfactant aggregates with prevalence of monocharged ones, in agreement with experimental results of ESI mass spectrometry. The quantitative analysis of the MD trajectories allows to evidence molecular details potentially useful in designing future ESI experimental conditions. Copyright © 2013 John Wiley & Sons, Ltd.