| Molecular Dynamics Studies on Free and Bound Targets of the Bovine Papillomavirus Type I E2 Protein: The Protein Binding Effect on DNA and the Recognition Mechanism Biophysical Journal, Volume 89, Issue 4, 1 October 2005, Pages 2542-2551 D. Djuranovic and B. Hartmann Abstract Molecular dynamics simulations of a total duration of 30ns in explicit solvent were carried out on the BPV-1-E2 protein complexed to a high-affinity DNA target containing the two hydrogen-bonded ACCG.CGGT half-sites separated by the noncontacted ACGT sequence. The analysis of the trajectories focuses on the DNA structure and on the dynamics. The data are compared to those issued from recent simulations made on three free targets that recognize E2 with different affinities. E2 does not drastically perturb the mechanic properties of the free DNA: the structural relationships between the BI/BII backbone substates and some helical parameters are preserved in the complex despite a severe slowing down of the phosphate group motions. The structures of both free and bound half-sites are very close to each other although the conformational space explored by these regions is narrowed when they are contacted by the protein. The enhanced plasticity found in the best free target spacers, mainly manifested through the backbone motions, allows a clear overlap between several free and bound global DNA features such as the base displacement. Furthermore, this flexibility is preserved in the complex. Our results support the hypothesis that E2 takes advantage of free predistorted structures that may minimize the DNA deformation cost. In addition, we observe that E2 is far from totally stiffening the DNA, suggesting that the entropic penalty inherent in the complex formation could be limited. Abstract | Full Text | PDF (300 kb) |
| The RNA Binding Protein Pub1 Modulates the Stability of Transcripts Containing Upstream Open Reading Frames Cell, Volume 101, Issue 7, 23 June 2000, Pages 741-751 Maria J Ruiz-Echevarría and Stuart W Peltz Summary The nonsense-mediated mRNA decay (NMD) pathway functions to degrade transcripts containing nonsense codons. Transcripts containing mutations that insert an upstream open reading frame (uORF) in the 5′-UTR are degraded through NMD. However, several naturally occurring uORF-containing transcripts are resistant to NMD. Here we demonstrate that the GCN4 and YAP1 mRNAs, which contain uORFs, harbor a stabilizer element (STE) that prevents rapid NMD by interacting with the RNA binding protein Pub1. Conversely, a uORF-containing mRNA that lacks an STE, such as CPA1, is degraded by the NMD pathway. These results indicate that uORFs can play a pivotal role regulating both translation and turnover and that the Pub1p is a critical factor that modulates the stability of uORF-containing transcripts. Summary | Full Text | PDF (393 kb) |
| Gcn4 Activator Targets Gcn5 Histone Acetyltransferase to Specific Promoters Independently of Transcription Molecular Cell, Volume 6, Issue 6, 1 December 2000, Pages 1309-1320 Min-Hao Kuo, Elmar vom Baur, Kevin Struhl and C.David Allis Summary Histone acetylation correlates well with transcriptional activity, and histone acetyltransferases (HATs) selectively regulate subsets of target genes by mechanisms that remain unclear. Here, we provide in vivo evidence that the yeast transcriptional activator Gcn4 recruits Gcn5 HAT complexes to selective promoters positioned in natural or ectopic locations, thereby creating local domains of histone H3 hyperacetylation and subsequent transcriptional activation. A significant portion of the Gcn4-targeted histone acetylation by Gcn5 is independent of transcriptional activity. These observations provide strong evidence for promoter-selective, targeted histone acetylation by Gcn5 that facilitates transcription in a causal fashion. In addition, Gcn5 also functions in an untargeted manner to acetylate H3 on a genome-wide scale. Summary | Full Text | PDF (643 kb) |
Copyright © 2000 The Biophysical Society. All rights reserved.
Biophysical Journal, Volume 79, Issue 2, 656-669, 1 August 2000
doi:10.1016/S0006-3495(00)76324-7
Biophysical Theory and Modeling
Sylvie Derreumaux and Serge Fermandjian
, 
Address reprint requests to Dr. Serge Fermandjian, Département de Biologie et Pharmacologie Structurales, UMR 8532 CNRS, Institut Gustave Roussy, 39 rue Camille Desmoulins, 94800 Villejuif, France. Tel.: +33-1-42-11-49-85; Fax: +33-1-42-11-52-76.Recognition of DNA by bZIP (basic leucine zipper) dimers is a delicate process, involving subtle structural and dynamical variations that are generally related to dramatic biological effects. Some bZIP dimers recognize the cAMP-responsive element (CRE) in promoters of ubiquitary genes and the 12–O tetradecanoylphorbol-β-acetate (TPA)-responsive element (TRE) in protooncogene promoters. Proteins of the ATF-CREB subfamily bind to CRE and very weakly to TRE, whereas the yeast transcription factor GCN4 (general control protein 4) and the proteins belonging to the FOS-JUN subfamily bind to both CRE and TRE (Ryseck and Bravo, 1991,Hai and Curran, 1991,Hurst, 1994,Kim and Struhl, 1995). Such cross-bindings could be explained by the resemblances between CRE, d(TGACGTCA)2, and TRE, d(TGACTCA)2 (same TGA:TCA half-sites) and among the bZIP proteins (similar DNA binding sequences). Complex formation will further depend on the relative concentrations of the protein components, which can vary with the cell physiology. Some of them will be response productive and others nonproductive or even detrimental to cell life.
The specific recognition of DNA by proteins is highly dependent on the local and sublocal deformations occurring in the helix. According to Paolella et al, bending of DNA is the major structural factor influencing protein-CRE/TRE recognition. Actually, results of biochemical methods indicate a straight structure for TRE and a structure bent slightly (∼10°) toward the major groove for CRE (Paolella et al,Paolella et al,Sitlani and Crothers, 1998,Hockings et al,Sloan and Schepartz, 1998). Bending of CRE shortens both the distance and the difference in orientation between the half-sites TGA:TCA, seemingly reducing the structural difference between CRE and TRE and allowing the recognition of both CRE and TRE by the same bZIP dimers (Paolella et al,Keller et al). The DNA bending accompanying the complex formation varies as a function of the bZIP dimer. It could be responsible for a change in the dimer orientation in the transcription machinery, influencing the promotion of the gene transcription (Paolella et al,Leonard and Kerppola, 1998).
Our previous NMR refined molecular modeling argued in favor of a straight structure for TRE included in a 15-mer DNA and a ∼30° bend toward the minor groove for CRE at the center of a 16mer DNA (Chaoui et al). A structural analysis of the NMR atomic coordinates of Conte et al reported in the Protein Data Bank provided the same bend toward the minor groove.
The results from biochemical methods (Paolella et al,Paolella et al,Sitlani and Crothers, 1998,Hockings et al,Sloan and Schepartz, 1998) and NMR restrained modeling, while they converge for TRE, are thus diverging for CRE. This motivates the present modeling study of CRE, in which no NMR constraints are applied. We believe this study to be more “realistic,” as it resorts to a molecular dynamics (MD) simulation and includes neutralizing Na+ and water explicitly. To our knowledge, it is the first MD study of the CRE DNA, allowing us to investigate its flexibility and its conformational substates. For the sake of comparison, the selected 16-mer DNA d(GAGATGACGTCATCTC)2 was kept the same as the one previously analyzed by NMR (Chaoui et al). A MD trajectory of 10ns was carried out, as a previous study revealed that simulations longer than 5ns were necessary to reach a full convergence of the trajectory for a system of that size (Young et al). Calculations were performed with the AMBER software and the particle mesh Ewald implementation for the treatment of electrostatics. Previous studies have revealed a good agreement with experimental data, particularly regarding the sequence-specific structure and bending of DNA (Cheatham and Kollman, 1996,Cheatham and Kollman, 1997,Young et al,Young and Beveridge, 1998,de Souza and Ornstein, 1998,Sprous et al). The recently modified version of the force field of Cornell et al., parm98, was used (Cheatham et al).
The results indicate that the CRE DNA, although flexible, is on the average slightly curved toward the major groove (7°-8°). The inherent flexibility of CRE, mostly concentrated at the central CpG step, must allow it to easily adjust its curvature in the different DNA-protein complexes. In the cocrystal with GCN4 (Keller et al), the largest helical deformation involves the central CpG, which exhibits a correlated increase in the roll and a decrease in the twist, bringing the TGA:TCA half-sites closer to each other (distance and orientation) and attenuating the structural differences between CRE and TRE. MD results fit our NMR experimental data but disagree with the refined NMR structures (Conte et al,Chaoui et al), highlighting the impact of the electrostatic parameterization on the refinement of the DNA structure and, more particularly, of the axis curvature.
The MD simulation was performed on the 16-mer duplex d(GAGATGACGTCATCTC)2, which contains the CRE consensus sequence d(TGACGTCA)2 of interest at its center. Three factors led to the choice of this oligonucleotide: it has the minimum size needed to achieve stable interactions with the bZIP transcription factors (Talanian et al); the crystal structure of its complex with GCN4 has been determined (Keller et al); and its NMR refined structure has been reported recently (Chaoui et al). We used the following numbering for the DNA sequence: 5′-(G1A2G3A4T5G6A7C8G9T10C11A12T13C14T15C16: G17A18G19A20T21G22A23C24G25T26C27A28T29C30T31C32)-3′.
The MD simulation was carried out with the suite of programs AMBER 4.1 (Pearlman et al) and the modified version of the Cornell et al. force field (parm95), parm98 (Cornell et al,Cheatham et al). This version provides DNA structures that are in better agreement with experimental data. Water was simulated with the TIP3P model (Jorgensen et al). Periodic boundary conditions were simulated with a rectangular unit cell. The initial box was truncated to achieve a minimum distance of ∼11Å beyond all DNA atoms in all directions, resulting in a box size of 49Å×49Å×87Å and 4752 water molecules. Thirty neutralizing Na cations were placed with a simple energy minimization algorithm (Leap module of AMBER 5.0). The starting DNA conformation was the AMBER-generated canonical B-DNA double helix. The energy of the system was minimized for 2250 cycles, using a combination of steepest descent and conjugate gradient algorithms. During the energy minimization, the oligomer conformation and the counterion positions were maintained with harmonic restraints (100kcal/mol Å2), with a progressive diminution of the force constants for a final 500 cycles without restraints. The system was then heated to 303K over 10ps, the velocities were rescaled up to 150K, and then the system was coupled to a heat bath with the weak coupling Berendsen algorithm (Berendsen et al). During heating, harmonic constraints (25kcal/mol Å2) were imposed on the atomic positions of the oligomer and of the sodium counterions. After an constant-mass, constant-volume, and constant-temperature (NVT) equilibration period of 5ps at 303K, the simulation was performed under constant-mass, constant-pressure, and constant-temperature (NPT) conditions (303K, 1atm.), using coupling constants of 0.2ps. The restraints were relaxed over a period of 25ps until a free system was achieved, and a final equilibration was carried out for 5ps before the beginning of the 10-ns production phase. Then a Maxwell distribution of the velocities was assigned to the system, followed by 10ps of free equilibration. This procedure was repeated before the beginning of the production phase. After 50ps and 100ps of production, a Maxwell distribution of the velocities was again assigned. The velocities were redistributed to attenuate the influence on the dynamics of the initial configuration of the system. All bond lengths involving hydrogen atoms were restrained using the SHAKE algorithm (tolerance=0.0005) (Ryckaert et al), and a time step of 2 fs was used. Electrostatic interactions were calculated using the particle mesh Ewald summation method (Darden et al,Essmann et al), with a 10−6 Ewald convergence tolerance and a charge grid spacing of ∼1Å for the interpolation in the reciprocal space. A 9-Å residue-based cutoff was applied to the van der Waals and direct space Ewald interactions. Using a modified version of the Sander module of AMBER 4.1 (Young and Beveridge, 1998), we removed the translations and rotations of the DNA every 100 steps of MD, scaling the resulting solute atom velocities to conserve kinetic energy. The calculations were performed on a Silicon Graphics Origin200 workstation.
CURVES version 5.2 was used to analyze the DNA structures stored every picosecond of MD (Lavery and Sklenar, 1988). CURVES fits a curved or a straight global helix axis that follows the double helix. It allows both global and local geometrical analysis, providing the values of the helical parameters and of the backbone dihedral angles. This work refers to the local parameters output by CURVES. To avoid end effects, only the 14 central base pairs were taken into account for the analysis. For the MD simulation (10ns), each parameter was measured 10,000 times.
MD results were compared to the experimental NMR and x-ray data on unbound and bound CRE, respectively. For the 16-mer duplex d(GAGATGACGTCATCTC)2 of interest, one NMR refined structure (Chaoui et al) and one crystal structure of its complex with the GCN4 transcription factor (2.2Å) (Keller et al) (PDB ID: 2dgc) are available. Another NMR study concerns a different oligomer with CRE at its center, d(CATGTGACGTCACATG)2 (Conte et al) (PDB ID: 1saa). Yet the loose NMR constraints used by these authors relative to the helical parameters (internucleotide distances) and phosphate positions (ϵ dihedral) result in a poor experimental assessment of these DNA morphological parameters. Thus comparisons were made with the NMR experimental structure of Chaoui et al, CRE_NMR; with the crystal structure (Keller et al), CRE_GCN4; and with the NMR data of Conte et al for the sugar phases. Concerning the DNA bending, data of reference were provided by biochemical, NMR, and crystallography experiments. We calculated the local parameters and global bending describing the experimental structures, using the published PDB coordinates and CURVES 5.2. The canonical B-DNA and A-DNA conformations cited in the text are the ones generated by the AMBER software (Arnott et al).
A MD trajectory of 10ns was simulated with the 16-mer duplex d(GAGATGACGTCATCTC)2, using the AMBER software and the recent parm98 force field (Cheatham et al). According to Young et al, simulations longer than 5ns are necessary to reach a full convergence of the trajectory for a system with DNA of the above-mentioned size and containing explicit monovalent cations. The atom coordinates were saved and analyzed every picosecond. The results take into account the trajectory after 500ps from its starting point (i.e., data collected through 9.5ns).
In Figure 1a are plotted the root mean square deviations (rmsds) as a function of time, calculated for the coordinates of heavy atoms in the 14 central residues of the DNA, between the instantaneous structures extracted from the molecular dynamics trajectory (MD), CRE_MD, and the canonical B DNA conformation and A DNA conformation. One can see that after 500ps, the DNA has globally relaxed from the starting conformation and follows its own dynamics.
The DNA conformation calculated with the parm98 force field looks like more the canonical B form than the canonical A form, with rmsd values of nearly 2.5Å and 5.5Å, respectively. In comparison, calculations performed with the parm95 force field (unpublished results) provide a DNA conformation intermediate between the canonical B and A DNAs, with rmsd values of nearly 4.5Å.
In the Appendix are reported the averaged values and the standard deviations over the MD trajectory (0.5–10ns) of the interpair, pair-axis, base-base, and backbone structural parameters of the CRE oligomer. The DNA is globally stabilized in a B form, with averaged twist, rise, inclination, x displacement, and sugar phase values of 33.7°, 3.4Å, 0.0°, −1.0Å, and 147°, respectively, confirming that the new parm98 force field leads to DNA conformations that are more consistent with the experimental data compared to the parm95 force field (Young et al; de Souza and Ornstein, 1998; Cheatham and Kollman, 1997; Cieplak et al).
Sugar puckers. The average sugar phase angle in the CRE_MD structure is 147°. This value fits the B-DNA canonical value (162°) (Saenger, 1984) better than the value provided by the elder parm95 force field does (∼120°-125°) (Cheatham et al). The phase angle as a function of the nucleotide (Fig. 2) follows the well-known sequence dependence, with values higher at purines (∼160°) than at pyrimidines (∼135°) (Poncin et al). The transitions toward the Northern state (C3′-endo, ∼18°) are not numerous, and when they occur, they merely last more than 200ps. Exceptions are the sugars of A12 (in the half-site TCA) and of C14 (outside the consensus sequence), which display transitions of nearly 500ps between 1200ps and 1700ps and of 700ps between 9000ps and 9700ps, respectively. These two nucleotides are the only ones surrounded by pyrimidines.
Phosphate conformations. Numerous coupled variations of the dihedral angles ϵ and ζ associated with BI/BII transitions in the phosphate groups are observed during the dynamics, all over the DNA sequence (BI: ϵ in trans and ζ in gauché; BII: ϵ in gauché and ζ in trans) (Grzeskowiak et al,El antri et al,Gorenstein, 1994). The angle values are ∼180° (280°) for ϵ and ∼270° (150°) for ζ in the BI (BII) state. The difference ϵ−ζ is nearly −90° for BI and 130° for BII and is convenient for following the transitions. The value of ϵ−ζ averaged over time is proportional to the BII population percentage (Gorenstein, 1994). The BII population percentage averaged over the CRE_MD structure (21%) is in very good agreement with the percentage provided by NMR for solution DNAs and by x-ray for crystal DNAs (20–30%) (Gorenstein, 1994,Winger et al). Again the results are significantly improved compared to the results obtained with the parm95 force field (∼10% of BII) (Cheatham et al).
The histogram of the BII population as a function of the nucleotide (Figure 3a) indicates that a period of 10ns is sufficient to render the sequence dependency (similar variations along the two palindromic strands). Nevertheless, a significant difference (up to 20%) is found between the BII population of some symmetrical phosphates (for example, between G6pA7 and G22pA23). A longer MD simulation is needed to sample the phosphate conformations compared with the helical parameters (∼10ns versus ∼2ns). This is due to the potential energy barrier between the BI and the BII states and to the dependence of the phosphate conformations on the diffusion of water molecules (Winger et al).
The largest BII percentages are found within the CRE consensus segment, that is, the half-sites TGA:TCA and the central CpG linker: 70% at G6pA7/G22pA23, 45% at C8pG9/C24pG25 and C11pA12/C27pA28. With exception of G6pA7/G22pA23 and T5pG6/T21pG22, where BII rates are unexpectedly very high and low, respectively, the hierarchy of BII percentages (Figure 3a) is similar to YpR>RpR>YpY>RpY (where Y and R stay for pyrimidine and purine, respectively), obtained from statistical analysis of NMR and crystallographic data (El antri et al,Gorenstein, 1994,Winger et al,Bertrand et al).
The twist angle of CRE_MD averaged over the whole sequence, 33.7°, is that expected for B-DNAs in solution (34°) (Ulyanov and James, 1995). The averaged rise ranges between 3.2Å (A4pT5) and 3.8Å (C8pG9), and the averaged twist varies from 27.2° (A4pT5) to 44.4° (C8pG9) (Appendix). The averaged inclination and x displacement vary slightly, with values ranging from −1.8° and 2.6°, and from −1.2Å and −0.7Å, respectively. The base-base buckle and propeller twist parameters are submitted to important variations, with averaged values between −13.9° (T15:A18) and 14.1° (A2:T31), and between −16.0° (C14:G19) and −3.7° (C11:G22), respectively. No disruption of hydrogen bonds is observed. In the phosphodiester backbone, most variable dihedral angles are ζ (Δ=100°), ϵ (Δ=70°), sugar phase (Δ=50°), and χ (Δ=30°). In contrast, the α, β, and γ angles are roughly constant, staying in g−, t, and g+ conformations, respectively, along the sequence. At last, both the minor and the major groove widths are typical of B-DNAs, with averaged values over CRE of 6.3Å and 12.1Å, respectively. Yet the grooves are narrower at the C8pG9 step in the center of CRE, with width values of nearly 5Å (minor) and 11Å (major).
There is a sampling of crescent-shaped DNA structures with helix curvatures mostly concentrated on the central CpG. The DNA global curvature is defined by CURVES 5.2 as the angle between the local helical axes of the penultimate base pairs A2:T31 and T15:A18 (Lavery and Sklenar, 1988) (Fig. 4a, left). The global curvature (Figure 4b) varies from 0° to 51° (14 bp), with an averaged value of 16.1° and a SD of 8.5°. The bending direction is given by the Atip component of the kink angle (Fig. 4a, right). The kink is defined as the angle between the two optimal linear axes for the DNA half-sides, fragments A2:T31 to C8:G25 and G9:C24 to T15:A18, and the Atip is the component of the kink angle in the plane perpendicular to the long base-pair axis of the central base-pair step. The Atip component is negative or positive for a bending of the oligomer toward the minor or major groove, respectively. The variation of the Atip as a function of time (Figure 4c) shows that the CRE curvature oscillates between the major and the minor grooves, with an intermediate straight conformation, revealing no apparent anisotropy (examples of structures are shown in Fig. 5, left and right). The averaged Atip value reported as a function of the simulation length (Figure 4d) reveals that a trajectory of more than 5ns is necessary to get a convergence of the DNA global bending. The distribution of the Atip values calculated over the entire trajectory (0.5–10ns) is roughly Gaussian (data not shown), indicating that a trajectory of 10ns is sufficient to sample the bending of a DNA 16-mer.
An understanding of why and how bending of CRE occurs requires examination of local perturbations in the double helix. Typical parameters responsible for DNA bending include the roll and, to a lesser extent, the tilt (EMBO Workshop, 1989,Young et al,Suzuki and Yagi, 1995,Dickerson, 1998). The roll is positive or negative for pinching toward the major or the minor groove of the DNA, respectively. The tilt is positive or negative for a pinching toward the second or the first strand backbone of the DNA, respectively. The MD structure shows loci of curvature, mainly in the central consensus sequence. The pyrimidine-purine steps T5pG6, C8pG9, and C11pA12 display averaged positive roll values of 8.9°, 5.7°, and 8.6°, respectively, and high standard deviation values of 6.9°, 7.3°, and 6.8°, respectively (Figure 6a). The major tilt deformations are found at the T5pG6 and G6pA7 steps (2.4°±5.5° and −3.2°±4.5°), at the symmetrical T10pC11 and C11pA12 steps (4.4°±4.7° and −3.4°±5.6°), and at C8pG9 (1.0°±5.7°) (Figure 6b). The CpG and TpG/CpA steps appear as the most malleable steps, which conforms to the crystallographic (Young et al,Suzuki and Yagi, 1995,Dickerson, 1998) and NMR (Ulyanov and James, 1995,El antri et al,Lefebvre et al) data; the effect is related to a lower stacking energy in these steps (Bertrand et al).
The small tilts at T5pG6+G6pA7 and T10pC11+C11pA12 are in phase and additive (Dickerson, 1998) with the roll at C8pG9, resulting in a total kink angle of 7.5° for the CRE consensus segment. Thus the CRE DNA segment is globally curved on the average by 7°-8° toward the major groove at its center, i.e., toward the dimer bZIP protein partner, with a global “C” form. The positive rolls (4.2°) at the A2pG3 and C14pT15 steps ending the DNA duplex are roughly in antiphase with the positive roll at the central C8pG9, as these steps are separated by 6 bp. As a consequence the Atip is reduced, explaining the apparent isotropy of the Atip curve (Figure 4c). In the binding half-sites, the high positive rolls at the T5pG6 and C11pA12 steps create local bends toward the major groove. Globally, the molecule is S-shaped in the plane perpendicular to the Atip one, with each half-site precurved toward one monomer of the bZIP protein.
The plot of rolls averaged over the MD structures and classified as a function of their Atip value (Atip>10°,<−10°, or ≅ 0°) reveals a strong correlation between the roll of the central CpG and the Atip (Fig. 7). The structure of CpG modulates the overall bending: an increase in roll changes a global bending toward the minor groove (Atip<−10°) to a global bending toward the major groove (Atip>10°); a null roll is related to a global bending toward the minor groove because of the presence of positive rolls at the A2pG3 and C14pT15 steps 6 bp from CpG. Regarding the tilts, no clear correlation with the Atip was found.
The consensus CRE segment in CRE_MD (7–8°) presents an averaged curvature toward the major groove very similar to that obtained (nearly 10°) by phasing analysis (Paolella et al), cyclization kinetics (Hockings et al), and ligation ladder (Sloan and Schepartz, 1998). It can be noted, however, that the DNA bending is not due to the sole rolls at the TpG/CpA steps, as suggested by Sloan and Schepartz, 1998, but rather to a favorable combination of roll and tilts at the central CpG and the TGA half-sites, respectively.
The CRE_NMR structure (Chaoui et al) has been calculated with the JUMNA 10.0 molecular mechanics software (Lavery, 1988,Lavery et al), using the JUMNA force field, Flex (Lavery et al,Lavery et al,Lavery et al). A set of different initial structures corresponding to local energy minima were generated by JUMNA. These were minimized with the NMR restraints on the phases, the differences ϵ−ζ, and the 1H-1H distances, introduced step by step (Chaoui et al). All of the initial structures converged toward a final unique structure, CRE_NMR, with a curvature directed toward the minor groove. When the force field Flex was replaced by parm95 (Cornell et al), one obtained a final conformation with a rmsd of 1.0Å to CRE_NMR (Chaoui et al).
In Figure 1b are plotted the rmsds as a function of time between the instantaneous MD structures, CRE_MD, and the available experimental structures. CRE_MD is clearly closer to CRE_GCN4 (rmsds ∼2Å) than to CRE_NMR (rmsds ∼3.5Å). CRE_NMR is found to be very different from 〈CRE_MD〉 and CRE_GCN4 (rmsd=3.6Å).
Backbone conformation. The sugar phase angles in CRE_MD were compared to the NMR data of Chaoui et al. There is a good agreement between calculation and experiments, with the same phase angle variation as a function of the nucleotide sequence. The phase angles averaged over the MD trajectory stay within the NMR uncertainties, except that of the C14 sugar (Fig. 2). The agreement between the current MD data and the NMR data of Conte et al related to the CRE segment is also good. Sugars at C14/C30 show a preference for the Eastern state (O4′-endo, ∼90°) in CRE_NMR and for the C1′-exo state (∼126°) in CRE_MD. Actually, during the MD trajectory, these two sugars keep the C2′-endo conformation most of the time, with phase angles of nearly 150°. As they transit from time to time into the Northern conformation, their averaged conformation lies in the range of C1′-exo puckers. The Eastern state provided by NMR could better correspond to an average of Southern and Northern states. Whatever the real conformation adopted by residues C14/C30, there is a disagreement between the MD calculations and the NMR measurements.
Regarding the phosphate conformations, the BII populations predicted by MD simulation are close to those estimated from NMR 31P chemical shifts, using the Gorenstein, 1994 relationship (Figure 3b). For the average over the sequence, both NMR and MD provide 21% of BII population, whereas for the BII population at each step, the absolute difference does not exceed 10%. Exceptions are G6pA7/G22pA23 and their complementary steps T10pC11/T26pC27. For the first steps, relative BII populations of 70% and only 25% are predicted by MD and NMR, respectively, while for the second steps MD predicts a smaller BII population (−20%). Globally, the sequence dependence of the BII population is found to be more pronounced by MD calculations than by NMR measurements (1 SD=21% versus 12%).
Helical parameters and bending propensity. The current MD simulation disagrees with the recently reported NMR restrained molecular modelings of both Conte et al and Chaoui et al. These show a curvature of 30° (16-mer) toward the minor groove and rolls negative at CpG (−19°/−10°) and null or negative at TpG:CpA (0°/−12°). Actually, the plot of rolls as a function of dinucleotide steps for CRE_NMR (Chaoui et al) is the opposite of that of CRE_MD (Figure 6a), but this also proves to be true for the tilts at G6pA7 and T10pC11. It seems that tilts adapt themselves to smooth the bending caused by rolls. In addition, the averaged twist in CRE_MD is smaller compared with that found in CRE_NMR (33.7° versus 35.4°). Yet the profiles of the twist sequence dependences are very similar (Figure 6c), with the greatest differences being observed at CpG and T5pG6 and C11pA12. The values of CRE_NMR are within the averaged values±1 SD of CRE_MD.
Influence of software and protocol on NMR structure determination. The curvature of CRE predicted by MD is similar to that derived from biochemical data but differs from that provided by NMR refinement methods.
DNA conformations resolved by NMR refinement depend on different factors such as measurements and methodology. The latter includes the modeling software and the refinement protocol. The helical parameters are determined by the internucleotide distances obtained from NOE, data which are refined through back-calculations of theoretical NOESY spectra, using the NOE Simulate routine of INSIGHT II (MSI) (Boelens et al,Lefebvre et al,Lefebvre et al,Chaoui et al). Moreover, the bending can be strongly influenced by the force field. Its impact can be tested by comparing, for instance, the rolls in the minimal energy structures provided by JUMNA, using Flex and parm95. Conditions are the same as the ones used for the CRE refinement by Chaoui et al (dielectric function and screening of the phosphate charges), but the NMR constraints are excluded. One finds that rolls are modulated by the force field (Figure 8a), especially at the CpG and TpG/CpA steps, where JUMNA-Flex provides negative rolls and JUMNA-parm95 provides positive ones. Replacing parm95 with parm98 (force-field parameters modified for three dihedral angles, no change of the partial charges) in JUMNA does not significantly alter the rolls: they differ by less than 2° at each step of the CRE sequence. As in most softwares dedicated to the NMR refinement of molecular structures, the DNA environment in JUMNA is simulated in an implicit way: the dielectric effect of water is modeled by a sigmoid function with a slope value of 0.16 and a plateau of 78, and the screening effects of counterions on phosphates are taken into account by reducing the charge of the phosphate groups to a value of −0.5e (Lavery et al). Using different slope values of the dielectric function (0.16 and 0.356) and varying the apparent charge of the phosphates (−0.5e, 0.0e, and −1.0e) for the minima provided by JUMNA-parm95 shows that the roll values increase parallel to the dielectric function and when the phosphate charge tends to zero. A MD simulation with the same parm95 force field (unpublished results) with an explicit consideration of the solvent provides a similar curve, but with more positive values of rolls (Fig. 8, b and c). However, we cannot conclude that it is possible to find an implicit electrostatic parameterization that would correctly simulate, whatever the DNA sequence, the influence of the solvent on the DNA structure. Actually, structure alterations caused by specific localization of water molecules or hydrated ions (Diekmann, 1987,Rouzina and Bloomfield, 1998,Lyubartsev and Laaksonen, 1998,Winger et al,McFail-Isom et al,Shui et al) cannot be simulated by varying the apparent phosphate charge or by using a dielectric function.
From the above results one can conclude that the electrostatic parameterization of the software exerts a significant effect on DNA bending, through both the force field (higher polarity of the atomic charges for parm95/parm98 than for Flex) and the solvent model. As far as CRE is concerned, a major factor could be the electronegativity of the major groove along the consensus d(TGACGTCA)2 helix. Actually, CRE displays a succession of intrastrand electronegative pockets, at 5′-T5pG6-3′, 3′-T26pG25-5′, 5′-G9pT10-3′, and 3′-G22pT21-5′, connected one to another by interstrand electronegative pockets 5′-G6pA7-3′, 5′-C8pG9-3′, and 5′-T10pC11-3′ (Young et al). The succession of thymine O4 and guanine O6 and N7 atoms creates a continuous negative channel along the major groove, which could push away subsequent base pairs and induce negative (or less positive) rolls. The repelling of base pairs would in this case decrease as the dielectric effect increases.
The NMR restraints introduced successively on the sugar phases and on the ϵ−ζ differences did not significantly modify the roll values. In contrast, the restraints on 1H-1H distances, particularly the H6/8-H6/8 distances, have both “pulled” the rolls of CpG and TpG/CpA toward strong negative values and bent CRE toward the minor groove. This occurs whatever the force field used, Flex or parm95 (Fig. 9) (Chaoui et al). Although the interbase H6/8-H6/8 distance between subsequent bases is related to the roll value, it is not sufficient to determine precisely the roll because, for the same H6/8-H6/8 distance, different positive or negative rolls are possible, depending on the twist value. Actually, both the H6/8-H6/8 distances and the H3′-H6/8 distances (which are not always amenable to 2D NMR experiments for signal overlapping reasons) are needed to determine the rolls univocally (Lefebvre et al). For instance, in the CRE NMR restrained structures obtained with JUMNA-Flex and JUMNA-parm95 by Chaoui et al, the roll values of the central CpG are different, although these structures present similar H6/8-H6/8 distances (Fig. 10). On the other hand, the roll values in the CRE structures determined by Conte et al and Chaoui et al (with JUMNA-parm95) are similar, although the H6/8-H6/8 distances are very different. In these structures the restraints on the H3′-H6/8 distances strongly related to the rolls are severely missed (Lefebvre et al). Perhaps, whatever the starting structure, the software reacts to the adjunction of the H6/8-H6/8 distance restraints by “pushing” the CRE conformation toward a region of the conformational space that is “more convenient” for it (i.e., abnormally negative rolls but correct twists).
Another obvious fault of the refinement procedure that could influence the rolls concerns the internal molecular motions, so that the refinement procedure provides erroneous internucleotide 1H-1H distances, especially at the most flexible steps (TpG:CpA and CpG).
The disagreement on the bending of CRE_NMR and CRE_MD could also result from the solvent conditions, which are different in the NMR and MD experiments, with regard to both the ionic strength (MD ≅ 0.5; NMR ≅ 0.1) and the nature of the ions (MD: Na+; NMR: buffer+purification+impurities → Na+, Li+, K+, and perhaps also traces of divalent cations like Mn2+, Mg2+, Zn2+). Different monovalent ions could induce different local DNA deformations, as crystallographic (Tereshko et al) and MD (Lyubartsev and Laaksonen, 1998) studies performed with various monovalent cations argue that they bind to different specific DNA sites (for instance, Na+ preferentially binds to the minor groove and Li+ binds to the phosphate oxygens). In a recent MD study, Young and Beveridge, 1998 have shown that the nature and the concentration of monovalent cations can have a significant impact on the DNA curvature: replacing neutralizing Na+ with K+ and KCl doubles the bending of a 25-bp DNA segment containing an A-tract. Moreover, gel retardation (Diekmann, 1987), theoretical investigations (Rouzina and Bloomfield, 1998), and crystallographic data (McFail-Isom et al) suggest that divalent cations undergo sequence-specific interactions with DNA that have a particular impact on the DNA bending. Thus the nonequivalent ionic conditions between the NMR and MD studies could in principle contribute to the appearance of diverging results regarding the bending of CRE. As the MD results agree with the NMR data regarding the conformation of the phosphodiester backbone (sugar phases and phosphate conformations), this would mean that the diverging ionic conditions would only influence the helical parameters, including the rolls. But it is interesting to note that the biochemical data on the bending of CRE support the MD results, although the experimental conditions diverge from the MD conditions. For instance, the measurements of cyclization kinetics (Hockings et al) have been made in the presence of KCl (50mM) and MgCl2 (5mM), and the ligation ladder experiments could have been performed in the presence of Mg2+ impurities (Sloan and Schepartz, 1998).
In the end, the influence of the phosphorus NMR restraints has also to be considered. Such restraints could force the steps involved in a BI/BII equilibrium to adopt time-averaged ϵ−ζ values devoid of physical sense. A forced backbone may strongly alter the shape of the helix because of the relation that exists between ϵ−ζ and the helical parameters, including the rolls that are strongly involved in the DNA curvature (El antri et al,Hartmann et al,Gorenstein, 1994,Lefebvre et al,Bertrand et al). Actually, one could generate different initial structures with each phosphate either in the BI or in the BII conformation, making a combinatorial search with all of the steps susceptible to adopting both BI and BII states. One could then perform an energy minimization without constraint on the ϵ and ζ angles on the different initial structures and obtain final structures with different patterns of BI/BII phosphates. Finally, one could compare the averaged ϵ−ζ values to the NMR data (Lefebvre et al,Tisné et al).
The rmsd between 〈CRE_MD〉 (averaged MD structure without further minimization) and CRE_GCN4 is only 1.3Å (Table 1 and Fig. 5, center). This is lower compared to the averaged rmsd of the instantaneous MD structures (∼2Å) (Figure 1b), probably because, at a first approximation, it does not take into account the vibrational motions of the DNA molecule. The good agreement between the simulated CRE structure and the crystal structure of CRE bound to GCN4 indicates that the binding of the protein does not significantly affect the DNA, in agreement with the biochemical results (Paolella et al,Paolella et al,Hockings et al,Sloan and Schepartz, 1998). The rmsd for the central octamer is higher compared with that for the entire 14-mer (1.4Å versus 1.3Å), translating a larger DNA deformation at the binding sequence. Actually, in CRE_GCN4 the octamer helix has drifted away from the canonical B form toward the canonical A form: the rmsds from the B DNA (A DNA) are 1.5Å (4.8Å) and 2.1Å (4.3Å) for 〈CRE_MD〉 and CRE_GCN4, respectively. Four sugars have smaller phase angles in CRE_GCN4 compared with the free oligomer (Fig. 2). For the other sugars, the low resolution of the crystal structure (2.2Å) makes the precise determination of phase angles difficult, preventing a comparison with the MD results. In the same way, in CRE_GCN4 the protein blocks all of the phosphates, except at the TpG:CpA steps, in the BI conformation (ϵ−ζ ≈ −90°) (Figure 3c), a characteristic of the A DNA helix. Such a preference of phosphate groups for the BI conformation has already been observed in DNA-protein complexes (Gorenstein, 1994). Note that previous mechanical studies have shown that the phosphate conformations correlate with the twists: the twist angle increases while the phosphate tends toward BII (Hartmann et al,Bertrand et al). This correlation is roughly observed in CRE_MD, but not at all in the crystal structure (Figure 3cc and Figure 6cc).
| Table 1 Rmsds (in Å) between the averaged MD, experimental, and canonical structures |
| 〈CRE_MD〉 | B DNA | A DNA | CRE_NMR | CRE_GCN4 | |||
|---|---|---|---|---|---|---|---|
| 〈CRE_MD〉 | 1.5 | 4.8 | 2.2 | 1.4 | |||
| B-DNA | 2.2 | 5.4 | 1.6 | 2.1 | |||
| A-DNA | 6.1 | 7.5 | 6.0 | 4.3 | |||
| CRE_NMR | 3.1 | 2.1 | 8.4 | 2.7 | |||
| CRE_GCN4 | 1.3 | 2.7 | 5.8 | 3.6 | |||
| Root mean square deviations (rmsd) between the structure of (AGATGACGTCATCT)2 averaged over the MD trajectory (9.5ns), 〈CRE_MD〉; the canonical B DNA and A DNA conformations; the NMR refined structure by Chaoui et al, CRE_NMR; and the conformation of the DNA in the cocrystal with GCN4 (Keller et al), CRE_GCN4. The upper right part of the table gives the rmsd values calculated for the CRE consensus sequence, i.e., the central octamer (TGACGTCA)2; and the lower left part of the table gives the rmsd values for the whole 14-mer. |
The crescent structure of CRE_MD resembles that displayed by CRE_GCN4, with a global curvature toward the major groove of 16.1° in the former and 15° in the latter. Rolls at the CpG and the TpG/CpA steps are positive, with, however, a higher value for CpG (12.3° versus 5.7°) and a slightly smaller value for TpG/CpA (7.3° versus 8.9° and 8.6°) (Figure 6a). A significant difference consists of rolls at G6pA7 and T10pC11 that are positive in the crystal and almost null in CRE_MD (7.1° versus 1.8° and 0.8°). However, the roll values of CRE_GCN4 stay within the averaged MD values±1 SD, suggesting that the small DNA bending deformation induced by GCN4 is accessible through thermal energy. The straight conformation adopted by CRE in the complex with BP1 (Paolella et al) could also be reached with little consumption of energy, as a null roll at CpG stays within the MD averaged value±1 SD. Rolls with averaged values and SD higher at the CpG and TpG/CpA steps reflect the malleability of the helix at both the linker and the half-site positions. The roll differences between CRE_MD and CRE_GCN4 are caused by a drift of the kinks from TpG/CpA to CpG induced by the protein, thus confirming the gel results (Sloan and Schepartz, 1998).
Differences in tilts are observed only for G3pA4 and the T13pC14 steps (Figure 6b). The twist angle of CRE_MD averaged over the whole sequence (33.7°) is similar to that found in CRE_GCN4 (33.8°). However, although similar for the residues lying outside the CRE consensus octamer, the twists are very different inside CRE, especially at the CpG step (29.9° in CRE_GCN4 versus 44.4° in CRE_MD). Such a twist variation at CpG (correlated to the roll) is not accessible through thermal energy (Figure 6c) and must be induced by work done during complexation.
The comparison of CRE with TRE seemed interesting, as both are stable partners of GCN4. As far as the distance between the half-sites TGA/TCA and their relative orientations are concerned, the structural difference between CRE and TRE is reduced upon protein binding, as TRE keeps a straight conformation (Ellenberger et al,Keller et al). The fact that a greater positive roll is mechanically correlated to a smaller twist (Hartmann et al,Bertrand et al) suggests that the coupled change in the central CpG roll and twist during the formation of the complex of CRE with GCN4 follows a general mechanism in DNA that allows a reduction of the structural difference between promoters with half-sites that diverge by the spacer length.
A long molecular dynamics simulation of 10ns was made on a 16mer DNA sequence containing the CRE consensus sequence. Results were compared with those provided by NMR and biochemical methods. The B-DNA conformation presented by CRE oscillates between a helix curved toward the major groove and a helix curved toward the minor groove, through an intermediate straight form. However, as a final result, the averaged helix is bent slightly toward the major groove. Although this finding is consistent with the biochemical results, it disagrees with the refined NMR structures, where the bending is directed toward the minor groove. This reversion is imputed to an improper electrostatic parameterization in the NMR structure refinement softwares, inasmuch as the NMR data agree with the MD simulation results. In contrast, the MD structure appears to be very similar to the crystal structure of CRE bound to GCN4. In particular, the rolls (often involved in the DNA curvature) are roughly conserved in the complex. These parameters could be the structural determinants of the CRE recognition by GCN4. A large part of the helix deformation induced by GCN4 is concentrated at the central CpG step that links the TGA:TCA half-sites. CpG undergoes a concerted increase in roll and a decrease in twist that reduces both the distance between the TGA:TCA half-sites and their difference in orientation. As a result, the structure of CRE bound to GCN4 is more like the structure of TRE, which remains straight in its complex with GCN4.
Thus the current MD simulation reveals the high flexibility of the CRE consensus sequence, reflecting its ability to accommodate itself to various bZIP proteins. The impact of the central CpG dinucleotide on the conformation of CRE is clearly identified: the inherent bending and malleability of CpG are the major contributions to the curving and flexibility of CRE. The absence of CpG at the center of TRE could be the structural explanation for the observation by biochemical methods that TRE is unable to bind to such varied bZIP partners and in such numbers as CRE does.
Arnott et al., 1972 (1972). Optimized parameters for A-DNA and B-DNA. Biochem. Biophys. Res. Commun. 47, 1504–1509. CrossRef | PubMed
Berendsen et al., 1984 (1984). Molecular dynamics with coupling to an external bath. J. Chem. Phys. 81, 3684–3690. CrossRef | PubMed
Bertrand et al., 1998 (1998). Flexibility of the B-DNA backbone: effects of local and neighbouring sequences on pyrimidine-purine steps. Nucleic Acids Res. 26, 1261–1267. CrossRef | PubMed
Boelens et al., 1989 (1989). Iterative procedure for structure determination from proton-proton NOEs using a full relaxation matrix approach. Application to a DNA octamer. J. Magn. Reson. 82, 290–308. PubMed
Chaoui et al., 1999 (1999). An intrinsic curvature towards the minor groove in the cAMP-responsive element DNA found by combined NMR and molecular modelling studies. Eur. J. Biochem. 259, 877–886. CrossRef | PubMed
Cheatham et al., 1999 (1999). A modified version of the Cornell et al. force field with improved sugar pucker phases and helical repeat. J. Biomol. Struct. Dyn. 16, 845–862. PubMed
Cheatham and Kollman, 1996 (1996). Observation of the A-DNA to B-DNA transition during unrestrained molecular dynamics in aqueous solution. J. Mol. Biol. 259, 434–444. CrossRef | PubMed
Cheatham and Kollman, 1997 (1997). Molecular dynamics simulations highlight the structural differences among DNA:DNA, RNA:RNA, and DNA:RNA hybrid duplexes. J. Am. Chem. Soc. 119, 4805–4825. CrossRef | PubMed
Cieplak et al., 1997 (1997). Molecular dynamics simulations find that 3′ phosphoramidate modified DNA duplexes undergo a B to A transition and normal DNA duplexes an A to B transition. J. Am. Chem. Soc. 119, 6722–6730. CrossRef | PubMed
Conte et al., 1997 (1997). Solution structure of the ATF-2 recognition site and its interaction with the ATF-2 peptide. Nucleic Acids Res. 25, 3808–3815. CrossRef | PubMed
Cornell et al., 1995 (1995). A second generation force field for the simulation of proteins, nucleic acids, and organic molecules. J. Am. Chem. Soc. 117, 5179–5197. CrossRef | PubMed
Darden et al., 1983 (1983). Particle mesh Ewald: an N·log(N) method for Ewald sums in large systems. J. Chem. Phys. 98, 10089–10092. CrossRef | PubMed
de Souza and Ornstein, 1998 (1998). Inherent DNA curvature and flexibility correlate with TATA box functionality. Biopolymers 46, 403–415. CrossRef | PubMed
Dickerson, 1998 (1998). DNA bending: the prevalence of kinkiness and the virtues of normality. Nucleic Acids Res. 26, 1906–1926. CrossRef | PubMed
Diekmann, 1987 (1987). Temperature and salt dependence of the gel migration anomaly of curved DNA fragments. Nucleic Acids Res. 15, 247–265. PubMed
El antri et al., 1993 (1993). Effect of distortions in the phosphate backbone conformation of six related octanucleotide duplexes on CD and 31P NMR spectra. Biochemistry 32, 7079–7088. PubMed
Ellenberger et al., 1992 (1992). The GCN4 basic region leucine zipper binds DNA as a dimer of uninterrupted α helices: crystal structure of the protein-DNA complex. Cell 71, 1223–1237. Abstract | | CrossRef | PubMed
EMBO Workshop, 1989 (1989). EMBO workshop. EMBO J. 8, 1–4. PubMed
Essmann et al., 1995 (1995). A smooth particle mesh Ewald method. J. Chem. Phys. 103, 8577–8593. CrossRef | PubMed
Gorenstein, 1994 (1994). Conformation and dynamics of DNA and protein-DNA complexes by 31P NMR. Chem. Rev. 94, 1315–1338. CrossRef | PubMed
Grzeskowiak et al., 1991 (1991). The structure of B-helical C-G-A-T-C-G-A-T-C-G and comparison with C-C-A-A-C-G-T-T-G-G. J. Biol. Chem. 266, 8861–8883. PubMed
Hai and Curran, 1991 (1991). Cross-family dimerization of transcription factors Fos/Jun and ATF/CREB alters DNA binding specificity. Proc. Natl. Acad. Sci. USA 88, 3720–3724. CrossRef | PubMed
Hartmann et al., 1993 (1993). BI-BII transitions in B-DNA. Nucleic Acids Res. 21, 561–568. PubMed
Hockings et al., 1998 (1998). Characterization of the ATF/CREB site and its complex with GCN4. Proc. Natl. Acad. Sci. USA 95, 1410–1415. CrossRef | PubMed
Hurst, 1994 (1994). Transcription factors 1: bZIP proteins. Protein Profile 1, 123–168. PubMed
Jorgensen et al., 1983 (1983). Comparison of simple potential functions for simulating liquid water. J. Chem. Phys. 79, 926–935. CrossRef | PubMed
Keller et al., 1995 (1995). Crystal structure of a bZIP/DNA complex at 2.2Å: determinants of DNA specific recognition. J. Mol. Biol. 254, 657–667. CrossRef | PubMed
Kim and Struhl, 1995 (1995). Determinants of half-site spacing preferences that distinguish AP-1 and ATF/CREB bZIP domains. Nucleic Acids Res. 23, 2531–2537. PubMed
Lavery, 1988 (1988). Junctions and bends in nucleic acids: a new theoretical modeling approach. In Olson, W.K., Sarma, M.H., Sundaralingam, M., eds. Structure and Expression, DNA Bending and Curvature 3, (New York: Adenine Press), pp. 191–211. PubMed
Lavery et al., 1986a (1986). A general approach to the optimization of the conformation of ring molecules with an application to valinomycin. J. Biomol. Struct. Dyn. 4, 443–461. PubMed
Lavery and Sklenar, 1988 (1988). The definition of generalized helicoidal parameters and of axis curvature for irregular nucleic acids. J. Biomol. Struct. Dyn. 6, 63–91. PubMed
Lavery et al., 1986b (1986). The flexibility of the nucleic acids. II. The calculation of internal energy and applications to mononucleotide repeat DNA. J. Biomol. Struct. Dyn. 3, 989–1014. PubMed
Lavery et al., 1984 (1984). Optimized monopole expansions for the representation of the electrostatic properties of the nucleic acids. J. Comp. Chem. 5, 363–373. PubMed
Lavery et al., 1995 (1995). JUMNA (junction minimisation of nucleic acids). Comp. Phys. Comm. 91, 135–158. PubMed
Lefebvre et al., 1997 (1997). Sensitivity of NMR internucleotide distances to B-DNA conformation: underlying mechanics. Nucleic Acids Res. 25, 3855–3862. CrossRef | PubMed
Lefebvre et al., 1995 (1995). The structural behaviour of the CpG step in two related oligonucleotides reflects its malleability in solution. Biochemistry 34, 12019–12028. PubMed
Lefebvre et al., 1996 (1996). Solution structure of the CpG containing d(CTTCGAAG)2 oligonucleotide: NMR data and energy calculations are compatible with a BI/BII equilibrium at CpG. Biochemistry 35, 12560–12569. PubMed
Leonard and Kerppola, 1998 (1998). DNA bending determines Fos-Jun heterodimer orientation. Nature Struct. Biol. 5, 877–881. CrossRef | PubMed
Lyubartsev and Laaksonen, 1998 (1998). Molecular dynamics simulations of DNA in solution with different counter-ions. J. Biomol. Struct. Dyn. 16, 579–592. PubMed
McFail-Isom et al., 1998 (1998). Divalent cations stabilize unstacked conformations of DNA and RNA by interacting with base π systems. Biochemistry 37, 17105–17111.