Periodic Box Of Water Research Articles

Microscopic Mechanism of Proton Transfer in Pure Water under Ambient Conditions.

Water molecules and the associated proton transfer (PT) are prevalent in chemical and biological systems and have been a hot research topic. Spectroscopic characterization and ab initio molecular dynamics (AIMD) simulations have previously revealed insights into acidic and basic liquids. Presumably, the situation in the acidic/basic solution is not necessarily the same as in pure water; in addition, the autoionization constant for water is only 10-14 under ambient conditions, making the study of PT in pure water challenging. To overcome this issue, we modeled periodic water box systems containing 1000 molecules for tens of nanoseconds based on a neural network potential (NNP) with quantum mechanical accuracy. The NNP was generated by training a dataset containing the energies and atomic forces of 17 075 configurations of periodic water box systems, and these data points were calculated at the MP2 level that considers electron correlation effects. We found that the size of the system and the duration of the simulation have a significant impact on the convergence of the results. With these factors considered, our simulations showed that hydronium (H3O+) and hydroxide (OH-) ions in water have distinct hydration structures, thermodynamic and kinetic properties, e.g., the longer-lasting and more stable hydrated structure of OH- ions than that of H3O+, as well as a significantly higher free energy barrier for the OH--associated PT than that of H3O+, leading the two to exhibit completely different PT behaviors. Given these characteristics, we further found that PT via OH- ions tends not to occur multiple times or between many molecules. In contrast, PT via H3O+ can synergistically occur among multiple molecules and prefers to adopt a cyclic pattern among three water molecules, while it occurs mostly in a chain pattern when more water molecules are involved. Therefore, our studies provide a detailed and solid microscopic explanation for the PT process in pure water.

Molecular modelling and quantum mechanical calculations of the hydration free energy of buckminsterfullerene

The hydration free energy ΔGs of fullerene C60 has been determined using the quantum mechanical polarizable continuum model, microscopic free energy perturbation molecular dynamics (FEP/MD) techniques and compilation of known experimental data. The standard model density and united atom model were employed in quantum calculations. The FEP/MD values of ΔGs(C60) were calculated from modelling solute insertion in periodic water box and solute annihilation in a water sphere that additionally employed estimation of the cavitation term. The computed hydration free energy varies by up to 80kJmol−1 depending on model parameters and acquires positive values when using parameterizations for fullerene.

Dynamics and mechanisms of interactions between ring-shaped heterohexameric TIP49a/b protein complexes and double-stranded DNA

TIP49a and TIP49b, highly conserved proteins belonging to the AAA+ superfamily of DNA-dependent ATPases, participate in numerous cell processes, such as chromatin remodeling, regulation of gene transcription and mitotic cell division, maintenance of genome stability, and snoRNP biogenesis, as well as in the formation of active DNA–telomerase complexes. It has been shown that they are involved in complex networks of protein–protein interactions and, in spite of their structural similarity, sometimes perform opposite functions. Although these proteins exhibit a wide range of activities, the mechanisms of their actions are still poorly understood. In this work, ring-shaped heterohexameric TIP49a/b complexes containing in their central channel short fragments of double-stranded DNA (dsDNA, 20 base pairs of different GC composition) were obtained for the first time using the molecular docking technique, while methods of molecular dynamics in a periodic water box were applied to investigate the conformational dynamics of these proteins and the mechanisms of their helicase activity. It was found that (1) interaction of a DNA helix with positively charged protein loops inside the central channel of the ring-shaped hexameric complex caused unwinding of the helix; (2) the unwinding occurred only inside the hexameric ring, whereas the tails of the helix, which lie outside, retained the initial classical B-form conformation throughout the 50 ns of molecular dynamics; and (3) the presence of ATP in TIP49a/b complexes affected the dynamics and the final structure of dsDNA, causing partial breaking of complementary bonds in GC-poor DNA sequences.

Correlated motion of protein subdomains and large-scale conformational flexibility of RecA protein filament

Based on X-ray crystallographic data available at Protein Data Bank, we have built molecular dynamics (MD) models of homologous recombinases RecA from E. coli and D. radiodurans. Functional form of RecA enzyme, which is known to be a long helical filament, was approximated by a trimer, simulated in periodic water box. The MD trajectories were analyzed in terms of large-scale conformational motions that could be detectable by neutron and X-ray scattering techniques. The analysis revealed that large-scale RecA monomer dynamics can be described in terms of relative motions of 7 subdomains. Motion of C-terminal domain was the major contributor to the overall dynamics of protein. Principal component analysis (PCA) of the MD trajectories in the atom coordinate space showed that rotation of C-domain is correlated with the conformational changes in the central domain and N-terminal domain, that forms the monomer-monomer interface. Thus, even though C-terminal domain is relatively far from the interface, its orientation is correlated with large-scale filament conformation. PCA of the trajectories in the main chain dihedral angle coordinate space implicates a co-existence of a several different large-scale conformations of the modeled trimer. In order to clarify the relationship of independent domain orientation with large-scale filament conformation, we have performed analysis of independent domain motion and its implications on the filament geometry.

Molecular dynamics simulations of the interaction of glucose with imidazole in aqueous solution

Molecular dynamics simulations were carried out on a concentrated system of β-d-glucopyranose and imidazole molecules in a periodic box of water at 298K. The purpose of the simulations was to determine whether or not there was any tendency for these two solutes to associate in an aqueous environment, as has previously been observed for other planar functional groups from amino acid side chains, such as the indole group of tryptophan or the phenolic group of tyrosine. A weak stacking interaction between β-d-glucopyranose, as a model for cellulose, and imidazole was indeed observed, with an energy of ∼0.25kcal/mol per pair, less than kT. Somewhat surprisingly, considerable imidazole self-association into small aggregates (dimers and trimers) was also observed, with binding energies of ∼0.4kcal/mol per pair, although still less than kT. Considerable non-stacked interactions between glucose and imidazole through hydrogen bonding were also found. These hydrogen bonds primarily involved the N3 atom of imidazole, because the N–H group of N1 was insufficiently polar to compete for water hydrogen bond partners.

Molecular Dynamics Simulation Studies of Caffeine Aggregation in Aqueous Solution

Molecular dynamics simulations were carried out on a system of eight independent caffeine molecules in a periodic box of water at 300 K, representing a solution near the solubility limit for caffeine at room temperature, using a newly developed CHARMM-type force field for caffeine in water. Simulations were also conducted for single caffeine molecules in water using two different water models (TIP3P and TIP4P). Water was found to structure in a complex fashion around the planar caffeine molecules, which was not sensitive to the water model used. As expected, extensive aggregation of the caffeine molecules was observed, with the molecules stacking their flat faces against one another like coins, with their methylene groups staggered to avoid steric clashes. A dynamic equilibrum was observed between large n-mers, including stacks with all eight solute molecules, and smaller clusters, with the calculated osmotic coefficient being in acceptable agreement with the experimental value. The insensitivity of the results to water model and the congruence with experimental thermodynamic data suggest that the observed stacking interactions are a realistic representation of the actual association mechanism in aqueous caffeine solutions.

Molecular dynamics of sialic acid analogues and their interaction with influenza hemagglutinin

Synthetic sialic acid analogues with multiple modifications at different positions(C-1/C-2/C-4/C-8/C-9) are investigated by molecular mechanics and molecular dynamics to determine their conformational preferences and structural stability to interact with their natural receptors. Sialic acids with multiple modifications are soaked in a periodic box of water as solvent. Molecular mechanics and a 2 nanosecond molecular dynamics are done using amber force fields with 30 picosecond equilibrium. Direct and water mediated hydrogen bonds existing in the sialic acid analogues, aiding for their structural stabilization are identified in this study. The accessible conformations of side chain linkages of sialic acid analogues holding multiple substituents are determined from molecular dynamics trajectory at every 1ps interval. Transitions between different minimum energy regions in conformational maps are also noticed in C-1, C-2, C-4, C-8 and C-9 substituents. Docking studies were done to find the binding mode of the sialic acid analogues with Influenza hemagglutinin. This finding provides stereo chemical explanation and conformational preference of sialic acid analogues which may be crucial for the design of sialic acid analogues as inhibitors for different sialic acid specific pathogenic proteins such as influenza toxins and neuraminidases.

Towards temperature-dependent coarse-grained potentials of side-chain interactions for protein folding simulations. I: Molecular dynamics study of a pair of methane molecules in water at various temperatures

By means of molecular dynamics simulations of a pair of methane molecules in a TIP3P periodic water box with the NVT scheme at six temperatures and, additionally, the NPT scheme at three temperatures ranging from T = 283 to 373 K, we determined the potential of mean force (PMF) of pairs of interacting methane molecules in water as functions of distance between the methane molecules. The PMFs converge to a single baseline only for r> 11 A at all temperatures. The curves of the dimensionless PMF obtained at different temperatures with the NVT scheme overlap almost perfectly in the region of the contact minimum and still very well in the regions of the desolvation maximum and the solvent-separated minimum, which suggests that the temperature-dependent hydrophobic interaction potentials at constant volume in united-residue force fields can be obtained by scaling the respective dimensionless potentials by RT, R being the universal gas constant. For the dimensionless potentials of mean force obtained with the NPT scheme, the depth of the contact minimum increases, whereas the height of the desolvation maximum and the depth of the solvent-separated minimum decrease with temperature, in agreement with results reported in the literature.

How hydrophobic buckminsterfullerene affects surrounding water structure.

The hydrophobic hydration of fullerenes in water is of significant interest as the most common Buckminsterfullerene (C60) is a mesoscale sphere; C60 also has potential in pharmaceutical and nanomaterial applications. We use an all-atom molecular dynamics simulation lasting hundreds of nanoseconds to determine the behavior of a single molecule of C60 in a periodic box of water, and compare this to methane. A C60 molecule does not induce drying at the surface; however, unlike a hard sphere methane, a hard sphere C60 solute does. This is due to a larger number of attractive Lennard-Jones interactions between the carbon atom centers in C60 and the surrounding waters. In these simulations, water is not uniformly arranged but rather adopts a range of orientations in the first hydration shell despite the spherical symmetry of both solutes. There is a clear effect of solute size on the orientation of the first hydration shell waters. There is a large increase in hydrogen-bonding contacts between waters in the C60 first hydration shell. There is also a disruption of hydrogen bonds between waters in the first and second hydration shells. Water molecules in the first hydration shell preferentially create triangular structures that minimize the net water dipole near the surface near both the methane and C60 surface, reducing the total energy of the system. Additionally, in the first and second hydration shells, the water dipoles are ordered to a distance of 8 A from the solute surface. We conclude that, with a diameter of approximately 1 nm, C60 behaves as a large hydrophobic solute.

Mapping the Binding Site on Small Ankyrin 1 for Obscurin

Small ankyrin 1 (sAnk1), an integral protein of the sarcoplasmic reticulum encoded by the ANK1 gene, binds with nanomolar affinity to the C terminus of obscurin, a giant protein surrounding the contractile apparatus in striated muscle. We used site-directed mutagenesis to characterize the binding site on sAnk1, specifically addressing the role of two putative amphipathic, positively charged helices. We measured binding qualitatively by blot overlay assays and quantitatively by surface plasmon resonance and showed that both positively charged sequences are required for activity. We showed further that substitution of a lysine or arginine with an alanine or glutamate located at the same position along either of the two putative helices has similar inhibitory or stimulatory effects on binding and that the effects of a particular mutation depended on the position of the mutated amino acid in each helix. We modeled the structure of the binding region of sAnk1 by homology with ankyrin repeats of human Notch1, which have a similar pattern of charged and hydrophobic residues. Our modeling suggested that each of the two positively charged sequences forms pairs of amphipathic, anti-parallel alpha-helices flanked by beta-hairpin-like turns. Most of the residues in homologous positions along each helical unit have similar, though not identical, orientations. CD spectroscopy confirmed the alpha-helical content of sAnk1, approximately 33%, predicted by the model. Thus, structural and mutational studies of the binding region on sAnk1 for obscurin suggest that it consists of two ankyrin repeats with very similar structures.

Towards accurate solvation dynamics of divalent cations in water using the polarizable amoeba force field: From energetics to structure

Molecular dynamics simulations were performed using a modified amoeba force field to determine hydration and dynamical properties of the divalent cations Ca2+ and Mg2+. The extension of amoeba to divalent cations required the introduction of a cation specific parametrization. To accomplish this, the Thole polarization damping model parametrization was modified based on the ab initio polarization energy computed by a constrained space orbital variation energy decomposition scheme. Excellent agreement has been found with condensed phase experimental results using parameters derived from gas phase ab initio calculations. Additionally, we have observed that the coordination of the calcium cation is influenced by the size of the periodic water box, a recurrent issue in first principles molecular dynamics studies.

Comparison of two approaches to potential of mean force calculations of hydrophobic association: particle insertion and weighted histogram analysis methods

We investigated the convergence behaviour of potential of mean force (PMF) calculations using the particle insertion method (PI) and the umbrella sampling with weighted histogram analysis method (WHAM). For PI, two approaches were tested: insertion of two particles at various separations into pure water, and insertion of a second particle into a heterogeneous system consisting of a first particle and water. The comparison of these methods is illustrated by a study of the PMF between two neon atoms and two methane molecules embedded in a periodic box of water. Two water models, TIP3P and TIP4P, were used in the two series of simulations. The choice of the method employed to estimate the PMF affects the efficiency of the calculation and, in turn, the amount of sampling necessary to attain proper convergence. For small particles, with radii similar to that of the neon atom, PI converges faster than WHAM calculations. For larger particles, umbrella sampling with WHAM converges faster, and the size of the methane molecule seems to be close to the practical limit of convergence of the PI method. The umbrella sampling/WHAM approach cannot provide the absolute value of the PMF, and depends on an estimation of the baseline at large distances; however, the PMF curves for two neon atoms in water, for which both the particle-insertion and the umbrella sampling/WHAM approaches give results of comparable quality, converge at large distances, which justifies the assumption that the two- and multi-body contribution to the PMF of clusters of hydrophobic solutes tends to zero at distances beyond the solvent-separated minimum.

Simulations of Structure and Vibrational Spectra of Deoxyoctanucleotides

Combined molecular dynamics and ab initio computations were applied for analysis of infrared absorption and vibrational circular dichroism spectra of deoxyoctanucleotides. Unlike for previous idealized models, molecular geometries of these shorter DNA fragments in solutions were obtained as dynamic averages from simulations in a periodic water box. Vibrational spectra for the whole octamers including hydrogen-bonded solvent molecules were simulated on the basis of density-functional calculations on small fragments and subsequent transfer of molecular property tensors. Explicit and continuum solvent models were compared. Apparently, the DNA segments retain an approximate B-conformation in the aqueous solutions, but the terminal base pairs significantly deviate from the planar arrangement and the vibrational circular dichroism spectrum for (CG)(2) nucleotide indicates a larger average helical twist. Sodium counterions moved freely around the molecule during the simulation and do not influence spectral intensities. Simulated absorption spectra faithfully reproduced the experimental signal of principal functional groups, while only qualitative agreement was obtained for the dependence on the basis sequence.

Molecular dynamics study of the influence of calcium ions on the conformation of gelsolin S2 domain

Gelsolin is an actin-severing protein whose action is promoted by Ca 2+ ions and inhibited by binding to lipid phosphoinositides incorporated in the inner leaflet of the plasma membrane inner lipid bilayer. In this study, we carried out molecular dynamics (MD) simulations to investigate the influence of calcium cations on the conformation of gelsolin S2 domain. First, gelsolin S2 domain taken from the crystal structure of apo-gelsolin (PDB code: 1D0N) was subjected to three 1100 ps MD simulations in a periodic water box with the amber 5.0 force field at T=298 K. In the first simulation (S2_Ca 2+) excess concentration of Ca 2+ was applied, in the second one (S2_phys) the concentration of Ca 2+ ions was physiological and in the third one (S2_w) no Ca 2+ ions were added. The results of MD simulations showed high conformational flexibility of the N-terminal part of the S2 domain. S2_w deviated from the starting structure considerably more that S2_phys and S2_Ca 2+ suggesting that Ca 2+ ions stabilize the conformation of the S2 domain of gelsolin.

Computer simulation of the hydrophobic hydration of concave surfaces

Hemispherical objects of radius from 6.5 to 12.2 Å were constructed from tessellated icosahedra where the vertices were hydrophobic ‘atoms’. These hemispheres were hydrated in a periodic box of water and molecular dynamics simulations in the NVE ensemble were performed on this solution at 320K. Results show that water molecules in the concave region exhibit a layered structure and have their dipoles preferentially oriented approximately orthogonal to the radial direction of the hemisphere. They possess fewer hydrogen bonds and near neighbours. Changes in the translational dynamics of these water molecules are not detected, but there is an acceleration of the re-orientational dynamics. This work thus represents the first systematic study of the hydration of concave hydrophobic surfaces.

Elastic molecular dynamics with self-consistent flexible constraints

A new algorithm for constrained molecular dynamics is proposed. In contrast to the standard approach, the constrained bond-length/bond-angle value is adjusted at each time step so that total energy is minimized with respect to the constrained distances. This can be viewed as modifying the equilibrium bond-length/bond-angle according to external and centripetal forces. Two approaches are constructed to implement the algorithm. Method I includes all energy terms, but it is neither holonomic nor symplectic. Method II neglects a rotational kinetic energy term, resulting in a more expensive symplectic integrator. Both integrators are reversible and well conserve total energy. Due to the iterative nature of the optimization, which requires an energy evaluation at each iteration, these methods are intended for use in conjunction with other integration techniques that require an optimization at each time step iteration, such as the self-consistent electronic polarization methods or iterative combined quantum mechanics/molecular mechanics (QM/MM) hybrid simulations. Simulation results for collisions between two diatomic molecules, two water molecules, and for a periodic water box are compared and contrasted with fully constrained and free dynamics. These results demonstrate the utility and efficacy of these new methods.

Alternative approaches to potential of mean force calculations: Free energy perturbation versus thermodynamic integration. Case study of some representative nonpolar interactions

We investigated the convergence behavior of potential of mean force (PMF) calculations using free energy perturbation (FEP), thermodynamic integration (TI), and “slow growth” (SG) techniques. The critical comparison of these alternative approaches is illustrated by the study of three different systems: two tagged argon atoms in a periodic box of argon, two methane molecules, and two benzene molecules maintained in a “T-shaped” conformation, both dimers embedded in a periodic box of water. The complete PMF simulations were carried out considering several protocols, in which the number of intermediate “λ” states, together with the amount of sampling per individual state, were varied. In most cases, as much as 1 ns of molecular dynamics (MD) sampling was used to derive each free energy profile. For the different systems examined, we find that FEP and TI unquestionably constitute robust computational methods leading to results of comparable accuracy. We also show that proper convergence of the free energy calculations, and further quantitative interpretation of the PMFs, requires total simulation times much higher than has been hitherto estimated. In some circumstances, the free energy profiles derived from FEP calculations tend to be slightly poorer than those obtained with TI, as a probable consequence of the greater sensitivity of FEP to the window spacing δλ. In the context of TI, and to a lesser extent FEP, simulations, it appears preferable to employ a limited number of “λ” points of the integrand involving extensive sampling, rather than numerous points with fewer samplings. Finally, we note that, at least in the case of nonpolar interactions, PMFs of reasonable quality can be generated using SG, and at a substantially lower cost than with either FEP or TI. © 1996 by John Wiley & Sons, Inc.

Free energy derivatives: A new method for probing the convergence problem in free energy calculations

AbstractThe convergence behavior of free energy calculations has been explored in more detail than in any previously reported work, using a model system of two neon atoms in a periodic box of water. We find that for thermodynamic integration‐type free energy calculations as much as a nanosecond or more molecular dynamics sampling is required to obtain a fully converged value for a single λ point of the integrand. The concept of “free energy derivatives” with respect to the individual parameters of the force field is introduced. This formalism allows the total convergence of the simulation to be deconvoluted into components. A determination of the statistical “sampling ratio” from these simulations indicates that for window‐type free energy calculations carried out in a periodic waterbox of typical size at least 0.6 ps of sampling should be performed at each window (0.7 ps if constraint contributions to the free energy are being determined). General methods to estimate and reduce the error in thermodynamic integration and free energy perturbation calculations are discussed. We show that the difficulty in applying such methods is determining a reliable estimate of the correlation length from a short series of data. © 1994 by John Wiley & Sons, Inc.

A molecular dynamics study of solvent behavior around a protein.

The solvent structure and behavior around a protein were examined by analyzing a trajectory of molecular dynamics simulation of the trp-holorepressor in a periodic box of water. The calculated self diffusion coefficient indicated that the solvent within 10 A of the protein had lower mobility. Examination of the solvent diffusion around different atoms of different kinds of residues showed no general tendency. This fact suggested that the solvent mobility is not influenced significantly by the kind of the atom or residue they solvated. Distribution analysis around the protein revealed two peaks of water oxygen: a sharp one at 2.8 A around polar and charged atoms and a broad one at approximately 3.4 A around apolar atoms. The former was stabilized by water-protein hydrogen bonds, and the latter was stabilized by water-water hydrogen bonds, suggesting the existence of a hydrophobic shell. An analysis of protein atom-water radial distribution functions confirmed these shell structures around polar or charged atoms and apolar ones.

Determining the contributions of constraints in free energy calculations: Development, characterization, and recommendations

We develop and/or characterize three methods for determining the contributions of constraints in free energy calculations. A new method for determining such contributions in thermodynamic integration (TI) calculations, the potential forces (PF) method, is developed and compared to a second method, the constraint forces (CF) method. Both methods are also compared to a previously described technique for calculating such contributions in free energy perturbation (FEP) simulations. We find that the TI/PF protocol is considerably more efficient than the TI/CF method, and is preferred except in cases where the constraints contributing to the free energy are part of a closed ring. Compared to TI/PF, the FEP method is shown to be relatively poor for generating potential of mean force (PMF) curves, though the FEP method is suitable for determining the ‘‘PMF bond contribution’’ in compositional free energy changes. PMF curves for a system of two neon atoms in a periodic box of water have been derived. The convergence behavior of the free energy derivative ∂G/∂R, where R is a distance constraint, has been examined in detail for this system. As much as a nanosecond of molecular dynamics sampling can be required to derive a fully converged value for this derivative at a single (λ) point. We have determined the sampling ratio for ∂G/∂R as a function of Ne–Ne distance for a 295 water (21 Å/side) periodic water box, and conclude that for free energy changes where long-distance constraint contributions are being determined, a minimum of 0.7 ps of sampling should be performed per window. When constraint contributions arise from short distances (such as when the PMF bond contribution is being calculated), correlation in the constraint derivative series dies out relatively quickly and the minimum sampling—then dictated by correlation in the nonbonded series—should be about 0.6 ps.