A molecular dynamics computer simulation of an eight-base-pair DNA fragment in aqueous solution: comparison with experimental two-dimensional NMR data

Molecular dynamics simulation or structure refinement of proteins: are solvent molecules required? A case study using hen lysozyme

In protein simulation or structure refinement based on values of observable quantities measured in (aqueous) solution, solvent (water) molecules may be explicitly treated, omitted, or represented by a potential of mean-solvation-force term, depending on protein coordinates only, in the force field used. These three approaches are compared for hen egg white lysozyme (HEWL). This 129-residue non-spherical protein contains a variety of secondary-structure elements, and ample experimental data are available: 1630 atom–atom Nuclear Overhauser Enhancement (NOE) upper distance bounds, 213 ³ J-couplings and 200 S² order parameters. These data are used to compare the performance of the three approaches. It is found that a molecular dynamics (MD) simulation in explicit water approximates the experimental data much better than stochastic dynamics (SD) simulation in vacuo without or with a solvent-accessible-surface-area (SASA) implicit-solvation term added to the force field. This is due to the missing energetic and entropic contributions and hydrogen-bonding capacities of the water molecules and the missing dielectric screening effect of this high-permittivity solvent. Omission of explicit water molecules leads to compaction of the protein, an increased internal strain, distortion of exposed loop and turn regions and excessive intra-protein hydrogen bonding. As a consequence, the conformation and dynamics of groups on the surface of the protein, which may play a key role in protein–protein interactions or ligand or substrate binding, may be incorrectly modelled. It is thus recommended to include water molecules explicitly in structure refinement of proteins in aqueous solution based on nuclear magnetic resonance (NMR) or other experimentally measured data.

On the use of intra-molecular distance and angle constraints to lengthen the time step in molecular and stochastic dynamics simulations of proteins

Computer simulation of proteins in aqueous solution at the atomic level of resolution is still limited in time span and system size due to limited computing power available and thus employs a variety of time‐saving techniques that trade some accuracy against computational effort. An example of such a time‐saving technique is the application of constraints to particular degrees of freedom when integrating Newton's or Langevin's equations of motion in molecular dynamics (MD) or stochastic dynamics (SD) simulations, respectively. The application of bond‐length constraints is standard practice in protein simulations and allows for a lengthening of the time step by a factor of three. Applying recently proposed algorithms to constrain bond angles or dihedral angles, it is investigated, using the protein trypsin inhibitor as test molecule, whether bond angles and dihedral angles involving hydrogen atoms or even stiff proper (torsional) dihedral angles as well as improper ones (maintaining particular tetrahedral or planar geometries) may be constrained without generating too many artificial side effects. Constraining the relative positions of the hydrogen atoms in the protein allows for a lengthening of the time step by a factor of two. Additionally constraining the improper dihedral angles and the stiff proper (torsional) dihedral angles in the protein does not allow for an increase of the MD or SD time step.

On the use of multiple-time-step algorithms to save computing effort in molecular dynamics simulations of proteins

Computer simulation of proteins in aqueous solution at the atomic level of resolution is still limited in time span and system size due to limited computing power available and thus employs a variety of time-saving techniques that trade some accuracy against computational effort. Examples of such time-saving techniques are the application of constraints to particular degrees of freedom or the use of a multiple-time-step (MTS) algorithm distinguishing between particular forces when integrating Newton's equations of motion. The application of two types of MTS algorithms to bond-stretching forces versus the remaining forces in molecular dynamics (MD) simulations of a protein in aqueous solution or of liquid water is investigated and the results in terms of total energy conservation and the influence on various other properties are compared to those of MD simulations of the same systems using bond-length, and for water bond-angle, constraints. At comparable computational effort, the use of bond-length constraints in proteins leads to better energy conservation and less distorted properties than the two MTS algorithms investigated.

On the use of 3J-coupling NMR data to derive structural information on proteins

Values of ³J-couplings as obtained from NMR experiments on proteins cannot easily be used to determine protein structure due to the difficulty of accounting for the high sensitivity of intermediate ³J-coupling values (4-8 Hz) to the averaging period that must cover the conformational variability of the torsional angle related to the ³J-coupling, and due to the difficulty of handling the multiple-valued character of the inverse Karplus relation between torsional angle and ³J-coupling. Both problems can be solved by using ³J-coupling time-averaging local-elevation restraining MD simulation. Application to the protein hen egg white lysozyme using 213 backbone and side-chain ³J-coupling restraints shows that a conformational ensemble compatible with the experimental data can be obtained using this technique, and that accounting for averaging and the ability of the algorithm to escape from local minima for the torsional angle induced by the Karplus relation, are essential for a comprehensive use of ³J-coupling data in protein structure determination.

On the Effect of the Various Assumptions and Approximations used in Molecular Simulations on the Properties of Bio-Molecular Systems: Overview and Perspective on Issues

Computer simulations of molecular systems enable structure-energy-function relationships of molecular processes to be described at the sub-atomic, atomic, supra-atomic or supra-molecular level and plays an increasingly important role in chemistry, biology and physics. To interpret the results of such simulations appropriately, the degree of uncertainty and potential errors affecting the calculated properties must be considered. Uncertainty and errors arise from (1) assumptions underlying the molecular model, force field and simulation algorithms, (2) approximations implicit in the interatomic interaction function (force field), or when integrating the equations of motion, (3) the chosen values of the parameters that determine the accuracy of the approximations used, and (4) the nature of the system and the property of interest. In this overview, advantages and shortcomings of assumptions and approximations commonly used when simulating bio-molecular systems are considered. What the developers of bio-molecular force fields and simulation software can do to facilitate and broaden research involving bio-molecular simulations is also discussed.

On the Use of Side-Chain NMR Relaxation Data to Derive Structural and Dynamical Information on Proteins: A Case Study Using Hen Lysozyme

Values of S2CH and S2NH order parameters derived from NMR relaxation measurements on proteins cannot be used straightforwardly to determine protein structure because they cannot be related to a single protein structure, but are defined in terms of an average over a conformational ensemble. Molecular dynamics simulation can generate a conformational ensemble and thus can be used to restrain S2CH and S2NH order parameters towards experimentally derived target values S2CH (exp) and S2NH (exp). Application of S2CH and S2NH order‐parameter restraining MD simulation to bond vectors in 63 side chains of the protein hen egg white lysozyme using 51 S2CH (exp) target values and 28 S2NH (exp) target values shows that a conformational ensemble compatible with the experimentally derived data can be obtained by using this technique. It is observed that S2CH order‐parameter restraining of C−H bonds in methyl groups is less reliable than S2NH order‐parameter restraining because of the possibly less valid assumptions and approximations used to derive experimental S2CH (exp) values from NMR relaxation measurements and the necessity to adopt the assumption of uniform rotational motion of methyl C−H bonds around their symmetry axis and of the independence of these motions from each other. The restrained simulations demonstrate that side chains on the protein surface are highly dynamic. Any hydrogen bonds they form and that appear in any of four different crystal structures, are fluctuating with short lifetimes in solution.

Dynamics of K+ counterions around DNA double helix in the external electric field: A molecular dynamics study

The structure of DNA double helix is stabilized by metal counterions condensed to a diffuse layer around the macromolecule. The dynamics of counterions in real conditions is governed by the electric fields from DNA and other biological macromolecules. In the present work the molecular dynamics study was performed for the system of DNA double helix with neutralizing K⁺ counterions and for the system of KCl salt solution in an external electric field of different strength (up to 32mV/Å). The analysis of ionic conductivities of these systems has shown that the counterions around the DNA double helix are slowed down compared with the KCl salt solution. The calculated values of ion mobility are within (0.05-0.4)mS/cm depending on the orientation of the external electric field relatively to the double helix. Under the electric field parallel to the macromolecule K⁺ counterions move along the grooves of the double helix staying longer in the places with narrower minor groove. Under the electric field perpendicular to the macromolecule the dynamics of counterions is less affected by DNA atoms, and starting with the electric field values about 30mV/Å the double helix undergoes a phase transition from a double-stranded to a single-strand state.

Density artefacts at interfaces caused by multiple time-step effects in molecular dynamics simulations

Background: Molecular dynamics (MD) simulations have become an important tool to provide insight into molecular processes involving biomolecules such as proteins, DNA, carbohydrates and membranes. As these processes cover a wide range of time scales, multiple time-step integration methods are often employed to increase the speed of MD simulations. For example, in the twin-range (TR) scheme, the nonbonded forces within the long-range cutoff are split into a short-range contribution updated every time step (inner time step) and a less frequently updated mid-range contribution (outer time step). The presence of different time steps can, however, cause numerical artefacts. Methods: The effects of multiple time-step algorithms at interfaces between polar and apolar media are investigated with MD simulations. Such interfaces occur with biological membranes or proteins in solution. Results: In this work, it is shown that the TR splitting of the nonbonded forces leads to artificial density increases at interfaces. The presence of the observed artefacts was found to be independent of the interface shape and the thermostatting method used. It is further shown that integration with an impulse-wise reversible reference system propagation algorithm (RESPA) only shifts the occurrence of density artefacts towards larger outer time steps. Using a single-range (SR) treatment of the nonbonded interactions, on the other hand, resolves the density issue for pairlist-update periods of up to 40 fs. Conclusion: A SR scheme avoids numerical artefacts and offers an interesting alternative to TR RESPA with respect to performance optimization.

Density artefacts at interfaces caused by multiple time-step effects in molecular dynamics simulations

Background: Molecular dynamics (MD) simulations have become an important tool to provide insight into molecular processes involving biomolecules such as proteins, DNA, carbohydrates and membranes. As these processes cover a wide range of time scales, multiple time-step integration methods are often employed to increase the speed of MD simulations. For example, in the twin-range (TR) scheme, the nonbonded forces within the long-range cutoff are split into a short-range contribution updated every time step (inner time step) and a less frequently updated mid-range contribution (outer time step). The presence of different time steps can, however, cause numerical artefacts. Methods: The effects of multiple time-step algorithms at interfaces between polar and apolar media are investigated with MD simulations. Such interfaces occur with biological membranes or proteins in solution. Results: In this work, it is shown that the TR splitting of the nonbonded forces leads to artificial density increases at interfaces for weak coupling and Nosé-Hoover (chain) thermostats. It is further shown that integration with an impulse-wise reversible reference system propagation algorithm (RESPA) only shifts the occurrence of density artefacts towards larger outer time steps. Using a single-range (SR) treatment of the nonbonded interactions or a stochastic dynamics thermostat, on the other hand, resolves the density issue for pairlist-update periods of up to 40 fs. Conclusion: TR schemes are not advisable to use in combination with weak coupling or Nosé-Hoover (chain) thermostats due to the occurrence of significant numerical artifacts at interfaces.

Density artefacts at interfaces caused by multiple time-step effects in molecular dynamics simulations

Background: Molecular dynamics (MD) simulations have become an important tool to provide insight into molecular processes involving biomolecules such as proteins, DNA, carbohydrates and membranes. As these processes cover a wide range of time scales, multiple time-step integration methods are often employed to increase the speed of MD simulations. For example, in the twin-range (TR) scheme, the nonbonded forces within the long-range cutoff are split into a short-range contribution updated every time step (inner time step) and a less frequently updated mid-range contribution (outer time step). The presence of different time steps can, however, cause numerical artefacts. Methods: The effects of multiple time-step algorithms at interfaces between polar and apolar media are investigated with MD simulations. Such interfaces occur with biological membranes or proteins in solution. Results: In this work, it is shown that the TR splitting of the nonbonded forces leads to artificial density increases at interfaces for weak coupling and Nosé-Hoover (chain) thermostats. It is further shown that integration with an impulse-wise reversible reference system propagation algorithm (RESPA) only shifts the occurrence of density artefacts towards larger outer time steps. Using a single-range (SR) treatment of the nonbonded interactions or a stochastic dynamics thermostat, on the other hand, resolves the density issue for pairlist-update periods of up to 40 fs. Conclusion: TR schemes are not advisable to use in combination with weak coupling or Nosé-Hoover (chain) thermostats due to the occurrence of significant numerical artifacts at interfaces.

Using Complementary NMR Data Sets To Detect Inconsistencies and Model Flaws in the Structure Determination of Human Interleukin-4

The derivation of protein structure from values of observable quantities measured in NMR experiments is a rather nontrivial task due to (i) the limited number of data compared to degrees of freedom of a protein, (ii) the uncertainty inherent to the function connecting an observable quantity to molecular structure, (iii) the finite quality of biomolecular models and force fields used in structure refinement, and (iv) the conformational freedom of a protein in aqueous solution, which requires extensive conformational sampling and appropriate conformational averaging when calculating or restraining to sets of NMR data. The protein interleukin-4 (IL-4) has been taken as a test case using NOE distances, S² order parameters, and ³J-couplings as test data and the former two types of data as restraints. It is shown that, by combining sets of different, complementary NMR data as restraints in MD simulations, inconsistencies in the data or flaws in the model and procedures used to derive protein structure from NMR data can be detected. This leads to an improved structural interpretation of such data particularly in more mobile loop regions.

Challenges in computational studies of enzyme structure, function and dynamics

In this review we give an overview of the field of Computational enzymology. We start by describing the birth of the field, with emphasis on the work of the 2013 chemistry Nobel Laureates. We then present key features of the state-of-the-art in the field, showing what theory, accompanied by experiments, has taught us so far about enzymes. We also briefly describe computational methods, such as quantum mechanics-molecular mechanics approaches, reaction coordinate treatment, and free energy simulation approaches. We finalize by discussing open questions and challenges.

The long and winding road to the structure of homo-DNA

Homo-DNA ((4'-->6')-linked oligo-2',3'-dideoxy-beta-D-glucopyranose nucleic acid) constitutes the earliest synthetic model system whose pairing properties have been studied within an etiology of nucleic acid structure. Its conception as part of a program directed at a rationalization of Nature's selection of pentoses over other candidates as the carbohydrate building block in the genetic material was motivated by the question: why pentose and not hexose? Homo-DNA forms an autonomous pairing system and its duplexes are entropically stabilized relative to DNA duplexes. Moreover, the base pairing priorities in homo-DNA duplexes differ from those in DNA. A deeper understanding of the particular properties of homo-DNA requires knowledge of its structure. Although diffraction data for crystals of a homo-DNA octamer duplex were available to medium resolution in the mid-1990s, it took another decade for the structure to be solved. In this tutorial Review we describe the odyssey from the crystallization to the final structure determination with its many failures and disappointments and the development of selenium chemistry to derivatize nucleic acids for crystallographic phasing. More than fifty years after the discovery of the DNA double helix, the story of homo-DNA also provides a demonstration of the limits of theoretical models and offers a fresh view of fundamental issues in regard to the natural nucleic acids, such as the origins of antiparallel pairing and helicality.

Comparing atomistic simulation data with the NMR experiment: how much can NOEs actually tell us?

Simulated molecular dynamics trajectories of proteins and nucleic acids are often compared with nuclear magnetic resonance (NMR) data for the purposes of assessing the quality of the force field used or, equally important, trying to interpret ambiguous experimental data. In particular, nuclear Overhauser enhancement (NOE) intensities or atom-atom distances derived from them are frequently calculated from the simulated ensembles because the distance restraints derived from NOEs are the key ingredient in NMR-based protein structure determination. In this study, we ask how diverse and nonnative-like an ensemble of structures can be and still match the experimental NOE distance upper bounds well. We present two examples in which simulated ensembles of highly nonnative polypeptide structures (an unfolded state ensemble of the villin headpiece and a high-temperature denatured ensemble of lysozyme) are shown to match fairly well the experimental NOE distance upper bounds from which the corresponding native structures were derived. For example, the unfolded ensemble of villin headpiece, which is on average 0.90 +/- 0.13 nm root-mean-square deviation away from the native NMR structure, deviates from the experimental restraints by only 0.027 nm on average. However, this artificially good agreement is largely a consequence of 1) the highly nonlinear effects of r(-6) (or r(-3)) averaging and 2) focusing only on the experimentally observed set of NOE bounds. Namely, in addition to the experimentally observed NOEs, both simulated ensembles (especially the villin ensemble) also predict a large number of NOEs, which are not seen in the experiment. If these are taken into account, the agreement between simulation and experiment gets markedly worse, as it should, given the nonnative nature of the underlying simulated ensembles. In light of the examples given, we conclude that comparing experimental NOE distance restraints with large simulated ensembles provides just by itself only limited information about the quality of simulation.

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 dynamics in aqueous solution: opening the double helix

The opening of a DNA base pair is a simple reaction that is a prerequisite for replication, transcription, and other vital biological functions. Understanding the molecular mechanisms of biological reactions is crucial for predicting and, ultimately, controlling them. Realistic computer simulations of the reactions can provide the needed understanding. To model even the simplest reaction in aqueous solution requires hundreds of hours of supercomputing time. We have used molecular dynamics techniques to simulate fraying of the ends of a six base pair double strand of DNA, [TCGCGA]2, where the four bases of DNA are denoted by T (thymine), C (cytosine), G (guanine), and A (adenine), and to estimate the free energy barrier to this process. The calculations, in which the DNA was surrounded by 2,594 water molecules, required 50 hours of CRAY-2 CPU time for every simulated 100 picoseconds. A free energy barrier to fraying, which is mainly characterized by the movement of adenine away from thymine into aqueous environment, was estimated to be 4 kcal/mol. Another fraying pathway, which leads to stacking between terminal adenine and thymine, was also observed. These detailed pictures of the motions and energetics of DNA base pair opening in water are a first step toward understanding how DNA will interact with any molecule.

Molecular dynamics simulation of nucleic acids

We review molecular dynamics simulations of nucleic acids, including those completed from 1995 to 2000, with a focus on the applications and results rather than the methods. After the introduction, which discusses recent advances in the simulation of nucleic acids in solution, we describe force fields for nucleic acids and then provide a detailed summary of the published literature. We emphasize simulations of small nucleic acids ( approximately 6 to 24 mer) in explicit solvent with counterions, using reliable force fields and modern simulation protocols that properly represent the long-range electrostatic interactions. We also provide some limited discussion of simulation in the absence of explicit solvent. Absent from this discussion are results from simulations of protein-nucleic acid complexes and modified DNA analogs. Highlights from the molecular dynamics simulation are the spontaneous observation of A B transitions in duplex DNA in response to the environment, specific ion binding and hydration, and reliable representation of protein-nucleic acid interactions. We close by examining major issues and the future promise for these methods.

DNA Interstrand Crosslink Formation by Mechlorethamine at a Cytosine–Cytosine Mismatch Pair: Kinetics and Sequence Dependence

Expansion of the triplet repeat DNA sequence d[CGG]n.d[CCG]n is a characteristic of Fragile X syndrome, a human neurodegenerative disease. Stable intrastrand conformations formed by both d[CGG]n and d[CCG]n, and involving G-G and C-C mismatch pairs, respectively, are believed to be of importance in the development of the disease. We have shown previously that C-C mismatch pairs can be crosslinked covalently by mechlorethamine, a nitrogen mustard alkylating agent, and hence this reaction may be of value as a probe for conformers of d[CCG]n. To characterize the mechlorethamine C-C crosslink reaction further, here we report the kinetics and sequence dependence of formation of the crosslink species, using a series of model duplexes. The rate of reaction depends on the base sequence proximal to the C-C mismatch pair. Hence, in 19mer duplexes containing a central d[M4M3M2M1Cn1n2n3n4].d[N4N3N2N1Cm1m2m3m4] sequence, where M-m and N-n are complementary base pairs, the amount of crosslink increased with increasing G-C content of the eight base pairs neighboring the C-C mismatch and with the proximity of the G-C pairs to the C-C mismatch. Molecular dynamics simulations of the solvated duplexes provided an explanation of these data. Hence, for a C-C pair flanked by G-C base pairs the mismatched cytosine bases remain stacked within the duplex, but for a C-C pair flanked by A-T base pairs, the simulations suggested local opening of the duplex around the C-C pair, making it a less effective target for mechlorethamine.

Molecular dynamics simulations of DNA in solutions with different counter-ions

Molecular dynamics simulations of the [d(ATGCAGTCAG]2 fragment of DNA, in water and in the presence of three different counter-ions (Li+, Na+ and Cs+) are reported. Three-dimensional hydration structure and ion distribution have been calculated using spatial distribution functions for a detailed picture of local concentrations of ions and water molecules around DNA. According to the simulations, Cs+ ions bind directly to the bases in the minor groove, Na+ ions bind prevailing to the bases in the minor groove through one water molecule, whereas Li+ ions bind directly to the phosphate oxygens. The different behavior of the counter-ions is explained by specific hydration structures around the DNA and the ions. It is proposed how the observed differences in the ion binding to DNA may explain different conformational behavior of DNA. Calculated self-diffusion coefficients for the ions agree well with the available NMR data.

DNA structure and fluctuations sensed from a 1.1ns molecular dynamics trajectory of a fully charged Zif268-DNA complex in water

Molecular dynamics simulations of the zinc finger domain of protein Zif268, in a complex with a high affinity DNA sequence, yields a globally stable system with small yet significant readjustments with persistence time of the order of 1.1ns. The results confirm the quality of the standard GROMOS87 force field with a corrected solvent-to-solute interaction that does not affect the water-water SPC interactions nor the intra-molecular cohesive forces. Specificity determinants are discussed. The simulations of DNA alone, with the same force field, showed the important role played by the solvent and the symmetry of the counterion distribution. (Tapia & Velázquez, J. Am. Chem. Soc., 119, 5934, 1997) In the present work, this feature was retained when appropriate. The results for root mean square deviations and temperature B-factors illustrate the reliability of this approach. The structure of DNA is held by its interactions with the zinc finger protein. This behavior is not much affected by the slow whithering away of finger-1 from DNA. The factors contributing to the molecular stability found in GROMOS' potential energy function appear to be sufficient to yield stable fluctuation patterns when surrounding medium effects are properly included.

The intrinsic curvature of a 51 bp K-DNA fragment of Leishmania tarentolae: a molecular model

DNA intrinsic structure and curvature is a subject of debate because of the importance of these attributes in processes such as DNA packaging, transcription, and gene regulation. X-ray crystallography of DNA single crystals has provided a wealth of information about the local, short range conformational features of DNA. On the other hand, gel electrophoresis analysis of DNA has not only uncovered the macroscopic curvature of DNA but it also provides most of the available data on DNA intrinsic curvature. However, gel electrophoresis can not identify features of DNA structure at the nucleotide or atomic level. In order to address the problem of DNA intrinsic curvature in an attempt to bridge the gap between X-ray crystallography and gel electrophoresis, we use the computational method of molecular dynamics (MD). In this study, we report the results of 2.0 ns MD simulations on a 51 bp fragment of the K-DNA of Leishmania tarentolae containing several A-tracts. The K-DNA double helix is very stable and remains in an intermediate state between the canonical A and B forms of the duplex. The magnitude of global curvature (75 degrees) agrees well with the experimental estimate (72 degrees) available. Analysis of local (every base triplet) and sublocal (every helix turn) curvature shows that the 51 bp K-DNA fragment has curvature features also present in the Wedge, Junction and Calladine's models of DNA intrinsic curvature. We further characterize the flexibility of individual nucleotides in the molecule and find the sugar flexibility within the A-tracts to be strongly correlated with the pattern of A-tract cleavage by the hydroxyl radical. Differential curvature and flexibility at the 5' and 3'junctions between A-tracts and general-sequence DNA are found to modulate the global curvature of the K-DNA fragment.

The trinucleotide repeat sequence d(GTC)15 adopts a hairpin conformation

The structure of a single-stranded (ss) oligonucleotide containing (GTC)15 [ss(GTC)15] was examined. As a control, parallel studies were performed with ss(CTG)15, an oligonucleotide that forms a hairpin. Electrophoretic mobility, KMnO4 oxidation and P1 nuclease studies demonstrate that, similar to ss(CTG)15, ss(GTC)15 forms a hairpin containing base paired and/or stacked thymines in the stem. Electrophoretic mobility melting profiles performed in approximately 1 mM Na+ revealed that the melting temperature of ss(GTC)15 and ss(CTG)15 were 38 and 48 degrees C respectively. The loop regions of ss(GTC)15 and ss(CTG)15 were cleaved by single-strand-specific P1 nuclease at the T25-C29 and G26-C27 phosphodiester bonds respectively (where the loop apex of the DNAs is T28). Molecular dynamics simulations suggested that in ss(GTC)15 the loop was bent towards the major groove of the stem, apparently causing an increased exposure of the T25-C29 region to solvent. In ss(CTG)15 guanine--guanine stacking caused a separation of the G26 and C27 bases, resulting in exposure of the intervening phosphodiester to solvent. The results suggest that ss(GTC)15 and ss(CTG)15 form similar, but distinguishable, hairpin structures.

Molecular dynamics simulations of the glucocorticoid receptor DNA-binding domain in complex with DNA and free in solution

Molecular dynamics simulations have been performed on the glucocorticoid receptor DNA binding domain (GR DBD) in aqueous solution as a dimer in complex with DNA and as a free monomer. In the simulated complex, we find a slightly increased bending of the DNA helix axis compared with the crystal structure in the spacer region of DNA between the two half-sites that are recognized by GR DBD. The bend is mainly caused by an increased number of interactions between DNA and the N-terminal extended region of the sequence specifically bound monomer. The recognition helices of GR DBD are pulled further into the DNA major groove leading to a weakening of the intrahelical hydrogen bonds in the middle of the helices. Many ordered water molecules with long residence times are found at the intermolecular interfaces of the complex. The hydrogen-bonding networks (including water bridges) on either side of the DNA major groove involve residues that are highly conserved within the family of nuclear receptors. Very similar hydrogen-bonding networks are found in the estrogen receptor (ER) DBD in complex with DNA, which suggests that this is a common feature for proper positioning of the recognition helix in ER DBD and GR DBD.

Molecular dynamics simulations of chlorambucil/DNA adducts. A structural basis for the 5'-GNC interstrand DNA crosslink formed by nitrogen mustards

The alkylation of DNA by chlorambucil has been studied using a computational approach. Molecular dynamics simulations were performed on the fully solvated non-covalent complex, two monoadducts and a crosslinked diadduct of chlorambucil with the d(CGG3G2CGC).-d(GCG1CCCG) duplex, in which the N7 atoms of G1, G2 and G3 are potential alkylation sites. The results provide a structural basis for the preference of nitrogen mustards to crosslink DNA duplexes at a 5'-GNC site (a 1,3 crosslink, G1-G3) rather than at a 5'-GC sites (a 1,2 crosslink, G1-G2). In the non-covalent complex simulation the drug reoriented from a non-interstrand crosslinking location to a position favorable for G1-G3 diadduct formation. It proved possible to construct a G1-G3 diadduct from a structure from the non-covalent simulation, and continue the molecular dynamics calculation without further disruption of the DNA structure. A crosslinked diadduct developed with four BII conformations on the 3' side of each alkylated guanine and of their respective complementary cytosine. In the first monoadduct simulation the starting point was the same DNA conformation used in the crosslinked diadduct simulation with alkylation at G1. In this simulation the DNA deformation was reduced, with the helix returning to a more canonical form. A second monoadduct simulation was started from a canonical DNA conformation alkylated at G3. Here, no significant motion towards a potential crosslinking conformation occurred. Collectively, the results suggest that crosslink formation is dependent upon the drug orientation prior to alkylation and the required deformation of the DNA to permit 1,3 crosslinking can largely be achieved in the non-covalent complex.

Grand canonical Monte Carlo molecular and thermodynamic predictions of ion effects on binding of an oligocation (L8+) to the center of DNA oligomers

Grand canonical Monte Carlo (GCMC) simulations are reported for aqueous solutions containing excess univalent salt (activities a +/- = 1.76-12.3 mM) and one of the following species: an octacationic rod-like ligand, L8+; a B-DNA oligomer with N phosphate charges (8 < or = N < or = 100); or a complex resulting from the binding of L8+ at the center of an N-mer (24 < or = N < or = 250). Simplified models of these multiply charged species are used in the GCMC simulations to predict the fundamental coulombic contributions to the following experimentally relevant properties: 1) the axial distance over which ligand binding affects local counterion concentrations at the surface of the N-mer; 2) the dependence on N of GCMC preferential interaction coefficients, gamma 32MC identical to delta C3/delta C2l a +/-, T, where C3 and C2 are, respectively, the molar concentrations of salt and the multiply charged species (ligand, N-mer or complex); and 3) the dependence on N of SaKobs identical to d in Kobs/d in a +/- = delta (magnitude of ZJ + 2 gamma 32J), where Kobs is the equilibrium concentration quotient for the binding of L8+ to the center of an N-mer and delta denotes the stoichiometric combination of terms, each of which pertains to a reactant or product J having magnitude of ZJ charges. The participation of electrolyte ions in the ligand binding interaction is quantified by the magnitude of SaKobs, which reflects the net (stoichiometrically weighted) difference in the extent of thermodynamic binding of salt ions to the products and reactants. Results obtained here from GCMC simulations yield a picture of the salient molecular consequences of binding a cationic ligand, as well as thermodynamic predictions whose applicability can be tested experimentally. Formation of the central complex is predicted to cause a dramatic reduction in the surface counterion (e.g., Na+) concentration over a region including but extending well beyond the location of the ligand binding site. For binding a cationic ligand, SaKobs is predicted to be negative, indicating net electrolyte ion release in the binding process. At small enough N, -SaKobs is predicted to decrease strongly toward zero with decreasing N. At intermediate N, -SaKobs appears to exceed its limiting value as N-->infinity.

Monte Carlo track structure studies of energy deposition and calculation of initial DSB and RBE

Estimation of exposure due to environmental and other sources of radiations of high-LET and low-LET is of interest in radiobiology and radiation protection for risk assessment. To account for the differences in effectiveness of different types of radiations various parameters have been used. However, the relative inadequacy of the commonly used parameters, including dose, fluence, linear energy transfer, lineal energy, specific energy and quality factor, has been made manifest by the biological importance of the microscopic track structure and primary modes of interaction. Monte Carlo track structure simulations have been used to calculate the frequency of energy deposition by radiations of high- and low-LET in target sizes similar to DNA and higher order genomic structure. Tracks of monoenergetic heavy ions and electrons were constructed by following the molecular interaction-by-interaction histories of the particles down to 10 eV. Subsequently, geometrical models of these assumed biological targets were randomly exposed to the radiation tracks and the frequency of energy depositions obtained were normalized to unit dose in unit density liquid water (l0(3) kg m-3). From these data and a more sophisticated model of the DNA, absolute yields of both single- and double-strand breaks expressed in number of breaks per dalton per Gray were obtained and compared with the measured yields. The relative biological effectiveness (RBE) for energy depositions in cylindrical targets has been calculated using 100 keV electrons as the reference radiation assuming the electron track-ends contribution is similar to that in 250 kV X-ray or Co60 gamma-ray irradiations.

Binding of echinomycin to d(GCGC)2 and d(CCGG)2: distinct stacking interactions dictate the sequence-dependent formation of Hoogsteen base pairs

Molecular dynamics simulations have been used to explore the behavior of the complexes of echinomycin with the DNA tetramers d(GCGC)2 and d(CCGG)2 in which the terminal bases have been paired according to either a Hoogsteen or a Watson-Crick hydrogen bonding scheme. The energy of the four resulting complexes has been monitored along the dynamics trajectories and the interaction energy between echinomycin and DNA has been decomposed into contributions arising from the planar aromatic systems and the depsipeptide part of the antibiotic. Our calculations predict a large increase in overall stabilization upon protonation of the terminal cytosines and subsequent Hoogsteen pair formation in the complex of echinomycin with d(GCGC)2 but not with d(CCGG)2, in agreement with the experimental evidence [Gao and Patel, Quart. Rev. Biophys. 22, 93-138 (1989)]. The conformational preferences appear to arise mainly from differential stacking interactions in which the electrostatic component is shown to play a dominant role. Differences in hydrogen bonding patterns are also found among the complexes and these are compared in relation to available crystal structures. The binding of echinomycin to DNA appears as a complex process involving many interrelated variables.

Dynamics and relative stabilities of parallel- and antiparallel-stranded DNA duplexes

The dynamics and stability of four DNA duplexes are studied by means of molecular dynamics simulations. The four molecules studied are combinations of 4, 15 bases long, single-stranded oligomers, F1, F2, F3, and F4. The sequence of these single strand oligomers are chosen such that F1-F2 and F3-F4 form parallel (ps) DNA double helices, whereas F1-F4 and F2-F3 form antiparallel-stranded (aps) DNA double helices. Simulations were done at low (100 K) and room (300 K) temperatures. At low temperatures the dynamics are quasi-harmonic and the analysis of the trajectories gives good estimates of the low frequency vibrational modes and density of states. These are used to estimate the linear (harmonic) contribution of local fluctuations to the configurational entropy of the systems. Estimates of the differences in enthalpy between ps and aps duplexes show that aps double helices are more stable than the corresponding ps duplexes, in agreement with experiments. At higher temperatures, the distribution of the fluctuations around the average structures are multimodal and estimates of the configurational entropy cannot be obtained. The multi-basin, nonlinear character of the dynamics at 300 K is established using a novel method which extracts large amplitude nonlinear motions from the molecular dynamics trajectories. Our analysis shows that both ps DNA exhibit much larger fluctuations than the two aps DNA. The large fluctuations of ps DNA are explained in terms of correlated transitions in the beta, epsilon, and zeta backbone dihedral angles.

Molecular dynamics simulations of oligonucleotides in solution: visualization of intrinsic curvature

We have undertaken molecular dynamics simulations on the d(CGCAAAAAAGCG).d(CGCTTTTTTGCG) dodecamer in solution. In this study, we focus on aspects of conformation and dynamics, including the possibility of cross-strand hydrogen bonds. We compare our results with those from crystallography as well as infrared, Raman and NMR spectroscopy and cyclization kinetics. Our method of analysis allows us to visualise the curvature of the helix as a function of time during the simulation. We find that the major distortions of the helix axis path occur at the junctions between the (essentially straight) A-tract and the CG- and GC-tracts, although at one junction this is due to hyperflexibility (i.e., regions of high flexibility with no preferred direction of curvature), while at the other junction a static curvature is found (i.e., a preferred, sustained direction of curvature).