On the Ewald artifacts in computer simulations. The test-case of the octaalanine peptide with charged termini

The treatment of electrostatic interactions in molecular simulations is of fundamental importance. Ewald and related methods are being increasingly used to the detriment of cutoff schemes, which are known to produce several artifacts. A potential drawback of the Ewald method is the spatial periodicity that is imposed to the system, which could produce artifacts when applied in the simulation of liquids. In this work we analyze the octaalanine peptide with charged termini in explicit solvent, for which severe effects due to the use of Ewald sums were predicted using continuum electrostatics. Molecular Dynamics simulations for a total of 158 nanoseconds were performed in cells of different sizes. From the comparison of the results of different system sizes, no significant periodicity-induced artifacts were observed. It is argued that in current biomolecular simulations, the incomplete sampling is likely to affect the results to a larger extent than the artifacts induced by the use of Ewald sums.     

On the Ewald Artifacts in Computer Simulations. The Test-Case of the Octaalanine Peptide With Charged Termini

Abstract

The treatment of electrostatic interactions in molecular simulations is of fundamental importance. Ewald and related methods are being increasingly used to the detriment of cutoff schemes, which are known to produce several artifacts. A potential drawback of the Ewald method is the spatial periodicity that is imposed to the system, which could produce artifacts when applied in the simulation of liquids. In this work we analyze the octaalanine peptide with charged termini in explicit solvent, for which severe effects due to the use of Ewald sums were predicted using continuum electrostatics. Molecular Dynamics simulations for a total of 158 nanoseconds were performed in cells of different sizes. From the comparison of the results of different system sizes, no significant periodicity-induced artifacts were observed. It is argued that in current biomolecular simulations, the incomplete sampling is likely to affect the results to a larger extent than the artifacts induced by the use of Ewald sums.     

Roles of Boundary Conditions in DNA Simulations: Analysis of Ion Distributions with the Finite-Difference Poisson-Boltzmann Method

The wide use of lattice-sum strategies in biomolecular simulations has raised many questions on potential artifacts in these strategies. One interesting question is the artifacts in the counterion distributions of highly charged systems. As one would anticipate, Coulombic interactions under the periodic boundary condition may deviate noticeably from those under the free boundary condition in the highly charged systems, significantly influencing their counterion distributions. On the other hand, the electrostatic screening due to water molecules and mobile ions may effectively damp the possible periodic distortions in Coulombic interactions. Therefore, the magnitude of periodicity-induced artifacts in counterion distributions is not straightforward to dissect without detailed analyses. In this study, we have developed a hybrid explicit counterion/implicit salt representation of mobile ions to address this question. We have chosen a well-studied DNA for easy validation of the minimal hybrid ion representation. Our detailed analysis of continuum ion distributions, explicit ion distributions, radial counterion distribution functions, and sequence-dependent counterion distributions, however, indicates that periodicity artifacts are not apparent at the surface of the tested DNA. Nevertheless, influence of boundary conditions does show up starting at the second solvation shell and becomes apparent at the cell boundary.     

Speed of Conformational Change: Comparing Explicit and Implicit Solvent Molecular Dynamics Simulations

Adequate sampling of conformation space remains challenging in atomistic simulations, especially if the solvent is treated explicitly. Implicit-solvent simulations can speed up conformational sampling significantly. We compare the speed of conformational sampling between two commonly used methods of each class: the explicit-solvent particle mesh Ewald (PME) with TIP3P water model and a popular generalized Born (GB) implicit-solvent model, as implemented in the AMBER package. We systematically investigate small (dihedral angle flips in a protein), large (nucleosome tail collapse and DNA unwrapping), and mixed (folding of a miniprotein) conformational changes, with nominal simulation times ranging from nanoseconds to microseconds depending on system size. The speedups in conformational sampling for GB relative to PME simulations, are highly system- and problem-dependent. Where the simulation temperatures for PME and GB are the same, the corresponding speedups are approximately onefold (small conformational changes), between ∼1- and ∼100-fold (large changes), and approximately sevenfold (mixed case). The effects of temperature on speedup and free-energy landscapes, which may differ substantially between the solvent models, are discussed in detail for the case of miniprotein folding. In addition to speeding up conformational sampling, due to algorithmic differences, the implicit solvent model can be computationally faster for small systems or slower for large systems, depending on the number of solute and solvent atoms. For the conformational changes considered here, the combined speedups are approximately twofold, ∼1- to 60-fold, and ∼50-fold, respectively, in the low solvent viscosity regime afforded by the implicit solvent. For all the systems studied, 1) conformational sampling speedup increases as Langevin collision frequency (effective viscosity) decreases; and 2) conformational sampling speedup is mainly due to reduction in solvent viscosity rather than possible differences in free-energy landscapes between the solvent models.     

Analysis and evaluation of channel models: simulations of alamethicin

Alamethicin is an antimicrobial peptide that forms stable channels with well-defined conductance levels. We have used extended molecular dynamics simulations of alamethicin bundles consisting of 4, 5, 6, 7, and 8 helices in a palmitoyl-oleolyl-phosphatidylcholine bilayer to evaluate and analyze channel models and to link the models to the experimentally measured conductance levels. Our results suggest that four helices do not form a stable water-filled channel and might not even form a stable intermediate. The lowest measurable conductance level is likely to correspond to the pentamer. At higher aggregation numbers the bundles become less symmetrical. Water properties inside the different-sized bundles are similar. The hexamer is the most stable model with a stability comparable with simulations based on crystal structures. The simulation was extended from 4 to 20 ns or several times the mean passage time of an ion. Essential dynamics analyses were used to test the hypothesis that correlated motions of the helical bundles account for high-frequency noise observed in open channel measurements. In a 20-ns simulation of a hexameric alamethicin bundle, the main motions are those of individual helices, not of the bundle as a whole. A detailed comparison of simulations using different methods to treat long-range electrostatic interactions (a twin range cutoff, Particle Mesh Ewald, and a twin range cutoff combined with a reaction field correction) shows that water orientation inside the alamethicin channels is sensitive to the algorithms used. In all cases, water ordering due to the protein structure is strong, although the exact profile changes somewhat. Adding an extra 4-nm layer of water only changes the water ordering slightly in the case of particle mesh Ewald, suggesting that periodicity artifacts for this system are not serious.     

Coarse-grained simulations of the salt dependence of the radius of gyration of polyelectrolytes as models for biomolecules in aqueous solution

The salt dependent radius of gyration of a polyelectrolyte in aqueous solution is calculated in an environment where the polyelectrolyte is surrounded by a permeable membrane that exchanges only solvent particles with the bulk. We obtain additionally the scaling exponent of the gyration radius as a function of the polymerization degree, and find that the polyelectrolyte retains a stretched conformation during the condensation and re-expansion process, indicating that these effects are of an electrostatic nature. The solvent quality is also shown to affect the polyelectrolyte conformation, especially for the poor solvent case. These results are obtained using a hybridized Monte Carlo technique with the coarse-grained, dissipative particle dynamics method with fluctuating number of solvent particles. The full range of the electrostatic interactions is included in the simulations, using the Ewald sum method, and the counterions and solvent molecules are included explicitly. In the complex systems mentioned above, the electrostatic interactions and the solvent quality play a key role in understanding phenomena that do not occur in uncharged systems. Our results are compared and validated with the behavior of some biomolecules under similar environments.     

Hydration effects on the electrostatic potential around tuftsin

The electrostatic potential and component dielectric constants from molecular dynamics (MD) trajectories of tuftsin, a tetrapeptide with the amino acid sequence Thr–Lys–Pro–Arg in water and in saline solution are presented. The results obtained from the analysis of the MD trajectories for the total electrostatic potential at points on a grid using the Ewald technique are compared with the solution to the Poisson–Boltzmann (PB) equation. The latter was solved using several sets of dielectric constant parameters. The effects of structural averaging on the PB results were also considered. Solute conformational mobility in simulations gives rise to an electrostatic potential map around the solute dominated by the solute monopole (or lowest order multipole). The detailed spatial variation of the electrostatic potential on the molecular surface brought about by the compounded effects of the distribution of water and ions close to the peptide, solvent mobility, and solute conformational mobility are not qualitatively reproducible from a reparametrization of the input solute and solvent dielectric constants to the PB equation for a single structure or for structurally averaged PB calculations. Nevertheless, by fitting the PB to the MD electrostatic potential surfaces with the dielectric constants as fitting parameters, we found that the values that give the best fit are the values calculated from the MD trajectories. Implications of using such field calculations on the design of tuftsin peptide analogues are discussed. © 1999 John Wiley & Sons, Inc. Biopoly 50: 133–143, 1999     

Effect of periodic box size on aqueous molecular dynamics simulation of a DNA dodecamer with particle-mesh Ewald method

The particle-mesh Ewald (PME) method is considered to be both efficient and accurate for the evaluation of long-range electrostatic interactions in large macromolecular systems being studied by molecular dynamics simulations. This method assumes "infinite" periodic boundary conditions resembling the symmetry of a crystal environment. Can such a "solid-state" method accurately portray a macromolecular solute such as DNA in solution? To address this issue, we have performed three 1500-ps PME molecular dynamics (MD) simulations, each with a different box size, on the d(CGCGA6CG)-(CGT6CGCG) DNA dodecamer. The smallest box had the DNA solvated by a layer of water molecules of at least 5 A along each orthogonal direction. The intermediate size box and the largest box had the DNA solvated by a layer of water molecules of at least 10 A and 15 A, respectively, along each orthogonal direction. The intermediate size box in the present study is similar to the box size currently chosen by most workers in the field. Based on a comparison of RMSDs and curvature for this single DNA dodecamer sequence, the larger two box sizes do not appear to afford any extra benefit over the smallest box. The implications of this finding are briefly discussed.     

Hydration and conformational equilibria of simple hydrophobic and amphiphilic solutes.

We consider whether the continuum model of hydration optimized to reproduce vacuum-to-water transfer free energies simultaneously describes the hydration free energy contributions to conformational equilibria of the same solutes in water. To this end, transfer and conformational free energies of idealized hydrophobic and amphiphilic solutes in water are calculated from explicit water simulations and compared to continuum model predictions. As benchmark hydrophobic solutes, we examine the hydration of linear alkanes from methane through hexane. Amphiphilic solutes were created by adding a charge of +/-1e to a terminal methyl group of butane. We find that phenomenological continuum parameters fit to transfer free energies are significantly different from those fit to conformational free energies of our model solutes. This difference is attributed to continuum model parameters that depend on solute conformation in water, and leads to effective values for the free energy/surface area coefficient and Born radii that best describe conformational equilibrium. In light of these results, we believe that continuum models of hydration optimized to fit transfer free energies do not accurately capture the balance between hydrophobic and electrostatic contributions that determines the solute conformational state in aqueous solution.     

Cutoff size does strongly influence molecular dynamics results on solvated polypeptides.

The behavior of a 17-residue model peptide is analyzed by means of molecular dynamics simulations including explicitly more than a thousand water molecules. On the basis of the charge-group concept, Coulomb interactions are truncated for three values of the cutoff radius: 0.6, 1.0, and 1.4 nm. It is found that the stability of an alpha-helix, which acts as a common starting configuration, is a function of the cutoff size. While the overall stability of the helix is conserved in a simulation using a cutoff of 1.0 nm, it is lost within a very short period of 100 ps when the cutoff is increased to 1.4 nm. This demonstrates that the commonly used cutoff size of 1.0 nm is inappropriate because it does not ensure the convergence of Coulomb interactions. In order to permit an independent judgment, we have performed a 225-ps simulation using the Ewald summation technique, which is more elaborate but circumvents the problem to find an appropriate cutoff value. In contrast to the 1.4-nm cutoff trajectory, the Ewald technique simulation conserves the helical character of the peptide conformation. This demonstrates that even 1.4 nm is too short a cutoff. Due to the fundamental uncertainty introduced by the use of a simple cutoff, this truncation scheme seems questionable for molecular dynamics simulations of solvated biomolecules.     

Reconciling the "old" and "new" views of protein allostery: a molecular simulation study of chemotaxis Y protein (CheY)

A combination of thirty-two 10-ns-scale molecular dynamics simulations were used to explore the coupling between conformational transition and phosphorylation in the bacteria chemotaxis Y protein (CheY), as a simple but representative example of protein allostery. Results from these simulations support an activation mechanism in which the beta4-alpha4 loop, at least partially, gates the isomerization of Tyr106. The roles of phosphorylation and the conserved Thr87 are deemed indirect in that they stabilize the active configuration of the beta4-alpha4 loop. The indirect role of the activation event (phosphorylation) and/or conserved residues in stabilizing, rather than causing, specific conformational transition is likely a feature in many signaling systems. The current analysis of CheY also helps to make clear that neither the "old" (induced fit) nor the "new" (population shift) views for protein allostery are complete, because they emphasize the kinetic (mechanistic) and thermodynamic aspects of allosteric transitions, respectively. In this regard, an issue that warrants further analysis concerns the interplay of concerted collective motion and sequential local structural changes in modulating cooperativity between distant sites in biomolecules.     

On Multiple Time-step Algorithms and the Ewald Sum

Abstract

We show how standard multiple time-step algorithms devised for systems with short-range potentials can be used successfully in simulations of periodic systems with long-range (Coulombic) potentials. Three strategies for incorporating the Ewald sum into a multiple time-step algorithm are considered. These are (i) evaluation of reciprocal space terms every time-step (ii) evaluating reciprocal space terms once every n time-steps and placing these terms in with the slowly varying forces and energies (iii) a modified form of the second strategy in which primary shell (close) electrostatic interactions are evaluated directly and the more distant interactions handled by the Ewald sum (once every n time-steps). Only the first and third approaches give satisfactory thermodynamic results. The third strategy is much more efficient than the first. With the third strategy substantial savings in cpu time are acheived in both the real space and, most importantly, the reciprocal space terms of the Ewald sum. This is achieved without significant loss of accuracy or stability. Overall execution time is decreased by a factor of between 2 and 3.     

Advanced nuclear magnetic resonance lanthanide probe analyses of short-range conformational interrelations controlling ribonucleic acid structures

An advanced method was developed for lanthanide-probe analyses of the conformations of flexible biomolecules such as nucleotides. The new method is to determine structure parameters (such as internal-rotation angles) and population parameters for local conformational equilibria of flexible sites, together with standard deviations of these parameters. As the prominent advantage of this method, the interrelations among local conformations of flexible sites may be quantitatively elucidated from the experimental data of lanthanide-induced shifts and relaxations and vicinal coupling constants. As a structural unit of ribonucleic acids, the molecular conformations and conformational equilibria of uridine 3'-monophosphate in aqueous solution were analyzed. The stable local conformers about the C3'-O3' bond are the G+ (phi' = 281 +/- 11 degrees) and G- (phi' = 211 +/- 8 degrees) forms. The internal rotation about the C3'-O3' bond and the ribose-ring puckering are interrelated; 97 +/- 5% of the C3'-endo ribose ring is associated with the G- form while 70 +/- 22% o the C2'-endo ribose ring is associated with the G+ form. An interdependency also exists between the internal rotation about the C4'-C5' bond and the ribose-ring puckering. These short-range conformational interrelations are probably important in controlling the dynamic aspects of ribonucleic acid structures.     

Molecular dynamics simulations of lipid bilayers: major artifacts due to truncating electrostatic interactions

We study the influence of truncating the electrostatic interactions in a fully hydrated pure dipalmitoylphosphatidylcholine (DPPC) bilayer through 20 ns molecular dynamics simulations. The computations in which the electrostatic interactions were truncated are compared to similar simulations using the particle-mesh Ewald (PME) technique. All examined truncation distances (1.8-2.5 nm) lead to major effects on the bilayer properties, such as enhanced order of acyl chains together with decreased areas per lipid. The results obtained using PME, on the other hand, are consistent with experiments. These artifacts are interpreted in terms of radial distribution functions g(r) of molecules and molecular groups in the bilayer plane. Pronounced maxima or minima in g(r) appear exactly at the cutoff distance indicating that the truncation gives rise to artificial ordering between the polar phosphatidyl and choline groups of the DPPC molecules. In systems described using PME, such artificial ordering is not present.     

Nucleotide-dependent conformational changes in HisP: molecular dynamics simulations of an ABC transporter nucleotide-binding domain

ATP-binding cassette (ABC) transporters mediate the movement of molecules across cell membranes in both prokaryotes and eukaryotes. In ABC transporters, solute translocation occurs after ATP is either bound or hydrolyzed at the intracellular nucleotide-binding domains (NBDs). Molecular dynamics (MD) simulations have been employed to study the interactions of nucleotide with NBD. The results of extended (approximately 20 ns) MD simulations of HisP (total simulation time approximately 80 ns), the NBD of the histidine transporter HisQMP2J from Salmonella typhimurium, are presented. Analysis of the MD trajectories reveals conformational changes within HisP that are dependent on the presence of ATP in the binding pocket of the protein, and are sensitive to the presence/absence of Mg ions bound to the ATP. These changes are predominantly confined to the alpha-helical subdomain of HisP. Specifically there is a rotation of three alpha-helices within the subdomain, and a movement of the signature sequence toward the bound nucleotide. In addition, there is considerable conformational flexibility in a conserved glutamine-containing loop, which is situated at the interface between the alpha-helical subdomain and the F1-like subdomain. These results support the mechanism for ATP-induced conformational transitions derived from the crystal structures of other NBDs.     

Monte-Carlo simulations of strongly interacting dipolar fluids between two conducting walls

We report Monte-Carlo simulation results for a strongly interacting dipolar soft sphere (DSS) fluid confined between two conducting, planar walls. The long-range dipolar interactions, including contributions from the “image dipoles” in the metal, are handled by mapping onto a problem with three-dimensional (3d) periodicity which can be treated by conventional Ewald summation methods. Considering two different wall separations our results indicate the occurence of wall-induced local and long-range ordering very similar to the related case of insulating walls. To understand this behavior we present ground-state (lattice) calculations on the basis of the Ewald sums, as well as simple macroscopic arguments appropriate for dipolar systems between conducting walls.     

Continuum solvent studies of the stability of RNA hairpin loops and helices

We apply continuum solvent models to investigate the relative stability of various conformational forms for two RNA sequences, GGAC(UUCG)GUCC and GGUG(UGAA)CACC. In the first part, we compare alternate hairpin conformations to explore the reliability of these models to discriminate between different local conformations. A second part looks at the hairpin-duplex conversion for the UUCG sequence, identifying major contributors to the thermodynamics of a much large scale transition. Structures were taken as snapshots from multi-nanosecond molecular dynamics simulations computed in a consistent fashion using explicit solvent and with long-range electrostatics accounted for using the Particle-Mesh Ewald procedure. The electrostatic contribution to solvation energies were computed using both a finite-difference Poisson-Boltzmann (PB) model and a pairwise Generalized Born model; non-electrostatic contributions were estimated with a surface-area dependent term. To these solvation free energies were added the mean solute internal energies (determined from a molecular mechanics potential) and estimates of the solute entropy (from a harmonic analysis). Consistent with experiment and with earlier solvated molecular dynamics simulations, the UUCG hairpin was found to prefer conformers close to a recent NMR structure determination in preference to those from an earlier NMR study. Similarly, results for the UGAA hairpin favored an NMR-derived structure over that to be expected for a generic GNRA hairpin loop. Experimental free energies are not known for the hairpin/duplex conversion, but must be close to zero since hairpins are seen in solution and duplexes in crystals; out calculations find a value near zero and illustrate the expected interplay of solvation, salt effects and entropy in affecting this equilibrium.     

An introduction to biomolecular simulations and docking

The biomolecules in and around a living cell – proteins, nucleic acids, lipids and carbohydrates – continuously sample myriad conformational states that are thermally accessible at physiological temperatures. Simultaneously, a given biomolecule also samples (and is sampled by) a rapidly fluctuating local environment comprising other biopolymers, small molecules, water, ions, etc. that diffuse to within a few nanometres, leading to inter-molecular contacts that stitch together large supramolecular assemblies. Indeed, all biological systems can be viewed as dynamic networks of molecular interactions. As a complement to experimentation, molecular simulation offers a uniquely powerful approach to analyse biomolecular structure, mechanism and dynamics; this is possible because the molecular contacts that define a complicated biomolecular system are governed by the same physical principles (forces and energetics) that characterise individual small molecules, and these simpler systems are relatively well-understood. With modern algorithms and computing capabilities, simulations are now an indispensable tool for examining biomolecular assemblies in atomic detail, from the conformational motion in an individual protein to the diffusional dynamics and inter-molecular collisions in the early stages of formation of cellular-scale assemblies such as the ribosome. This text introduces the physicochemical foundations of molecular simulations and docking, largely from the perspective of biomolecular interactions.