Effect of Warmup Protocol and Sampling Time on Convergence of Molecular Dynamics Simulations of a DNA Dodecamer Using AMBER 4.1 and Particle-Mesh Ewald Method

Abstract

This report describes one 3000 ps and two 1500 ps molecular dynamic simulations on a TATA box containing dodecamer DNA duplex in a periodic box of TIP3P water molecules, using the AMBER 4.1 implementation of the particle-mesh Ewald method. We compare the effect of warmup protocol and simulation time length on the root-mean square deviation (RMSD) parameter. For the longer simulation, the RMSD computed for the 500–1000 ps time interval is representative of longer time intervals, including 500–3000 ps. The various warmup protocols do not appear to have a significant effect on the simulation results. Based on the present results, DNA sequence-dependent differences in RMSD, or related properties, should exceed two standard deviations before being attributed to non-simulation factors, such as warmup protocol and sampling time effects; we recommend a minimum criterion of at least a three standard deviation difference with a sampling period of at least 500–1000 ps. In addition, while end effects appear negligible there is a consistent dependence of RMSD on DNA helix length.     

Lipoplex thermodynamics: determination of DNA-cationic lipoid interaction energies

An experimental study of the cationic lipid-DNA binding affinity is presented. The binding free energy was determined by monitoring lipoplex dissociation under conditions of increasing salt concentration. The primary procedure was based on the extent of quenching by energy transfer of fluorophores on DNA molecules by fluorophore on a lipid as these molecules came into close association in the lipoplex. Titration calorimetry on the Dickerson dodecamer was also done, with results that were in agreement with the fluorescence data. Measurements on short oligonucleotides allowed estimation of the binding energy per nucleotide. The binding free energy is approximately 0.6 kcal/mole nucleotide for the Dickerson dodecamer and declines for longer oligonucleotides. The entropy gained upon complex formation is approximately 1 entropy unit per released counterion. The method was applied to long DNA molecules (herring and lambda-phage DNA) and revealed that complete dissociation occurs at 750 mM NaCl. Likely contributions of macromolecular desolvation and DNA flexibility to the binding energy are discussed.     

Molecular dynamics simulations of B-DNA: An analysis of the role of initial molecular configuration,randomly assigned velocity distribution,long integration times,and nonconstrained termini

Molecular dynamics simulations of three DNA sequences using the AMBER 3.0 force field were performed with implicit inclusion of water through a distance-dependent dielectric constant and solvated counterions. Simulations of the self-complementary DNA dodecamer d(CGCGAATTCGCG) were started from a regular B-DNA structure and the x-ray single crystal B-DNA structure. Although mean convergence during the 89-ps calculation was confirmed, localized differences in backbone torsionals and base-pair helicoidals were observed. A nanosecond simulation of the nonself-complementary 14 base-pair DNA d(GGCGGAATTGGCGG) indicates that most structural parameters stabilize within the first 100–200 ps, while isolated features show low-frequency oscillations throughout the calculation. The lack of harmonic constraints on the ends of the molecules was shown not to perturb the structural dynamics of the internal oligonucleotide beyond the external 2 base pairs. Comparison of three simulations of the nonself-complementary 14 base-pair DNA d(GGCGAAATTCGCGG), identical in all respects other than the assignment of initial Maxwellian atomic velocity distributions, revealed the inherent systematic variability. The three calculations result in nearly superimposable global structures, with localized variations in torsionals and helicoidals. Our results provide a basis for performing a comparative analysis of the effect of DNA sequence on localized structure. © 1993 John Wiley & Sons, Inc.     

Probing the role of structural water in a duplex oligodeoxyribonucleotide containing a water-mimicking base analog.

Molecular dynamics simulations were performed on models of the dodecamer DNA double-stranded segment, [d(CGCGAATTCGCG)](2), in which each of the adenine residues, individually or jointly, was replaced by the water-mimicking analog 2'-deoxy-7-(hydroxy-methyl)-7-deazaadenosine (hm(7)c(7)dA) [Rockhill, J.K., Wilson,S.R. and Gumport,R.I. (1996) J. Am. Chem. Soc.,118, 10065-10068]. The simulations, when compared with those of the dodecamer itself, show that incorporation of the analog affects neither the overall DNA structure nor its hydrogen-bonding and stacking interactions when it replaces a single individual base. Furthermore, the water molecules near the bases in the singly-substituted oligonucleotides are similarly unaffected. Double substitutions lead to differences in all the aforementioned parameters with respect to the reference sequence. The results suggest that the analog provides a good mimic of specific 'ordered' water molecules observed in contact with DNA itself and at the interface between protein and DNA in specific complexes.     

A Simple Solvation Model Along with A Multibody Dynamics Strategy (MBO(N)D) Produces Stable DNA Simulations that are Faster than Traditional Atomistic Methods

Abstract

We are developing solvation strategies that complement the speed advantage of MBO(N)D (a multibody simulation approach developed by Moldyn) for simulating biomolecular systems. In this report we propose to approximate the effect of bulk waters on DNA by using only a thin layer of waters proximate to the surface of DNA (which we will call the ‘thin shell approach’ or TSA). We will show that the TSA combined with substructuring (the grouping of atoms into rigid or flexible bodies) of the Dickerson dodecamer produces good comparisons with standard atomistic methods (over a nanosecond trajectory) as judged by a variety of DNA specific geometric (e.g., CURVES output) and dynamics (power spectra) properties. The MBO(N)D method, however, was faster than atomistic by a factor of six using the same solvation strategy and factor of 70 when compared to fully solvated atomistic system. The key to the speed of MBO(N)D is in its ability to use large time steps during dynamics. By keeping only a shell of molecules of water proximate to the dodecamer, we limit artifacts due to surface tension at the water-vacuum interface. These proximate waters are fairly immobile as compared to those in bulk and therefore do not severely limit the time step in the simulation. The strengths and limitations of this solvation approach, and future directions, will also be discussed.     

Backbone dynamics of DNA containing 8-oxoguanine: importance for substrate recognition by base excision repair glycosylases.

Except for the functional groups sited within the major or minor grooves, the bases of B-DNA are quite protected from the external environment. Enzymes that modify the bases often "flip out" the target into an extrahelical position before the chemistry step is carried out. Examples of this mechanism are the base excision repair glycosylases and the restriction enzyme methylases. The question arises about the mechanism of substrate recognition for these enzymes and how closely it is linked to the base flipping step. Molecular dynamics simulations (AMBER, PME electrostatics) of fully solvated, cation neutralized, DNA sequences containing 8-oxoguanine (8OG) and of appropriate normal (control) DNAs have been carried out. The dynamics trajectories were analyzed to identify those properties of the DNA structure in the vicinity of the altered base, or its dynamics, that could contribute to molecular discrimination between substrate and non-substrate DNA sites. The results predict that the FPG enzyme should flip out the cytosine base paired with the scissile 8OG, not the target base itself.     

Restrained molecular dynamics of solvated duplex DNA using the particle mesh Ewald method

Restrained and unrestrained aqueous solution molecular dynamics simulations applying the particle mesh Ewald (PME) method to DNA duplex structures previously determined via in vacuo restrained molecular dynamics with NMR-derived restraints are reported. Without experimental restraints, the DNA decamer, d(CATTTGCATC)d(GATGCAAATG) and trisdecamer, d(AGCTTGCCTTGAG)d(CTCAAGGCAAGCT), structures are stable on the nanosecond time scale and adopt conformations in the B-DNA family. These free DNA simulations exhibit behavior characteristic of PME simulations previously performed on DNA sequences, including a low helical twist, frequent sugar pucker transitions, BI- BII(–) transitions and coupled crankshaft (–) motion. Refinement protocols similar to the original in vacuo restrained molecular dynamics (RMD) refinements but in aqueous solution using the Cornell et al. force field [Cornell et al. (1995) J. Am. Chem. Soc., 117, 5179–5197] and a particle mesh Ewald treatment produce structures which fit the restraints very well and are very similar to the original in vacuo NMR structure, except for a significant difference in the average helical twist. Figures of merit for the average structure found in the RMD PME decamer simulations in solution are equivalent to the original in vacuo NMR structure while the figures of merit for the free MD simulations are significantly higher. The free MD simulations with the PME method, however, lead to some sequence-dependent structural features in common with the NMR structures, unlike free MD calculations with earlier force fields and protocols. There is some suggestion that the improved handling of electrostatics by PME improves long-range structural aspects which are not well defined by the short-range nature of NMR restraints.     

Molecular dynamics simulations of charged and neutral lipid bilayers: treatment of electrostatic interactions

Molecular dynamics (MD) simulations complement experimental methods in studies of the structure and dynamics of lipid bilayers. The choice of algorithms employed in this computational method represents a trade-off between the accuracy and real calculation time. The largest portion of the simulation time is devoted to calculation of long-range electrostatic interactions. To speed-up evaluation of these interactions, various approximations have been used. The most common ones are the truncation of long-range interactions with the use of cut-offs, and the particle-mesh Ewald (PME) method. In this study, several multi-nanosecond cut-off and PME simulations were performed to establish the influence of the simulation protocol on the bilayer properties. Two bilayers were used. One consisted of neutral phosphatidylcholine molecules. The other was a mixed lipid bilayer consisting of neutral phosphatidylethanolamine and negatively charged phosphatidylglycerol molecules. The study shows that the cut-off simulation of a bilayer containing charge molecules generates artefacts; in particular the mobility and order of the charged molecules are vastly different from those determined experimentally. In the PME simulation, the bilayer properties are in general agreement with experimental data. The cut-off simulation of bilayers containing only uncharged molecules does not generate artefacts, nevertheless, the PME simulation gives generally better agreement with experimental data.     

Bending and flexibility of methylated and unmethylated EcoRI DNA

We used cyclization kinetics experiments and Monte Carlo simulations to determine a structural model for a DNA decamer containing the EcoRI restriction site. Our findings agree well with recent crystal and NMR structures of the EcoRI dodecamer, where an overall bend of seven degrees is distributed symmetrically over the molecule. Monte Carlo simulations indicate that the sequence has a higher flexibility, assumed to be isotropic, compared to that of a "generic" DNA sequence. This model was used as a starting point for the investigation of the effect of cytosine methylation on DNA bending and flexibility. While methylation did not affect bend magnitude or direction, it resulted in a reduction in bending flexibility and under-winding of the methylated nucleotides. We demonstrate that our approach can augment the understanding of DNA structure and dynamics by adding information about the global structure and flexibility of the sequence. We also show that cyclization kinetics can be used to study the properties of modified nucleotides.     

Structure prediction and molecular simulation of gases diffusion pathways in hydrogenase

Although hydrogen is considered to be one of the most promising future energy sources and the technical aspects involved in using it have advanced considerably, the future supply of hydrogen from renewable sources is still unsolved. The [Fe]- hydrogenase enzymes are highly efficient H(2) catalysts found in ecologically and phylogenetically diverse microorganisms, including the photosynthetic green alga, Chlamydomonas reinhardtii. While these enzymes can occur in several forms, H(2) catalysis takes place at a unique [FeS] prosthetic group or H-cluster, located at the active site. 3D structure of the protein hydA1 hydrogenase from Chlamydomonas reinhardtti was predicted using the MODELER 8v2 software. Conserved region was depicted from the NCBI CDD Search. Template selection was done on the basis NCBI BLAST results. For single template 1FEH was used and for multiple templates 1FEH and 1HFE were used. The result of the Homology modeling was verified by uploading the file to SAVS server. On the basis of the SAVS result 3D structure predicted using single template was chosen for performing molecular simulation. For performing molecular simulation three strategies were used. First the molecular simulation of the protein was performed in solvated box containing bulk water. Then 100 H(2) molecules were randomly inserted in the solvated box and two simulations of 50 and 100 ps were performed. Similarly 100 O(2) molecules were randomly placed in the solvated box and again 50 and 100 ps simulation were performed. Energy minimization was performed before each simulation was performed. Conformations were saved after each simulation. Analysis of the gas diffusion was done on the basis of RMSD, Radius of Gyration and no. of gas molecule/ps plot.     

Sliding Box Docking: a new stand-alone tool for managing docking-based virtual screening along the DNA helix axis

Sliding Box Docking is a program that manages simulations of ligand docking at different defined positions of a three-dimensional DNA structure. The procedure is similar to inverse docking, which is a method that performs docking simulations of a single ligand in the active sites of different targets. Sliding Box Docking manages docking simulations of one ligand into a box that slides along the DNA helix axis in regular steps. For each box position a score is calculated using the separate Autodock Vina software, and the results are automatically plotted. The evaluation of ligand interaction at different DNA locations can highlight the specificity of ligands for different DNA- sequences. When assessing the affinity between ligans AT base pairs, results for docking simulations with a test set that included berenil, distamycin, hoechst 33258, and netropsin were as expected, agreeing well with affinities previously described in the literature.

Availability

Binaries are freely available at https://sourceforge.net/projects/slidingboxdocki     

Open-Boundary Molecular Dynamics of a DNA Molecule in a Hybrid Explicit/Implicit Salt Solution

The composition and electrolyte concentration of the aqueous bathing environment have important consequences for many biological processes and can profoundly affect the behavior of biomolecules. Nevertheless, because of computational limitations, many molecular simulations of biophysical systems can be performed only at specific ionic conditions: either at nominally zero salt concentration, i.e., including only counterions enforcing the system’s electroneutrality, or at excessive salt concentrations. Here, we introduce an efficient molecular dynamics simulation approach for an atomistic DNA molecule at realistic physiological ionic conditions. The simulations are performed by employing the open-boundary molecular dynamics method that allows for simulation of open systems that can exchange mass and linear momentum with the environment. In our open-boundary molecular dynamics approach, the computational burden is drastically alleviated by embedding the DNA molecule in a mixed explicit/implicit salt-bathing solution. In the explicit domain, the water molecules and ions are both overtly present in the system, whereas in the implicit water domain, only the ions are explicitly present and the water is described as a continuous dielectric medium. Water molecules are inserted and deleted into/from the system in the intermediate buffer domain that acts as a water reservoir to the explicit domain, with both water molecules and ions free to enter or leave the explicit domain. Our approach is general and allows for efficient molecular simulations of biomolecules solvated in bathing salt solutions at any ionic strength condition.     

Binding of a Hoechst dye to d(CGCGATATCGCG) and its influence on the conformation of the DNA fragment

Hoechst dye 33258 is a planar drug molecule that binds to the minor groove of DNA, especially where there are a number of A.T base pairs. We have solved the structure of the Hoechst dye bound to the DNA dodecamer d(CGCGATATCGCG) at 2.3 A. This structure is compared to that of the same dodecamer with the minor-groove-binding drug netropsin bound to it, as well as to structures that have been solved for this Hoechst dye bound to a DNA dodecamer containing the central four base pairs with the sequence AATT. We find that the position of the Hoechst drug in this dodecamer is quite different from that found in the other dodecamer since it has an opposite orientation compared to the other two structures. The drug covers three of the four A.T base pairs and extends its piperazine ring to the first G.C base pair adjacent to the alternating AT segment. Furthermore, the drug binding has modified the structure of the DNA dodecamer. Other DNA dodecamers with alternating AT sequences show an alternation in the size of the helical twist between the ApT step (small twist) and the TpA step (large twist). In this structure the alternation is reversed with larger twists in the ApT steps than in the TpA step. In addition, there is a rotation of one of the thymine bases in the DNA dodecamer that is associated with hydrogen bonding to the Hoechst drug. This structure illustrates the considerable plasticity found in the DNA molecule when it binds to different planar molecules inserted into the minor groove.     

Cell cycle dependent morphology changes and associated mitochondrial DNA redistribution in mitochondria of human cell lines

Mitochondria of osteosarcoma cells (143B) in culture have variable morphologies, classified according to the shape and size of the organelle as reticular, fragmented or intermediate. Synchronization and release from G0 has shown that the morphology of mitochondria oscillates between the reticular and fragmented state in a cell cycle dependent manner. Cells in G1 have reticular mitochondria while those in S phase have fragmented mitochondria. By using a novel method of fluorescence in situ hybridization, the morphology of mitochondria was correlated with mitochondrial DNA distribution. MtDNA molecules were seen in clusters of two to four along mitochondrial filaments. In the fully fragmented state, each mitochondrion contained at least one cluster. We discuss the importance of fission and fusion events in regulating the morphology of mitochondria, segregation of mtDNA and maintenance of the organelle's functional unity.     

Structural studies on the second Mycobacterium smegmatis Dps: invariant and variable features of structure, assembly and function

A second DNA binding protein from stationary-phase cells of Mycobacterium smegmatis (MsDps2) has been identified from the bacterial genome. It was cloned, expressed and characterised and its crystal structure was determined. The core dodecameric structure of MsDps2 is the same as that of the Dps from the organism described earlier (MsDps1). However, MsDps2 possesses a long N-terminal tail instead of the C-terminal tail in MsDps1. This tail appears to be involved in DNA binding. It is also intimately involved in stabilizing the dodecamer. Partly on account of this factor, MsDps2 assembles straightway into the dodecamer, while MsDps1 does so on incubation after going through an intermediate trimeric stage. The ferroxidation centre is similar in the two proteins, while the pores leading to it exhibit some difference. The mode of sequestration of DNA in the crystalline array of molecules, as evidenced by the crystal structures, appears to be different in MsDps1 and MsDps2, highlighting the variability in the mode of Dps-DNA complexation. A sequence search led to the identification of 300 Dps molecules in bacteria with known genome sequences. Fifty bacteria contain two or more types of Dps molecules each, while 195 contain only one type. Some bacteria, notably some pathogenic ones, do not contain Dps. A sequence signature for Dps could also be derived from the analysis.     

Pentacoordinate oxyphosphorane intermediate always exists in aqueous solution.

Gas-phase ab initio calculations indicate that dianionic pentacoordinate oxyphosphoranes do not have a kinetically meaningful intermediate. The simplest oxyphosphorane PO5H3(2-) has the least tendency to have a pentacoordinate intermediate. However, it does have a pentacoordinate intermediate when it is solvated with six water molecules. These results support the hypothesis that the phosphoryl transfer reactions take place via pentacoordinate intermediate not only in acidic but also in basic media.     

Mediation of the A/B-DNA helix transition by G-tracts in the crystal structure of duplex CATGGGCCCATG

The crystal structure of the DNA dodecamer duplex CATGGGCCCATG lies on a structural continuum along the transition between A- and B-DNA. The dodecamer possesses the normal vector plot and inclination values typical of B-DNA, but has the crystal packing, helical twist, groove width, sugar pucker, slide and x-displacement values typical of A-DNA. The structure shows highly ordered water structures, such as a double spine of water molecules against each side of the major groove, stabilizing the GC base pairs in an A-like conformation. The different hydration of GC and AT base pairs provides a physical basis for solvent-dependent facilitation of the A↔B helix transition by GC base pairs. Crystal structures of CATGGGCCCATG and other A/B-DNA intermediates support a ‘slide first, roll later’ mechanism for the B→A helix transition. In the distribution of helical parameters in protein–DNA crystal structures, GpG base steps show A-like properties, reflecting their innate predisposition for the A conformation.     

Macromolecular design, nucleic acid junctions, and crystal formation

The simplest form of macromolecular design involves the ligation of nucleic acids. Recent results on the concatenation of nucleic acid junctions show that these molecules can act as fairly rigid macromolecular valence clusters on the nanometer scale. These clusters can be joined to form closed stick figures in which each edge is double helical DNA or RNA and each vertex is a nucleic acid junction. The geometrical criteria for forming discrete-closed and periodic structures from these components are established. The helicity of each edge limits the possible structures that can be formed. The formation of a periodic array from nucleic acid junction building blocks is compared with the crystallization of molecular systems. This comparison leads to a new interpretation of the nature of order in the solid state for molecular crystals. The suggestion is made that the structure of a solid molecular system described by the fewest unique orthogonal (Fourier) components is the one which will be entropically favored, since it contains the least information. This is the crystalline state, with a small number of molecules per asymmetric unit. The free energy from the proposed entropic driving force responsible for this behavior is available, in principle, to correct small deviations from ideality in forming covalent crystals from nucleic acid junction components, as well as in non-bonded molecular systems. Nucleic acid junction periodic arrays provide an appropriate vehicle with which to test this interpretation.     

Imaging of kinked configurations of DNA molecules undergoing orthogonal field alternating gel electrophoresis by fluorescence microscopy

The dynamics of individual DNA molecules undergoing orthogonal field alternating gel electrophoresis (OFAGE) have been studied by use of T2 DNA molecules labeled with a dye and visualized with a fluorescence microscope. The mechanism of reorientation used by a molecule to align itself in the direction of the new orthogonal field depends on the degree of extension of the chain immediately before the application of this field. The formation of kinks is promoted when time is allowed between the application of the two orthogonal fields so that the molecule attains a partially relaxed configuration. In this case, the chain appears bunched up in domains moving along the contour of the molecule. These regions are found to be the locations where the kinks are formed upon application of the second field perpendicular to the chain. The formation of kinks provides a significant retardation of the reorientation of the molecules, relative to molecules that do not form kinks, and appears to play an important role in the fractionation attained with OFAGE. A classification of various reorientation mechanisms observed in molecules that form kinks is presented.     

Validation of the 53A6 GROMOS force field

The quality of biomolecular dynamics simulations relies critically on the force field that is used to describe the interactions between particles in the system. Force fields, which are generally parameterized using experimental data on small molecules, can only prove themselves in realistic simulations of relevant biomolecular systems. In this work, we begin the validation of the new 53A6 GROMOS parameter set by examining three test cases. Simulations of the well-studied 129 residue protein hen egg-white lysozyme, of the DNA dodecamer d(CGCGAATTCGCG)₂, and a proteinogenic ³-dodecapeptide were performed and analysed. It was found that the new parameter set performs as well as the previous parameter sets in terms of protein (45A3) and DNA (45A4) stability and that it is better at describing the folding–unfolding balance of the peptide. The latter is a property that is directly associated with the free enthalpy of hydration, to which the 53A6 parameter set was parameterized.