Collective motion artifacts arising in long-duration molecular dynamics simulations

Molecular Modeling of Water-in-Salt Electrolytes: A Comprehensive Analysis of Polarization Effects and Force Field Parameters in Molecular Dynamics Simulations

Accurate modeling of highly concentrated aqueous solutions, such as water-in-salt (WiS) electrolytes in battery applications, requires proper consideration of polarization contributions to atomic interactions. Within the force field molecular dynamics (MD) simulations, the atomic polarization can be accounted for at various levels. Nonpolarizable force fields implicitly account for polarization effects by incorporating them into their van der Waals interaction parameters. They can additionally mimic electron polarization within a mean-field approximation through ionic charge scaling. Alternatively, explicit polarization description methods, such as the Drude oscillator model, can be selectively applied to either a subset of polarizable atoms or all polarizable atoms to enhance simulation accuracy. The trade-off between simulation accuracy and computational efficiency highlights the importance of determining an optimal level of accounting for atomic polarization. In this study, we analyze different approaches to include polarization effects in MD simulations of WiS electrolytes, with an example of a Na-OTF solution. These approaches range from a nonpolarizable to a fully polarizable force field. After careful examination of computational costs, simulation stability, and feasibility of controlling the electrolyte properties, we identify an efficient combination of force fields: the Drude polarizable force field for salt ions and non-polarizable models for water. This cost-effective combination is sufficiently flexible to reproduce a broad range of electrolyte properties, while ensuring simulation stability over a relatively wide range of force field parameters. Furthermore, we conduct a thorough evaluation of the influence of various force field parameters on both the simulation results and technical requirements, with the aim of establishing a general framework for force field optimization and facilitating parametrization of similar systems.

Anomalous Effects of Velocity Rescaling Algorithms: The Flying Ice Cube Effect Revisited

The flying ice cube effect is a molecular dynamics simulation artifact in which the use of velocity rescaling thermostats sometimes causes violation of the equipartition theorem, affecting both structural and dynamic properties. The reason for this artifact and the conditions under which it occurs have not been fully understood. Since the flying ice cube effect was first demonstrated, a new velocity rescaling algorithm (the CSVR thermostat) has been developed and become popular without its effects on the equipartition theorem being truly known. Meanwhile, the use of simple velocity rescaling and Berendsen (weak coupling) thermostat algorithms has not abated but has actually continued to grow. Here, we have calculated the partitioning of the kinetic energy between translational, rotational, and vibrational modes in simulations of diatomic molecules to explicitly determine whether the equipartition theorem is violated under different thermostats and while rescaling velocities to different kinetic energy distributions. We have found that the underlying cause of the flying ice cube effect is a violation of balance leading to systematic redistributions of kinetic energy under simple velocity rescaling and the Berendsen thermostat. When velocities are instead rescaled to the canonical ensemble's kinetic energy distribution, as is done with the CSVR thermostat, the equipartition theorem is not violated, and we show that the CSVR thermostat satisfies detailed balance. The critical necessity for molecular dynamics practitioners to abandon the use of popular yet incorrect velocity rescaling algorithms is underscored with an example demonstrating that the main result of a highly cited study is entirely due to artifacts resulting from the study's use of the Berendsen thermostat.

The good, the bad and the user in soft matter simulations

Molecular dynamics (MD) simulations have become popular in materials science, biochemistry, biophysics and several other fields. Improvements in computational resources, in quality of force field parameters and algorithms have yielded significant improvements in performance and reliability. On the other hand, no method of research is error free. In this review, we discuss a few examples of errors and artifacts due to various sources and discuss how to avoid them. Besides bringing attention to artifacts and proper practices in simulations, we also aim to provide the reader with a starting point to explore these issues further. In particular, we hope that the discussion encourages researchers to check software, parameters, protocols and, most importantly, their own practices in order to minimize the possibility of errors. The focus here is on practical issues. This article is part of a Special Issue entitled: Biosimulations edited by Ilpo Vattulainen and Tomasz Róg.

The "Hot-Solvent/Cold-Solute" Problem Revisited

The temperature steers the equilibrium and nonequilibrium conformational dynamics of macromolecules in solution. Therefore, corresponding molecular dynamics simulations require a strategy for temperature control which should guarantee that the experimental statistical ensemble is also sampled in silico. Several algorithms for temperature control have been proposed. All these thermostats interfere with the macromolecule's "natural" dynamics as given by the Newtonian mechanics. Furthermore, using a single thermostat for an inhomogeneous solute-solvent system can lead to stationary temperature gradients. To avoid this "hot solvent/cold solute" problem, two separate thermostats are frequently applied, one to the solute and one to the solvent. However, such a separate temperature control will perturb the dynamics of the macromolecule much more strongly than a global one and, therefore, can introduce large artifacts into its conformational dynamics. Based on the concept that an explicit solvent environment represents an ideal thermostat concerning the magnitude and time correlation of temperature fluctuations of the solute, we propose a temperature control strategy that, on the one hand, provides a homogeneous temperature distribution throughout the system together with the correct statistical ensemble for the solute molecule while, on the other hand, minimally perturbing its dynamics.

Mechanism of C-Terminal Fragments of Amyloid β-Protein as Aβ Inhibitors: Do C-Terminal Interactions Play a Key Role in Their Inhibitory Activity?

Targeting the early oligomerization of amyloid β protein (Aβ) is a promising therapeutic strategy for Alzheimer's disease (AD). Recently, certain C-terminal fragments (CTFs) derived from Aβ42 were shown to be potent inhibitors of Aβ-induced toxicity. The shortest peptide studied, Aβ(39-42), has been shown to modulate Aβ oligomerization and inhibit Aβ toxicity. Understanding the mechanism of these CTFs, especially Aβ(39-42), is of significance for future therapeutic development of AD and peptidomimetic-based drug development. Here we used ion mobility spectrometry-mass spectrometry to investigate the interactions between two modified Aβ(39-42) derivatives, VVIA-NH2 and Ac-VVIA, and full-length Aβ42. VVIA-NH2 was previously shown to inhibit Aβ toxicity, whereas Ac-VVIA did not. Our mass spectrometry analysis revealed that VVIA-NH2 binds directly to Aβ42 monomer and small oligomers while Ac-VVIA binds only to Aβ42 monomer. Ion mobility studies showed that VVIA-NH2 modulates Aβ42 oligomerization by not only inhibiting the dodecamer formation but also disaggregating preformed Aβ42 dodecamer. Ac-VVIA also inhibits and removes preformed Aβ42 dodecamer. However, the Aβ42 sample with the addition of Ac-VVIA clogged the nanospray tip easily, indicating that larger aggregates are formed in the solution in the presence of Ac-VVIA. Molecular dynamics simulations suggested that VVIA-NH2 binds specifically to the C-terminal region of Aβ42 while Ac-VVIA binds dispersedly to multiple regions of Aβ42. This work implies that C-terminal interactions and binding to Aβ oligomers are important for C-terminal fragment inhibitors.

An adaptive finite-element method for large-scale ab initio molecular dynamics simulations

A snapshot of ab initio molecular dynamics simulations for a polymer electrolyte membrane at low hydration.

Dynamics of lipids, cholesterol, and transmembrane α-helices from microsecond molecular dynamics simulations

Extensive all-atom molecular dynamics (∼24 μs total) allowed exploration of configurational space and calculation of lateral diffusion coefficients of the components of a protein-embedded, cholesterol-containing model bilayer. The three model membranes are composed of an ∼50/50 (by mole) dipalmitoylphosphatidylcholine (DPPC)/cholesterol bilayer and contained an α-helical transmembrane protein (HIV-1 gp41 TM). Despite the high concentration of cholesterol, normal Brownian motion was observed and the calculated diffusion coefficients (on the order of 10(-9) cm(2)/s) are consistent with experiments. Diffusion is sensitive to a variety of parameters, and a temperature difference of ∼4 K from thermostat artifacts resulted in 2-10-fold differences in diffusion coefficients and significant differences in lipid order, membrane thickness, and unit cell area. Also, the specific peptide sequence likely underlies the consistently observed faster diffusion in one leaflet. Although the simulations here present molecular dynamics (MD) an order of magnitude longer than those from previous studies, the three systems did not approach ergodicity. The distributions of cholesterol and DPPC around the peptides changed on the microsecond time scale, but not significantly enough to thoroughly explore configurational space. These simulations support conclusions of other recent microsecond MD in that even longer time scales are needed for equilibration of model membranes and simulations of more realistic cellular or viral bilayers.

Molecular modeling of nucleic acid structure: energy and sampling

An overview of computer simulation techniques as applied to nucleic acid systems is presented. This unit discusses methods used to treat energy and to sample representative configurations. Emphasis is placed on molecular mechanics and empirical force fields.

Elucidating the evolutionary conserved DNA-binding specificities of WRKY transcription factors by molecular dynamics and in vitro binding assays

WRKY transcription factors constitute a large protein family in plants that is involved in the regulation of developmental processes and responses to biotic or abiotic stimuli. The question arises how stimulus-specific responses are mediated given that the highly conserved WRKY DNA-binding domain (DBD) exclusively recognizes the 'TTGACY' W-box consensus. We speculated that the W-box consensus might be more degenerate and yet undetected differences in the W-box consensus of WRKYs of different evolutionary descent exist. The phylogenetic analysis of WRKY DBDs suggests that they evolved from an ancestral group IIc-like WRKY early in the eukaryote lineage. A direct descent of group IIc WRKYs supports a monophyletic origin of all other group II and III WRKYs from group I by loss of an N-terminal DBD. Group I WRKYs are of paraphyletic descent and evolved multiple times independently. By homology modeling, molecular dynamics simulations and in vitro DNA-protein interaction-enzyme-linked immunosorbent assay with AtWRKY50 (IIc), AtWRKY33 (I) and AtWRKY11 (IId) DBDs, we revealed differences in DNA-binding specificities. Our data imply that other components are essentially required besides the W-box-specific binding to DNA to facilitate a stimulus-specific WRKY function.

Binding of Congo red to amyloid protofibrils of the Alzheimer Aβ(9-40) peptide probed by molecular dynamics simulations

Congo red (CR) is a commonly used histological amyloid dye and a weak amyloid inhibitor. There is currently no experimentally available structure of CR bound to an amyloid fibril and the binding modes, and the mechanisms governing its inhibitory and optical properties are poorly understood. In this work, we present the first, to our knowledge, atomistically detailed picture of CR binding to protofibrils of the Alzheimer Aβ(9-40) peptide. We identify three major binding modes, with the primary mode residing in the grooves formed by the β-sheets, and observe a restriction of the torsional rotation of the CR molecule upon binding. Our simulations reveal a novel, to our knowledge, electrostatic steering mechanism that plays an important role in the initial recognition and binding of CR to the positively charged surface residues of the fibril. Our simulations provide new, to our knowledge, insights into the striking spectrophotometric and inhibitory properties of CR. In particular, we show that birefringence upon CR binding is due to the anisotropic orientation of the CR dipoles resulting from the spatial ordering of these molecules in the grooves along the fibril axis. The fluorescent enhancement of the bound CR, in turn, is associated with the torsional restriction of this molecule upon binding.

Early stage intercalation of doxorubicin to DNA fragments observed in molecular dynamics binding simulations

The intercalation mode between doxorubicin (an anticancer drug) and two 6-base-pair DNA model fragments (d(CGATCG)₂ and d(CGTACG)₂) has been well studied by X-ray crystallography and NMR experimental methods. Yet, the detailed intercalation pathway at molecular level remains elusive. In this study, we conducted molecular dynamics binding simulations of these two systems using AMBER DNA (parmbsc0) and drug (GAFF) force fields starting from the unbound state. We observed outside binding (minor groove binding or end-binding) in all six independent binding simulations (three for each DNA fragment), followed by the complete intercalation of a drug molecule in two simulations (one for each DNA fragment). First, our data directly supported that the minor groove binding is the dominant pre-intercalation step. Second, we observed that the opening and flipping of a local base pair (A3-T10 for d(CGATCG)₂ and C1-G12 for d(CGTACG)₂) in the two intercalation trajectories. This locally cooperative flipping-intercalation mechanism was different from the previously proposed rise-insertion mechanism by which the distance between two neighboring intact base pairs increases to create a space for the drug insertion. Third, our simulations provided the first set of data to support the applicability of the AMBER DNA and drug force fields in drug-DNA atomistic binding simulations. Implications on the kinetics pathway and drug action are also discussed.

Aβ(39-42) modulates Aβ oligomerization but not fibril formation

Recently, certain C-terminal fragments (CTFs) of Aβ42 have been shown to be effective inhibitors of Aβ42 toxicity. Here, we examine the interactions between the shortest CTF in the original series, Aβ(39-42), and full-length Aβ. Mass spectrometry results indicate that Aβ(39-42) binds directly to Aβ monomers and to the n = 2, 4, and 6 oligomers. The Aβ42:Aβ(39-42) complex is further probed using molecular dynamics simulations. Although the CTF was expected to bind to the hydrophobic C-terminus of Aβ42, the simulations show that Aβ(39-42) binds at several locations on Aβ42, including the C-terminus, other hydrophobic regions, and preferentially in the N-terminus. Ion mobility-mass spectrometry (IM-MS) and electron microscopy experiments indicate that Aβ(39-42) disrupts the early assembly of full-length Aβ. Specifically, the ion-mobility results show that Aβ(39-42) prevents the formation of large decamer/dodecamer Aβ42 species and, moreover, can remove these structures from solution. At the same time, thioflavin T fluorescence and electron microscopy results show that the CTF does not inhibit fibril formation, lending strong support to the hypothesis that oligomers and not amyloid fibrils are the Aβ form responsible for toxicity. The results emphasize the role of small, soluble assemblies in Aβ-induced toxicity and suggest that Aβ(39-42) inhibits Aβ-induced toxicity by a unique mechanism, modulating early assembly into nontoxic hetero-oligomers, without preventing fibril formation.

Alzheimer's disease drug candidates stabilize A-β protein native structure by interacting with the hydrophobic core

Deposition of amyloid fibrils, consisting primarily of Aβ(40) and Aβ(42) peptides, in the extracellular space in the brain is a major characteristic of Alzheimer's disease (AD). We recently developed new (to our knowledge) drug candidates for AD that inhibit the fibril formation of Aβ peptides and eliminate their neurotoxicity. We performed all-atom molecular-dynamics simulations on the Aβ(42) monomer at its α-helical conformation and a pentamer fibril fragment of Aβ(42) peptide with or without LRL and fluorene series compounds to investigate the mechanism of inhibition. The results show that the active drug candidates, LRL22 (EC(50) = 0.734 μM) and K162 (EC(50) = 0.080 μM), stabilize hydrophobic core I of Aβ(42) peptide (residues 17-21) to its α-helical conformation by interacting specifically in this region. The nonactive drug candidates, LRL27 (EC(50) > 10 μM) and K182 (EC(50) > 5 μM), have little to no similar effect. This explains the different behavior of the drug candidates in experiments. Of more importance, this phenomenon indicates that hydrophobic core I of the Aβ(42) peptide plays a major mechanistic role in the formation of amyloid fibrils, and paves the way for the development of new drugs against AD.

Molecular structures of quiescently grown and brain-derived polymorphic fibrils of the Alzheimer amyloid abeta9-40 peptide: a comparison to agitated fibrils

The presence of amyloid deposits consisting primarily of Amyloid-β (Aβ) fibril in the brain is a hallmark of Alzheimer's disease (AD). The morphologies of these fibrils are exquisitely sensitive to environmental conditions. Using molecular dynamics simulations combined with data from previously published solid-state NMR experiments, we propose the first atomically detailed structures of two asymmetric polymorphs of the Aβ_9-40 peptide fibril. The first corresponds to synthetic fibrils grown under quiescent conditions and the second to fibrils derived from AD patients' brain-extracts. Our core structure in both fibril structures consists of a layered structure in which three cross-β subunits are arranged in six tightly stacked β-sheet layers with an antiparallel hydrophobic-hydrophobic and an antiparallel polar-polar interface. The synthetic and brain-derived structures differ primarily in the side-chain orientation of one β-strand. The presence of a large and continually exposed hydrophobic surface (buried in the symmetric agitated Aβ fibrils) may account for the higher toxicity of the asymmetric fibrils. Our model explains the effects of external perturbations on the fibril lateral architecture as well as the fibrillogenesis inhibiting action of amphiphilic molecules.

Amyloid diseases are characterized by the presence of amyloid fibrils on organs and tissue in the body. Alzheimer's disease, Parkinson's diseases and Type II Diabetes are all examples of amyloid diseases. Determining the structure of amyloid fibrils is critical for understanding the mechanism of fibril formation as well as for the design of inhibitor molecules that can prevent aggregation. In the case of the Alzheimer Amyloid-β (Aβ) peptide, the structure of fibrils grown under conditions of mechanical agitation has been elucidated from a combination of simulation and experiments. However, the structures of the asymmetric quiescent Aβ fibrils (grown under conditions akin to physiological conditions) and of Alzheimer's brain–derived fibrils are not known. In this paper, we propose the first atomically detailed structures of these two fibrils, using molecular dynamics simulations combined with data from previously published experiments. In additions, we suggest a unifying lateral growth mechanism that explains the increased toxicity of quiescent Aβ fibrils, the effects of external perturbations on fibril lateral architecture and the inhibition mechanism of the small molecule inhibitors on fibril formation.

Novel procedure for thermal equilibration in molecular dynamics simulation

We describe a simple novel procedure for achieving thermal equilibration between a protein and a surrounding solvent during molecular dynamics (MD) simulation. The method uniquely defines the length of simulation time required to achieve thermal equilibrium over a broad range of parameters, thus removing ambiguities associated with the traditional heuristic approaches. The proposed protocol saves simulation time and avoids bias introduced by the inclusion of non-equilibrium events. The key element of the procedure involves coupling only the solvent atoms to a standard heat bath. Measuring progress towards thermal equilibration involves simply monitoring the difference in temperature between the solvent and the protein. Here, we report that the results of MD simulations using the above procedure are measurably improved relative to the traditional approaches in terms of root-mean-square deviations and principal components analysis both indicating significantly less undesirable divergence.

WHAT MAKES MOLECULAR DYNAMICS WORK?

The equations of motion for deterministic molecular dynamics (MD) are chaotic, creating problems for their numerical treatment due to the exponential growth of error with time. Indeed, modeling and computational errors overwhelm numerical trajectories in typical simulations. Consequently, accuracy is expected only in a statistical sense, based on random initial conditions. Of great interest then is the relationship between errors in the dynamics and their effects on the accuracy of statistical quantities, specifically, expectations. This article provides a formula for the effect of a perturbation on an ensemble average, which explains the accuracy of such calculations. It also provides a formula for the effect of a perturbation on a time correlation function, which, however, fails to explain accuracy for these calculations. Additionally, this article clarifies the relationships among various dynamical properties of MD and provides an extension to a theory of non-Hamiltonian MD.

The binding of thioflavin T and its neutral analog BTA-1 to protofibrils of the Alzheimer's disease Abeta(16-22) peptide probed by molecular dynamics simulations

Thioflavin T (ThT) is a fluorescent dye commonly used to stain amyloid plaques, but the binding sites of this dye onto fibrils are poorly characterized. We present molecular dynamics simulations of the binding of ThT and its neutral analog BTA-1 [2-(4'-methylaminophenyl)benzothiazole] to model protofibrils of the Alzheimer's disease Abeta(16-22) (amyloid beta) peptide. Our simulations reveal two binding modes located at the grooves of the beta-sheet surfaces and at the ends of the beta-sheet. These simulations provide new insight into recent experimental work and allow us to characterize the high-capacity, micromolar-affinity site seen in experiment as binding to the beta-sheet surface grooves and the low-capacity, nanomolar-affinity site seen as binding to the beta-sheet extremities of the fibril. The structure-activity relationship upon mutating charged ThT to neutral BTA-1 in terms of increased lipophilicity and binding affinity was studied, with calculated solvation free energies and binding energies found to be in qualitative agreement with the experimental measurements.

Phenol red interacts with the protofibril-like oligomers of an amyloidogenic hexapeptide NFGAIL through both hydrophobic and aromatic contacts

Amyloid-associated diseases affect millions of people worldwide. Phenol red exhibits modest inhibition toward fibril formation of human Islet amyloid polypeptide (hIAPP) and its toxicity, which is associated with type II diabetes mellitus. However, the molecular level mechanisms of interactions remain elusive. The binding of phenol red molecules to the protofibrils of an amyloidogenic fragment (NFGAIL) of hIAPP has been investigated by molecular dynamics simulations with explicit solvent. The phenol red molecules were observed to bind primarily along either beta-sheet stacking or beta-strand directions. Through its three aromatic rings, the phenol red molecule preferentially interacted with the hydrophobic side chains of Phe, Leu, and Ile; and the polar sulfone and hydroxyl groups were mainly exposed in solvent. Thus, phenol red improves the solubility of the early protofibrils and represses further growth. Interestingly, there was no obvious preference toward the aromatic Phe residue in comparison to the hydrophobic Leu or Ile residues. The lack of binding along the hydrogen bond direction indicates that phenol red does not directly block the beta-sheet extension. Further free energy analysis suggested that a phenol red analog may potentially improve the binding affinity.

Large-amplitude fast motions in double-stranded DNA driven by solvent thermal fluctuations

The nature of the internal dynamics of double-stranded DNA in aqueous environment remains to be established. We consider the motions to stem from thermal fluctuations/dissipations of the harmonic modes of beads (bases and sugars) in a cylindrical geometry that are tracked through the stochastic Langevin trajectories; these are characterized by parameters obtained from published data. The present approach has allowed a comparative study of the dynamics for DNA lengths in the range of 20-600 base pairs. For this range, we find that rotational motions about directions parallel to the helix axis (opening, twist) and perpendicular to it (propeller-twist, roll) contribute significantly to the dynamics. For a 20-mer at a solvent viscosity of 1 cP, the calculated fluorescence anisotropy profile exhibits a fast decay in the subnanosecond range due to large-amplitude fluctuations at the mesoscopic level. This feature reproduces the experimental behavior well, and suggests a possible way for the initiation of biological processes: they may be suddenly triggered on this scale through the occurrence of favorable thermal fluctuations. This analysis also reveals that, as is the case for a 20-mer, the dynamics of longer N-mers are dominated by internal motions, and are modulated by the viscosity of the solvent, in agreement with our previous experimental observations. Moreover, the model indicates that occurrence of partially concerted rotations of the bases due to thermal fluctuations can possibly be sustained over a DNA length of the order of 100 A at 1 ns, suggesting a possible mechanism for action-at-a-distance in transcription.

An improved nucleic acid parameter set for the GROMOS force field

Over the past decades, the GROMOS force field for biomolecular simulation has primarily been developed for performing molecular dynamics (MD) simulations of polypeptides and, to a lesser extent, sugars. When applied to DNA, the 43A1 and 45A3 parameter sets of the years 1996 and 2001 produced rather flexible double-helical structures, in which the Watson-Crick hydrogen-bonding content was more limited than expected. To improve on the currently available parameter sets, the nucleotide backbone torsional-angle parameters and the charge distribution of the nucleotide bases are reconsidered based on quantum-chemical data. The new 45A4 parameter set resulting from this refinement appears to perform well in terms of reproducing solution NMR data and canonical hydrogen bonding. The deviation between simulated and experimental observables is now of the same order of magnitude as the uncertainty in the experimental values themselves.

DNA basepair step deformability inferred from molecular dynamics simulations

The sequence-dependent DNA deformability at the basepair step level was investigated using large-scale atomic resolution molecular dynamics simulation of two 18-bp DNA oligomers: d(GCCTATAAACGCCTATAA) and d(CTAGGTGGATGACTCATT). From an analysis of the structural fluctuations, the harmonic potential energy functions for all 10 unique steps with respect to the six step parameters have been evaluated. In the case of roll, three distinct groups of steps have been identified: the flexible pyrimidine-purine (YR) steps, intermediate purine-purine (RR), and stiff purine-pyrimidine (RY). The YR steps appear to be the most flexible in tilt and partially in twist. Increasing stiffness from YR through RR to RY was observed for rise, whereas shift and slide lack simple trends. A proposed measure of the relative importance of couplings identifies the slide-rise, twist-roll, and twist-slide couplings to play a major role. The force constants obtained are of similar magnitudes to those based on a crystallographic ensemble. However, the current data have a less complicated and less pronounced sequence dependence. A correlation analysis reveals concerted motions of neighboring steps and thus exposes limitations in the dinucleotide model. The comparison of DNA deformability from this and other studies with recent quantum-chemical stacking energy calculations suggests poor correlation between the stacking and flexibility.

Formation pathways of a guanine-quadruplex DNA revealed by molecular dynamics and thermodynamic analysis of the substates

The formation of a cation-stabilized guanine quadruplex (G-DNA) stem is an exceptionally slow process involving complex kinetics that has not yet been characterized at atomic resolution. Here, we investigate the formation of a parallel stranded G-DNA stem consisting of four strands of d(GGGG) using molecular dynamics simulations with explicit inclusion of counterions and solvent. Due to the limitations imposed by the nanosecond timescale of the simulations, rather than watching for the spontaneous formation of G-DNA, our approach probes the stability of possible supramolecular intermediates (including two-, three-, and four-stranded assemblies with out-of-register base pairing between guanines) on the formation pathway. The simulations suggest that "cross-like" two-stranded assemblies may serve as nucleation centers in the initial formation of parallel stranded G-DNA quadruplexes, proceeding through a series of rearrangements involving trapping of cations, association of additional strands, and progressive slippage of strands toward the full stem. To supplement the analysis, approximate free energies of the models are obtained with explicit consideration of the integral cations. The approach applied here serves as a prototype for qualitatively investigating other G-DNA molecules using molecular dynamics simulation and free-energy analysis.

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.

Structure and dynamics of liquid water with different long-range interaction truncation and temperature control methods in molecular dynamics simulations

We have used molecular dynamics simulations to study the physical properties of modified TIP3P water model included in the CHARMM program, using four different methods-the Ewald summation technique, and three different spherical truncation methods-for the treatment of the long-range interactions. Both the structure and dynamics of the liquid water model were affected by the methods used to truncate the long-range interactions. For some of the methods artificial structuring of the model liquid was observed around the cutoff radius. The model liquid properties were also affected by the commonly applied temperature control methods. Four different methods for controlling the temperature of the system were studied, and the effects of these methods on the bulk properties for liquid water were analyzed. The system size was also found to change the dynamics of the model liquid water. Two control simulations with the SPC/E water model were carried out. The self-diffusion coefficient (D), the radial distribution function (g(OO)), the distance dependent Kirkwood G-factor [G(k)(r)] and the intermolecular potential energy (E(pot)) were determined from the different trajectories and compared with the experimental data.

Critical effect of the N2 amino group on structure, dynamics, and elasticity of DNA polypurine tracts

Unrestrained 5-20-ns explicit-solvent molecular dynamics simulations using the Cornell et al. force field have been carried out for d[GCG(N)11GCG]2 (N, purine base) considering guanine*cytosine (G*C), adenine*thymine (A*T), inosine*5-methyl-cytosine (I*mC), and 2-amino-adenine*thymine (D*T) basepairs. The simulations unambiguously show that the structure and elasticity of N-tracts is primarily determined by the presence of the amino group in the minor groove. Simulated A-, I-, and AI-tracts show almost identical structures, with high propeller twist and minor groove narrowing. G- and D-tracts have small propeller twisting and are partly shifted toward the A-form. The elastic properties also differ between the two groups. The sequence-dependent electrostatic component of base stacking seems to play a minor role. Our conclusions are entirely consistent with available experimental data. Nevertheless, the propeller twist and helical twist in the simulated A-tract appear to be underestimated compared to crystallographic studies. To obtain further insight into the possible force field deficiencies, additional multiple simulations have been made for d(A)10, systematically comparing four major force fields currently used in DNA simulations and utilizing B and A-DNA forms as the starting structure. This comparison shows that the conclusions of the present work are not influenced by the force field choice.

Molecular dynamics simulation of nucleic acids: successes, limitations, and promise

In the last five years we have witnessed a significant increase in the number publications describing accurate and reliable all-atom molecular dynamics simulations of nucleic acids. This increase has been facilitated by the development of fast and efficient methods for treating the long-range electrostatic interactions, the availability of faster parallel computers, and the development of well-validated empirical molecular mechanical force fields. With these technologies, it has been demonstrated that simulation is not only capable of consistently reproducing experimental observations of sequence specific fine structure of DNA, but also can give detailed insight into prevalent problems in nucleic acid structure, ion association and specific hydration of nucleic acids, polyadenine tract bending, and the subtle environmental dependence of the A-DNA-B-DNA duplex equilibrium. Despite the advances, there are still issues with the methods that need to be resolved through rigorous controlled testing. In general, these relate to deficiencies of the underlying molecular mechanical potentials or applied methods (such as the imposition of true periodicity in Ewald simulations and the need for energy conservation), and significant limits in effective conformational sampling. In this perspective, we provide an overview of our experiences, provide some cautionary notes, and provide recommendations for further study in molecular dynamics simulation of nucleic acids.

Molecular dynamics simulation of the human U2B" protein complex with U2 snRNA hairpin IV in aqueous solution

A 2200-ps molecular dynamics (MD) simulation of the U2 snRNA hairpin IV/U2B" complex was performed in aqueous solution using the particle mesh Ewald method to consider long-range electrostatic interactions. To investigate the interaction and recognition process between the RNA and protein, the free energy contributions resulting from individual amino acids of the protein component of the RNA/protein complex were calculated using the recently developed glycine-scanning method. The results revealed that the loop region of the U2 snRNA hairpin IV interacted mainly with three regions of the U2B" protein: 1) beta 1-helix A, 2) beta 2-beta 3, and 3) beta 4-helix C. U2 snRNA hairpin IV bound U2B" in a similar orientation as that previously described for U1 snRNA with the U1A' protein; however, the details of the interaction differed in several aspects. In particular, beta 1-helix A and beta 4-helix C in U2B" were not observed to interact with RNA in the U1A' protein complex. Most of the polar and charged residues in the interacting regions had larger mutant free energies than the nonpolar residues, indicating that electrostatic interactions were important for stabilizing the RNA/protein complex. The interaction was further stabilized by a network of hydrogen bonds and salt bridges formed between RNA and protein that was maintained throughout the MD trajectory. In addition to the direct interactions between RNA and the protein, solvent-mediated interactions also contributed significantly to complex stability. A detailed analysis of the ordered water molecules in the hydration of the RNA/protein complex revealed that bridged water molecules reside at the interface of RNA and protein as long as 2100 ps in the 2200-ps trajectory. At least 20 bridged water molecules, on average, contributed to the instantaneous stability of the RNA/protein complex. The stabilizing interaction energy due to bridging water molecules was obtained from ab initio Hartree-Fock and density functional theory calculations.

Molecular modeling of nucleic acid structure: energy and sampling

An overview of computer simulation techniques as applied to nucleic acid systems is presented. This unit expands an accompanying overview unit (UNIT 7.5) by discussing methods used to treat the energy and sample representative configurations. Emphasis is placed on molecular mechanics and empirical force fields.