| Stability of the I-motif Structure Is Related to the Interactions between Phosphodiester Backbones Biophysical Journal, Volume 84, Issue 6, 1 June 2003, Pages 3838-3847 Thérèse E. Malliavin, Jocelyne Gau, Karim Snoussi and Jean-Louis Leroy Abstract The i-motif DNA tetrameric structure is formed of two parallel duplexes intercalated in a head-to-tail orientation, and held together by hemiprotonated cytosine pairs. The four phosphodiester backbones forming the structure define two narrow and wide grooves. The short interphosphate distances across the narrow groove induce a strong repulsion which should destabilize the tetramer. To investigate this point, molecular dynamics simulations were run on the [d(C2)]4 and [d(C4)]4 tetramers in 3′E and 5′E topologies, for which the interaction of the phosphodiester backbones through the narrow groove is different. The analysis of the simulations, using the Molecular Mechanics Generalized Born Solvation Area and Molecular Mechanics Poisson-Boltzmann Solvation Area approaches, shows that it is the van der Waals energy contribution which displays the largest relative difference between the two topologies. The comparison of the solvent-accessible area of each topology reveals that the sugar-sugar interactions account for the greater stability of the 3′E topology. This stresses the importance of the sugar-sugar contacts across the narrow groove which, enforcing the optimal backbone twisting, are essential to the base stacking and the i-motif stability. Tighter interactions between the sugars are observed in the case of N-type sugar puckers. Abstract | Full Text | PDF (275 kb) |
| Solution Structure of a Substrate for the Archaeal Pre-tRNA Splicing Endonucleases: The Bulge-Helix-Bulge Motif Molecular Cell, Volume 1, Issue 6, 1 May 1998, Pages 883-894 John L Diener and Peter B Moore Summary The structure of the bulge-helix-bulge motif that constitutes the intron/exon splice site in pre-tRNA has been determined by NMR spectroscopy. The conformations of the two 3 nt bulges, where the pre-tRNA is cleaved, are stabilized by stacking interactions between bulge nucleotides and bases in the adjacent Watson-Crick helices and by a network of backbone hydrogen bonds. Both bulges are presented on the same minor groove face of the central 4 bp helix, and the overall structure has approximate two-fold symmetry, which makes it well-suited for attack by archaeal splicing endonucleases, which are symmetric dimers. Summary | Full Text | PDF (2842 kb) |
| The loop E–loop D region of Escherichia coli 5S rRNA: the solution structure reveals an unusual loop that may be important for binding ribosomal proteins Structure, Volume 5, Issue 12, 15 December 1997, Pages 1639-1653 Anne Dallas and Peter B Moore Summary This structure rationalizes all the biochemical and chemical protection data available for the loop E–loop D arm of intact 5S rRNA. While the molecule is double helical over its entire length, the geometry of its internal loop is highly irregular, and its irregularities may explain why the loop E–loop D arm of 5S rRNA interacts specifically with ribosomal protein L25 in . Summary | Full Text | PDF (1297 kb) |
Copyright © 2007 The Biophysical Society. All rights reserved.
Biophysical Journal, Volume 92, Issue 8, 2647-2665, 15 April 2007
doi:10.1529/biophysj.106.092601
Biophysical Theory and Modeling
Thomas Créty and Thérèse E. Malliavin
, 
Address reprint requests to Thérèse E. Malliavin, Laboratoire de Biochimie Théorique, CNRS UPR 9080, Institut de Biologie Physico-Chimique, 13 rue P. et M. Curie, 75 005 Paris, France. Tel.: 33-1-58-41-51-68; Fax: 33-1-58-41-50-26.During the formation of the ribosome 30S subunit, the protein S15 is among the first proteins, along with S18, S6, and S11, interacting with the 30S central domain 1. S15 binds to a phylogenetically conserved three-way junction formed by the intersection of helices H20, H21, and H22 of eubacterial 16S RNA (Fig. 1). The free 16S RNA undergoes an equilibrium 2 between an unfolded conformation, in which the junction is planar with ∼120° interhelical angles, and a folded conformation, in which two helices, H21 and H22, become collinear and the third, helix H20, forms a 60° angle with respect to H22. From the observations made in the literature, it was concluded 2 that the Mg2+ ions as well as the protein stabilize the RNA folded conformation.
The ribosomal protein S15 from Thermus thermophillus was cocrystallized with a 57-nucleotide 16S rRNA fragment containing shortened helices H20, H21, and H22 3. The two helices H21 and H22 are capped by UUCG tetraloops, which are called TL1 and TL2 in the following 4. The 2.8Å resolution of the crystal structure of the complex (Protein Data Bank (PDB) id: 1dk1, Figure 2a) shows 3 that S15 interacts in the minor groove with the three-way junction (first binding site) and a G-U/G-C motif (second binding site). The protein structure consists of four α-helices, the helices α1 (residues 3–15), α2 (residues 24–45), α3 (residues 49–73), and α4 (residues 73–81). The S15 structure in the complex is similar to the structure of the free protein previously determined by NMR (PDB id: 1ab3, Figure 2b) 5. Nevertheless, in the crystallographic structure of the free S15 protein from Bacillus stearothermophilus (PDB id: 1a32, Figure 2c) 6, the N-terminal helix α1 protrudes from the body of the molecule to make contacts with a neighboring protein in the crystal lattice.
Nine Mg2+ ions were positioned in the S15 RNA complex structure according to the observation of strong peaks of electronic density 3, and three of them are supposed to contribute to the junction stability. The Mg2+ ions were previously shown in the literature to play an important role in the mechanism of ribozyme catalysis 7,8, in the folding of RNA 9,10,11,12,13, in the stability of the ribosome 14, and in the interaction between RNA and proteins 15,16 and ligands 17. Molecular dynamics (MD) simulations were performed to study the influence of Mg2+ ions on the conformation of an RNA kissing loop 18, on the conformation of the 5S rRNA loop E motif 19, and on the conformation of the RNA-protein complex between E. coli loop E/helix IV rRNA and L25 protein 20.
The Mg2+ ions are supposed to stabilize the folded form of the RNA three-way junction and thus to enhance the S15 binding 2 by increasing the bimolecular association rate. However, the effect of the Mg2+ ions on the protein-RNA interaction at the atomic level has not yet been investigated. We analyze here MD simulations of the S15-RNA complex in the absence and in the presence of Mg2+ ions to study the influence of Mg2+ on the protein-RNA interaction. MD trajectories of the protein S15 are also analyzed to explore the conformational landscape of the free protein and relate it to the protein behavior in the complex. As said previously, the interaction of S15 with 16S RNA is one of the first steps in the 30S subunit formation, and a detailed inspection of the S15-RNA interaction features may thus be important for targeting the ribosome formation 21. Indeed, S15 was detected in a crystal structure of the 70S ribosome 22 in a bridge spanning the ribosomal subunit interface and seems to be essential for an efficient formation of the ribosome from the 30S and 50S subunits 23.
The use of MD simulations has proved to be a valuable tool for analyzing the behavior of systems including RNA (see McDowell et al. 24 and references cited therein). Nevertheless, the electrostatics of the Mg2+ ions is modeled here using a single point charge without including the polarizability. This approximation is certainly crude but cannot be avoided because introduction of polarizability 25 would result in an excessive computational load for a several-nanosecond trajectory. Another limitation of the MD simulations of nucleic acids is the difficulty of modeling the conformation of the phosphodiester backbone 26. But, there is not much drift of the phosphodiester backbone angles from the crystallographic structure during the time interval studied in the work presented here.
The starting point of the MD trajectories recorded on the protein S15 was the first conformer of the PDB NMR structure 5 of S15 (PDB id: 1ab3) and the smallest-energy conformer of the NMR structure recalculated in the work presented here using the CNS-ARIA protocol. The starting conformation for the simulations of the S15-RNA complex was the x-ray crystallographic structure 3 of the S15-RNA complex (PDB id: 1dk1). Two molecular conformations are present for residues 27–29 of RNA in the 1dk1 file: the conformation A corresponding to the maximum number of formed basepairs was selected.
The parm98 parameter set 27 and the TIP3P model for water 28 were used in simulations. The motion equations were integrated with versions 6 and 7 of the program AMBER 29, respectively, for the isolated protein and for the complex. The cutoff of the Lennard-Jones interaction was 9Å, and the long-range electrostatic interactions were calculated using the Particle Mesh Ewald protocol 30. The Lennard-Jones parameters of the Mg2+ ions were taken from Aqvist 31 and correspond to an integration radius of 0.7926Å and a well depth of 0.8947 kcal/mol. Trajectories of 15ns were recorded for each studied system: the 1ab3 NMR conformer of S15, the 2fkx conformer of S15 calculated in the work presented here, and the x-ray structure of the S15-RNA complex in the presence and the absence of Mg2+ ions, with a temperature of 298K and a pressure of 1atm. The atom coordinates were saved each 1ps.
In case of the simulation of S15-RNA in the presence of Mg2+ ions, the solute includes the nine crystallographic Mg2+ ions but not the three Na+ and K+ ions observed in the crystal. The counterions were added in the vicinity of the charged groups of the solute using a Coulombic potential on a grid, using the command AddIons of leap. The charge of the isolated protein was neutralized with chloride ions, whereas sodium ions were added to neutralize the complex S15-RNA. For the simulation in absence of Mg2+ ions, the Na+ ions were not placed at the positions of removed Mg2+. The resulting modeling is thus related to the presence or absence of a counterion, more than the electronic properties (polarizalibility) of the ions. This counterbalances the limitations of the treatment of electrostatics in the system and the uncertainty of the determination of the ion types in crystallographic structures.
The concentration of Mg2+ ions is ∼40mM, which is a value of the same order as the experimental concentration used to study the S15-RNA interaction 2. On the other hand, the concentrations of Na+ ions in the S15-RNA simulations are ∼140mM in the presence of Mg2+ ions and ∼220mM in absence of these ions. The removal of the Mg2+ ions thus introduces only a small perturbation of the electrostatic field around the complex, and the effect produced by their disappearance is not produced by a nonspecific electrostatic effect.
The solute conformations were hydrated by a box of water molecules in a pseudoequilibrated configuration. The water molecules farther than 9Å from any solute molecule were discarded. The simulation parameters (box dimensions, numbers of water molecules and atoms, type and number of counterions) are given in Table 1, along with a name for each trajectory.
| Table 1 MD simulations run on the S15 isolated molecule, and on the S15-RNA complex |
| Simulation name | Solute | Number of water molecules | Number of atoms | Box size (Å, Å, Å) | Duration (ns) | Ions | ||
|---|---|---|---|---|---|---|---|---|
| s15-pdb | S15 (1ab3) | 6,130 | 19,905 | 73, 66, 55 | 15 | 8 Cl− | ||
| s15-aria | S15 (2fkx) | 6,405 | 20,730 | 75, 67, 53 | 15 | 8 Cl− | ||
| s15-rna-mg | S15-RNA (1dk1) | 11,219 | 37,018 | 79, 67, 69 | 15 | 9 Mg2+, 31 Na+ | ||
| s15-rna | S15-RNA (1dk1) | 11,018 | 36,424 | 79, 66, 69 | 15 | 49 Na+ | ||
| The PDB ids of the structures providing the starting points are given in parentheses after the solute name. |
Heating and system equilibration were performed as follows: 10ps of heating up to 298K at constant volume was done while solute atom locations were restrained with a harmonic restraint of 25 kcal/mol Å2, followed by 5ps of MD at constant volume and four MD rounds of 2.5ps at constant pressure. The restraints on the solute atoms were progressively reduced to zero during these six rounds. No restraint was applied on the positions of Mg2+, Na+, or Cl−. The isothermal-isobaric unrestrained simulations were performed using the Berendsen algorithm 32 for temperature bath coupling. The equations of motion were integrated using a time step of 2fs. All covalent bonds involving hydrogen atoms were kept rigid using SHAKE 33.
The analysis was mainly performed with the program ptraj 29. The analysis of the conformational drift during the trajectories of the s15-rna complex and the monitoring of the S15, RNA, and complex electrostatic energies (data not shown) show that the s15-rna trajectories are stabilized in the 11–15-ns interval. The comparison of the s15-rna-mg and s15-rna trajectories was then based on mean and standard deviation values calculated on this interval. On the other hand, the trajectories s15-pdb and s15-aria were included integrally into the analysis, as the purpose of their analysis was to compare the conformational drift and the conformational landscapes of the structures 1ab3 and 2fkx.
During the simulation of the complex, the protein residues are numbered from 1 to 86 (corresponding to residues 101–186 of the 1dk1 file), the RNA residues are numbered from 87 to 143 (corresponding to residues 1–57 of the 1dk1 file), and the Mg2+ ions are numbered from 144 up to 152 (corresponding to residues 200–206 and 208–209 of the 1dk1 file). The residue numbers given in the analysis are those used in the MD simulation; the corresponding numbers of the 16S RNA, if they exist, are underlined in parentheses.
The calculation of S15 conformers under NMR restraints was performed using the simulated annealing protocol available in the CNS-ARIA protocol 34, which was recently used for the recalculation of NMR structures in the RECOORD database 35. The scripts were downloaded from http://www.ebi.ac.uk/msd-srv/docs/NMR/recoord/scripts.html. The nonbonded force field was PROLSQ, and the covalent geometry was parallhdg5.3.pro 36. The conformers of the NMR structure and the conformations of S15 observed in the trajectories were thus obtained in the presence of two different force fields, which alleviates the possible bias induced by the use of a unique force field.
The NMR restraints were 965 nuclear Overhauser effect (NOE) distance restraints and 38 dihedral restraints provided by H. Berglund and already applied in Berglund et al. 5. These restraints were supplemented by adding 76 hydrogen bond restraints (i, i+4) between amide hydrogens and carboxyl oxygens in α-helices. The definition of the α-helices was taken from the header of the 1ab3 file. The 25 best-energy conformers were selected from the 100 conformers calculated using the CNS-ARIA protocol. Among them, conformers displaying outlier orientations of the helices α1 and α4 were removed. The set of S15 conformers and the restraints used are deposited at the Protein Data Bank (id: 2fkx).
The NMR structures were analyzed using PROCHECK v.3.5.4 37 and WHATIF 8.1 38. The root mean-square deviation (RMSD) values between the conformer coordinates were calculated on the backbone atoms. The NOE distance restraints were examined in the NMR structures and in the trajectories by calculating, for each restraint, all distances ri separating the hydrogen pairs involved in the NOE. The mean distance value R was then determined as
![]() | (1) |
The variation of S15 tertiary structure was monitored by calculating the angles and the distances between the helix axes. The points belonging to each helix axis were determined as the middles of the backbone atom segments (N(i), N(i+2)), (Cα(i), Cα(i+2)), and (C′(i), C′(i+2)), where the residue numbers i vary along the secondary structure element. The angle between two axes was defined as the angle between the vectors linking the first and the last point of each axis. The distance between two helix axes was calculated as the minimum distance between the points defining the axes.
The comparison of 2fkx and 1ab3 conformers (Fig. 3) shows a large increase in 2fkx of the mean percentage of residues located in the core Ramachandran regions determined by PROCHECK (Table 2, part a). Similarly, the percentages of residues in the allowed and in the generously allowed regions decrease. Thus, the 2fkx conformers display better Ramachandran parameters than the 1ab3 conformers. The mean number of bad contacts observed by PROCHECK is also considerably reduced if the conformers are calculated with the CNS-ARIA scripts.
| Table 2 Structural quality of the NMR conformers of S15 from Berglund et al. 5 (PDB id: 1ab3) and recalculated using the CNS-ARIA scripts (PDB id: 2fkx) |
| (a) | PROCHECK Core (%) | PROCHECK Allowed (%) | PROCHECK Generously allowed (%) | ||
|---|---|---|---|---|---|
| 1ab3 | 59.4±2.9 | 32.6±3.2 | 7.5±1.9 | ||
| 2fkx | 82.8±3.6 | 14.3±2.9 | 1.5±1.5 | ||
| PROCHECK Nonallowed (%) | PROCHECK Mean numbers of bad contacts | |||
|---|---|---|---|---|
| 1ab3 | 0.5±0.6 | 18.1±3.0 | ||
| 2fkx | 1.4±1.1 | 0.1±0.3 | ||
| (b) | WHATIF Z score Second-generation packing quality | WHATIF Z score Ramachandran plot appearance | WHATIF Z score χ-1/χ-2 rotamer normality | ||
|---|---|---|---|---|---|
| 1ab3 | –4.7±0.3 | –6.4±0.4 | –4.6±0.3 | ||
| 2fkx | –0.6±0.7 | –3.6±0.7 | –0.1±0.6 | ||
| WHATIF Z score Backbone conformation | Mean numbers of bumps | |||
|---|---|---|---|---|
| 1ab3 | –7.6±1.0 | 104.1 | ||
| 2fkx | –3.4±1.4 | 12.0 | ||
| (c) | Number of conformers | Backbone RMSD (Å) (residues 24–70) | NOE violations | Dihedral violations | ||
|---|---|---|---|---|---|---|
| 1ab3 | 26 | 3.1 | 48 | 13 | ||
| 2fkx | 18 | 2.1 | 26 | 18 | ||
| The analysis was performed using PROCHECK v.3.5.4 37 and WHATIF 38. The violations are consistent violations, i.e., observed in >50% of the conformers. |
The WHATIF Z-scores (Table 2, part b) show a picture similar to the PROCHECK analysis. The mean Z-score values are improved, and all Z-scores are in the −4/4 range for the 2fkx conformers. The number of bumps detected by WHATIF is also much reduced in the 2fkx conformers. The improved quality of 2fkx comes from the optimization of the conformers in the torsion angle space 39,40 as well as from the use of a final refinement step in water 41.
The number of violated nuclear Overhauser effect (NOE) restraints is divided by two in the 2fkx conformers with respect to the 1ab3 conformers (Table 2, part c): this decrease comes from a much smaller number of long-range violations as well as a smaller number of (i, i+1) violations (Table 3). Nevertheless, the number of violated dihedral restraints slightly increases in 2fkx (Table 2, part c).
| Table 3 Numbers of NOE restraint violations in the NMR structures (1ab3 and 2fkx) and in the MD trajectories (s15-pdb and s15-aria) recorded on the isolated protein |
| NOE restraints | 1ab3 | 2fkx | s15-pdb | s15-aria | ||
|---|---|---|---|---|---|---|
| i,i | 9 | 10 | 15 | 21 | ||
| i,i+1 | 15 | 8 | 8 | 10 | ||
| i,i+2 | 0 | 0 | 1 | 1 | ||
| i,i+3 | 5 | 3 | 9 | 9 | ||
| |i–j|>3 | 19 | 6 | 37 | 27 | ||
| Total | 48 | 26 | 70 | 68 | ||
| The violations are consistent violations, i.e., observed in >50% of the conformers. The trajectory names are those given in Table 1. |
The relative orientations of the α-helices show different patterns in 1ab3 and 2fkx (Table 4). The mean values of the angles between α1 and α2 and between α3 and α4 are smaller in 2fkx than in 1ab3: the helices α1 and α4 are thus less parallel to the protein core formed by α2 and α3. The standard deviations of angles between the axes of α-helices in 1ab3 are much smaller than those observed in 2fkx. A similar trend is observed for the distances between the helix axes and particularly between the axes of α1 and α2. Thus, the helix α1 and, to a lesser extent, the helix α4 display a much wider range of orientations in 2fkx than in 1ab3. Because of the variability of the α1 and α4 orientations, larger RMSD values between the conformers of the complete protein are observed in the 2fkx than in the 1ab3 structure. But backbone RMSD values calculated by superimposing only the helices α2 and α3 of the protein (residues 24–70) are smaller for 2fkx than for 1ab3 (Table 2, part c).
| Table 4 Geometry of the S15 tertiary structure in the NMR structures (1ab3 and 2fkx), in the x-ray structures (1a32 and 1dk1), and in the structures of S15 in two conformations of the ribosome of E. coli (2avy and 2aw7) |
| Angle α1-α2 (degrees) | Angle α2-α3 (degrees) | Angle α3-α4 (degrees) | |||
|---|---|---|---|---|---|
| 1ab3 | 140.9±6.6 | 159.3±6.0 | 147.8±5.3 | ||
| 2fkx | 104.8±39.2 | 154.9±14.5 | 139.7±17.0 | ||
| 1a32 | 62.5 | 150.5 | 145.2 | ||
| 1dk1 | 140.0 | 153.0 | 149.1 | ||
| 2avy | 140.8 | 150.4 | 145.1 | ||
| 2aw7 | 141.8 | 150.1 | 143.1 | ||
| Distance α1-α2 (Å) | Distance α2-α3 (Å) | Distance α3-α4 (Å) | |||
| 1ab3 | 9.2±0.2 | 8.2±0.4 | 7.3±0.1 | ||
| 2fkx | 8.5±1.2 | 9.7±0.8 | 7.1±0.3 | ||
| 1a32 | 21.1 | 9.3 | 7.4 | ||
| 1dk1 | 8.6 | 9.3 | 7.6 | ||
| 2avy | 8.2 | 9.0 | 7.3 | ||
| 2aw7 | 8.2 | 9.0 | 7.2 | ||
The large variability observed for the α1 and α4 orientations arises from the sparsity of the NOE restraints network. Indeed, among the 91 long-range NOE restraints, 8 are applied between α1 and α2, 6 between α1 and α4, and 4 between α2 and α4. In α1, only one residue, Phe-14, located at the C-terminal end of the helix, is involved in interhelix interactions, whereas in α4, five residues, Arg-76, Tyr-77, Leu-80, Lys-83, and Leu-84, are involved in restraints with other helices. It is thus not surprising that the orientation of α4 with respect to the protein core is better defined than the orientation of α1. Furthermore, between α2 and α3, which have the best-defined relative orientation, 21 restraints exist, and they involve 5 residues in each helix.
The mean orientations of the α helices were also compared (Table 4) in the NMR (1ab3 and 2fkx) and in the x-ray (1a32 and 1dk1) structures of S15 as well as in the S15 structures present in the 3.5Å resolution x-ray structures (2avy and 2aw7) of the 70S ribosome of E. coli42. The α1-α2 and α2-α3 angles of the 2fkx conformers are between the values observed in the structures 1a32 and 1ab3. The angles and the distances between the α-helices in 1dk1, 2avy, and 2aw7 are close to those observed in 1ab3.
The backbone angles ϕ and ψ of residues 15–23, which form the loop connecting the helices α1 and α2, exhibit quite different mean values in the NMR structures (1ab3, 2fkx), as well as in the x-ray structures (1a32, 1dk1, 2avy, 2a7w) (Figure 4ab). Various conformations of the loop are thus observed in the different S15 structures.
The 2fkx conformers display a better structure quality and a lighter internal strain than those of the 1ab3 structure. Nevertheless, the convergence of the structure is worse: the orientations of the helices α1 and α4 with respect to the protein core are more variable in the 2fkx conformers. This feature is in agreement with the sparse network of long-range NOE restraints and with the x-ray structure of S15 (PDB id: 1a32), where the first helix is extended far apart from the other protein helices 6.
A name for each MD trajectory is given in Table 1, along with the main simulation parameters. In the two MD simulations recorded on the protein S15, the variation of the coordinate RMSD with respect to the starting point, shows different trends (Figure 5a). RMSD values larger than 7.0Å are observed for s15-pdb, whereas values of about 4.5Å are observed in the case of s15-aria. The protein conformation in s15-aria thus displays a better stability than that in s15-pdb, which is probably because of the lighter internal strain of the structure, observed in the previous section.
The large RMSD values observed in the trajectory s15-pdb mainly arise from a conformational drift observed for the loop connecting the helices α2 and α3 (Figure 5b). RMSD values larger than 2.0Å are also observed on the loop connecting α1 and α2, in the 2- to 4-ns interval (data not shown). The whole protein, except a few residues (residues 3–15 (α1) and residues 35–40 (α2)) consequently displays larger atomic fluctuations in s15-pdb than in s15-aria (Figure 5c). The differences in fluctuations are particularly large in the 15–25 and in the 42–52 regions, which correspond to the loops between helices α1 and α2 and between helices α2 and α3.
The angles between the helices α1 and α2 (Figure 5d) and α3 and α4 (Figure 5f) vary in the same range (90–140°) during both trajectories. This variability of the orientations of α1 and α4 with respect to the protein core is in agreement with the observations made above on 2fkx. On the other hand, the angle between α2 and α3 (Figure 5e) undergoes a transition in the nonstabilized part of the trajectory s15-pdb but is stable in the trajectory s15-aria.
The secondary structure elements have been monitored along the trajectories (data not shown) by calculating the mean and the standard deviation values of the backbone NH…O hydrogen bond distances in the α helices. In both trajectories, the helix α1 (residues 3–15) is not well defined; nor are the C-terminal part of α2 and the N-terminal part of α3. The C-terminal regions of α3 and α4 are less stable in s15-pdb than in s15-aria.
The numbers of violations for the (i, i+1), (i, i+2), (i, i+3) NOE restraints are similar (Table 3) for both trajectories. On the other hand, s15-aria displays 7 more violations for the intraresidue (i, i) NOEs, whereas 10 more violations are observed in s15-pdb for the long-range NOEs. The visualization of the violated NOEs on a contact map (Fig. 6) demonstrates that very long-range NOE restraints (i, j) such that |i–j|>30 are observed in s15-pdb but not in s15-aria. On the other hand, most of the (i, j) NOE violations observed in s15-aria are such that |i–j|<10. This distribution of the violations is in agreement with the larger conformational drift observed in s15-pdb than in s15-aria.
The observations quoted above show that the trajectory starting from a 2fkx conformer is more stable than the one recorded from a 1ab3 conformer. This is probably because of the calculation of the 2fkx structure in the torsion angle space as well as the use of a final refinement step in water, as it was shown 35,39,41 that these two ingredients improve the quality of the structure and thus the stability in MD simulation. On the other hand, the better convergence of 1ab3, which is the sign of a larger internal strain of the structure, is correlated to a larger conformational drift in the s15-pdb simulation 43.
The variable orientations observed for the helices in S15 may arise from the bundle helix structure of the protein S15. Indeed, in p8MTPC144, a protein consisting of three α-helices and having a similar tertiary structure, the third helix displays internal motion and large displacements in MD simulations 45, and this behavior is in agreement with dipolar coupling measurements 46.
The mobility of the partners of S15-RNA interaction was analyzed during the trajectories s15-rna-mg and s15-rna, by calculating separately on S15 and on RNA the atomic fluctuations by residues from the 11–15-ns interval of the trajectories. In the protein, a slight increase of fluctuations is observed in absence of Mg2+ (Figure 7a) for residues 37–50 in helices α2 and α3 and in the loop connecting them.
In absence of Mg2+, an increase of flexibility of RNA is observed in four regions (Figure 7b): residues 92–100 (588–650) in the TL1 tetraloop, residues 102–104 (652–654) and 133–134 (748–749) near the three-way junction, residue 115 (665) close to the G-U/G-C motif, and residues 121–123 in the TL2 tetraloop. Among these residues, A103 (653), A115 (665), and C133 (748) are not involved in canonical basepairs in the structure 1dk1. The increase of RNA fluctuations is in agreement with the positions of the Mg2+ ions in the crystal, close to the G-U/G-C motif and to the three-way junction motif. The increases of flexibility observed in protein helices α2 and α3, in the α2-α3 loop (Figure 7a), and in the TL2 tetraloop (Figure 7b), are consistent because these regions interact.
In the presence of Mg2+ ions, the RMSD of S15 atomic coordinates from the initial conformation converges to a plateau of 3Å after 11ns (Figure 7c). The trajectory recorded in the absence of Mg2+ ions reaches the same value after 3ns (Figure 7c). The variation of the coordinate RMSD of the protein S15 in the presence of Mg2+ displays a jump in the 7.8- to 11.5-ns range (Figure 7a) before attaining the plateau value. The conformations corresponding to the minimum RMSD values in this interval (7881ps: 2.31Å, 9895ps: 2.41Å) were compared (Figure 8a) with the conformations corresponding to the maximum RMSD values (8837ps: 3.11Å, 11,347ps: 3.29Å). The conformation of the α1-α2 loop changes, and residues 15–19 undergo the largest displacement (Figure 8b). Residues 3–7 move simultaneously (Figure 8c), and the backbone hydrogen bond between Lys-4 and Gln-8 is disrupted. The simultaneous variation of conformation for residues 3–7 and 15–19 is consistent with a mechanical link between the loop connecting helices α1 and α2 and helix α1.
Because the α1-α2 loop is located at the intermolecular interface of the S15-RNA complexes in the crystal, the loop conformation is probably disturbed by crystal packing effects. The loop drifts during all MD trajectories, recorded on the S15-RNA complex and on S15, as similar mean ϕ and ψ values are calculated from the whole MD trajectories in the loop connecting α1 and α2 (Figure 4cd), whereas different values are observed in the different NMR and x-ray structures (Figure 4ab). Transitions of the ϕ and ψ values are thus observed (Fig. 9) in residues located in the α1-α2 loop, in the presence or in the absence of Mg2+ ions.
The angles between the S15 helices and the distances separating them vary <10° and 1Å, respectively (data not shown), during the trajectories s15-rna-mg and s15-rna: the presence of RNA freezes the variability of the helices’ orientations, observed for the isolated protein. Nevertheless, the angle between α2 and α3 is increased by 6° if the Mg2+ ions are removed, and because the helices are antiparallel, they are thus further apart in the presence of Mg2+ ions. Also, when the Mg2+ ions are removed, the distance between α1 and α2 is reduced by 1Å after 10ns.
The variation of the coordinate RMSD of the RNA in presence of Mg2+ displays two jumps at around 4 and 6ns (Figure 7b). The conformations exhibiting the maximum RMSD values (3812ps: 3.82Å, 6212ps: 4.22Å) were compared to conformations exhibiting the minimum RMSD values (3261ps: 2.20Å, 4782ps: 2.01Å, 6347ps: 2.28Å). The 3′ and 5′ extremities spanning the 87–90 (584–586) and 140–143 (755–757) residues exhibit different conformations.
The interphosphate distances between phosphate groups facing each other in the RNA backbone were calculated on the 11–15-ns interval of the trajectories s15-rna-mg and s15-rna. All distance differences are smaller than 1.5Å, except the four following ones. The pairs of residues A113(663)-G127(742), G114(664)-G126(741), G116(666)-U125(740), and U119-G122 show variations of distance of 2.0Å, 2.1Å, 2.0Å, and 4.6Å. These residues are located in the G-U/G-C region and in tetraloop TL2. The phosphate groups thus move apart in the absence of Mg2+ ions, which is in agreement with the increase of internal flexibility in this RNA region (Figure 7b).
The values of the RNA phosphodiester angles were compared in the presence and in the absence of Mg2+ ions, by calculating their mean and standard deviation values in the 11–15-ns interval (Fig. 10). Few angles display a difference of mean values significant with respect to the standard deviation values. These differences concern, for the angle α, 11 residues, for the angle β, 1 residue, for the angle γ, 7 residues, for the angle δ, 6 residues, for the angle ɛ, 6 residues, and for the angle ζ, 4 residues. These residues cluster in the UUCG TL2 tetraloop, in the three-way junction, and in the helix H22, close to the positions of the Mg2+ ions. These phosphodiester angle variations are thus in agreement with the stabilization of the RNA structure by the Mg2+ ions 2. But the hydrogen bond parameters (data not shown) involved in a base triple among G104(654), G137(752), and C139(754) are not modified between s15-rna and s15-rna-mg.
Comparison of the mean angle values from the trajectories with the corresponding angles measured in 1dk1 shows (Fig. 10) that, for the majority of the residues, the trajectories’ angles are in a similar angular region than the values observed in the structure. The x-ray angles were compared to the interval of variation of the trajectory angles, and the residues for which the x-ray angle is more than 90° outside of the interval are the following ones: G92 (588), U102 (652), G107(657), A113(663), G114(664), A115(665), U120, G122, G123, U125(740), G127(742), C133(748), C134(749), G135(750), A138(753), and C143.
The angle values calculated on the trajectories and measured on the 1dk1 structure were also compared to the rotameric classification introduced recently 47 to describe the RNA backbone. The mean absolute difference value between the RNA angles and the rotameric angles was calculated over the 53 sugar-to-sugar angle suites present in the RNA. For each pair of angle values p and q, the difference was the smallest absolute value of the three numbers: p–q, p-q+360° and p–q–360°. The RX structure of S15-RNA displays 18 suites with means difference larger than 20°, and the maximum mean difference observed is 79.9°. On the other hand, the s15-rna-mg and s15-rna trajectories have 10 and 13 such suites, and their maximum mean differences are 46.4° and 73.8° The conformations explored during the 11–15-ns trajectory intervals thus seem to be in reasonable agreement with the analysis of rotamers performed on 132 RNA crystallographic structures 47.
The simulation of UUCG loops was used as a test case in MD simulations 48,49; the first NMR structure of the UUCG loop 50 was incorrect and was corrected 51 by discovering a misassignment. The use of combined locally enhanced sampling and Particle Mesh Ewald was shown 49 to allow the simulation to find the correct structure starting from the incorrect one. However, several simulations 48 of the two UUCG structures concluded with the formation of the hydrogen bond U1(O2′)-U2(O5′) and not U1(O2′)-G4(O6).
The conformation of the two UUCG tetraloops was compared (Fig. 11) to experimental structures of UUCG (1dk1, 1f7y, 1hlx, 1i6u, 1kuq). The structure 1hlx is the NMR structure of the P1 helix from group I self-splicing introns 51, in which only one UUCG loop is present. The others are the S15-RNA 3,4 and S18-RNA 52 x-ray structures in which a three-way-junction RNA is present, with both tetraloops TL1 and TL2. The values of the phosphodiester backbone angles are similar in structures and in MD trajectories. The angles displaying the largest variations between s15-rna-mg and s15-rna also vary among the experimental structures.
A set of hydrogen bonds already analyzed in a MD simulation of UUCG 48 was compared in the MD trajectories and in x-ray and NMR UUCG structures (Table 5). The comparison between TL1 and TL2 inside the MD trajectories shows that in s15-rna-mg, over nine distances, five show a difference larger than 1Å, whereas in s15-rna, only two distances show a difference larger than 1Å. This reduction in the number of differences is in agreement with the homogenization of the environment of the UUCG loops obtained by removing the Mg2+ions. The loop TL1 shows the same pattern of formed hydrogen bonds in both MD trajectories, but the percentage of hydrogen bond formation generally decreases in absence of Mg2+.
| Table 5 Parameters of hydrogen bonds in UUCG loops calculated on the 11–15-ns time interval of the trajectories s15-rna-mg and s15-rna and on PDB structures (1dk1 3, 1f7y 4, 1hlx 51, 1i6u 52, 1kuq) |
| Hydrogen bond | Distance (Å) | Angle (degrees) | Formation (%) | ||
|---|---|---|---|---|---|
| s15-rna-mg | |||||
| TL1 U1(O2) G4(N2/H21) | 2.4±0.4 | 139.5±12.8 | 45.9 | ||
| TL1 U2(O2P) C3(N4/H41) | 3.8±0.4 | 37.8±6.7 | 0 | ||
| TL1 U1(O2) G4(N2/H22) | 4.0±0.4 | 25.6±7.7 | 0 | ||
| TL1 U2(O2P) C3(N4/H42) | 2.1±0.4 | 155.3±13.2 | 77.3 | ||
| TL1 U1(HO′2) G4(O6) | 3.9±0.9 | 57.3±42.4 | 12.9 | ||
| TL1 U1(HO′2) U2(O5′) | 2.5±0.72 | 40.1±37.5 | 46.4 | ||
| TL1 U1(H3) G4(O6) | 6.3±0.4 | 74.1±6.1 | 0 | ||
| TL1 U1(O2) G4(N1) | 1.9±0.2 | 153.3±12.1 | 90.1 | ||
| TL1 U1(O4) G4(N1) | 6.3±0.3 | 147.4±10.3 | 0 | ||
| TL2 U1(O2) G4(N2/H21) | 6.6±0.7 | 72.4±6.8 | 0 | ||
| TL2 U2(O2P) C3(N4/H41) | 3.7±0.2 | 37.8±5.2 | 0 | ||
| TL2 U1(O2) G4(N2/H22) | 7.1±0.7 | 41.6±7.3 | 0 | ||
| TL2 U2(O2P) C3(N4/H42) | 2.0±0.2 | 162.2±8.7 | 88.5 | ||
| TL2 U1(HO′2) G4(O6) | 6.1±0.9 | 53.3±41.7 | 0 | ||
| TL2 U1(HO′2) U2(O5′) | 3.1±0.6 | 59.7±28.9 | 1.2 | ||
| TL2 U1(H3) G4(O6) | 6.3±0.4 | 74.1±6.1 | 0 | ||
| TL2 U1(O2) G4(N1) | 4.9±0.6 | 62.1±7.1 | 0 | ||
| TL2 U1(O4) G4(N1) | 8.0±0.4 | 45.4±5.8 | 0 | ||
| s15-rna | |||||
| TL1 U1(O2) G4(N2/H21) | 2.1±0.3 | 148.9±13.7 | 75.4 | ||
| TL1 U2(O2P) C3(N4/H41) | 3.7±0.4 | 38.2±6.6 | 0 | ||
| TL1 U1(O2) G4(N2/H22) | 3.8±0.3 | 31.0±9.3 | 0 | ||
| TL1 U2(O2P) C3(N4/H42) | 2.0±0.4 | 159.6±11.1 | 86.8 | ||
| TL1 U1(HO′2) G4(O6) | 4.5±0.6 | 49.7±19.1 | 2.0 | ||
| TL1 U1(HO′2) U2(O5′) | 2.4±0.4 | 34.4±17.2 | 26.9 | ||
| TL1 U1(H3) G4(O6) | 6.6±0.5 | 79.1±6.5 | 0 | ||
| TL1 U1(O2) G4(N1) | 2.2±0.4 | 144.7±12.8 | 62.5 | ||
| TL1 U1(O4) G4(N1) | 6.6±0.4 | 136.8±11.3 | 0 | ||
| TL2 U1(O2) G4(N2/H21) | 2.6±0.5 | 131.6±11.8 | 23.7 | ||
| TL2 U2(O2P) C3(N4/H41) | 4.0±0.8 | 39.7±6.9 | 0 | ||
| TL2 U1(O2) G4(N2/H22) | 4.3±0.4 | 20.0±7.2 | 0 | ||
| TL2 U2(O2P) C3(N4/H42) | 2.3±0.8 | 162.4±9.8 | 67.2 | ||
| TL2 U1(HO′2) G4(O6) | 7.5±0.5 | 31.5±14.4 | 0 | ||
| TL2 U1(HO′2) U2(O5′) | 2.4±0.3 | 30.8±16.5 | 28.1 | ||
| TL2 U1(H3) G4(O6) | 3.8±0.3 | 132.7±9.3 | 0 | ||
| TL2 U1(O2) G4(N1) | 1.9±0.1 | 158.6±10.6 | 96.8 | ||
| TL2 U1(O4) G4(N1) | 5.8±0.3 | 149.5±8.5 | 0 | ||
| 1dk1 (Å) | 1f7y (Å) | 1i6u (Å) | 1kuq (Å) | |||
|---|---|---|---|---|---|---|
| TL1 U1(O2) G4(N2) | 3.3 | 3.3 | 3.7 | 2.8 | ||
| TL1 U2(O2P) C3(N4) | 2.9 | 3.0 | 2.9 | 3.4 | ||
| TL1 U1(O2′) G4(O6) | 2.5 | 2.4 | 2.8 | 2.5 | ||
| TL1 U1(O2′) U2(O5′) | 3.7 | 3.7 | 3.6 | 4.4 | ||
| TL1 U1(N3) G4(O6) | 5.9 | 5.9 | 5.8 | 6.3 | ||
| TL1 U1(O2) G4(N1) | 2.8 | 2.9 | 2.9 | 2.8 | ||
| TL1 U1(O4) G4(N1) | 7.0 | 7.0 | 6.9 | 7.1 | ||
| TL2 U1(O2) G4(N2) | 4.8 | 3.2 | 3.5 | 7.5 | ||
| TL2 U2(O2P) C3(N4) | 2.6 | 2.7 | 2.9 | 7.6 | ||
| TL2 U1(O2′) G4(O6) | 4.6 | 3.0 | 2.6 | 3.7 | ||
| TL2 U1(O2′) U2(O5′) | 3.8 | 3.9 | 5.0 | 4.0 | ||
| TL2 U1(N3) G4(O6) | 3.2 | 6.5 | 5.9 | 5.5 | ||
| TL2 U1(O2) G4(N1) | 3.0 | 3.2 | 3.0 | 5.3 | ||
| TL2 U1(O4) G4(N1) | 7.0 | 7.3 | 7.1 | 9.7 | ||
| 1hlx | Distance (Å) | Angle (degrees) | Formation (%) | ||
|---|---|---|---|---|---|
| TL1 U1(O2) G4(N2/H21) | 3.9±0. | 27.6±3.2 | 0 | ||
| TL1 U2(O2P) C3(N4/H41) | 3.8±1.1 | 131.5±5.5 | 0 | ||
| TL1 U1(O2) G4(N2/H22) | 2.1±0.2 | 145.5±5.8 | 65 | ||
| TL1 U2(O2P)C3(N4/H42) | 5.2±1.1 | 40.5±2.9 | 0 | ||
| TL1 U1(HO′2) G4(O6) | 2.7±0.3 | 72.5±18.9 | 15 | ||
| TL1 U1(HO′2) U2(O5′) | 3.0±0.4 | 84.0±26.8 | 10 | ||
| TL1 U1(N3) G4(O6) | 6.6±0.2 | 61.3±2.1 | 0 | ||
| TL1 U1(O2) G4(N1) | 2.0±0.1 | 153.7±6.8 | 90 | ||
| TL1 U1(O4) G4(N1) | 6.5±0.1 | 158.8±7.0 | 0 | ||
| The hydrogen bond mean distance is given as well as the donor-hydrogen-acceptor mean angle when available. The formation is the percentage of trajectory time or of NMR conformers in which the mean distance is <2.2Å. |
The comparison between the MD trajectories and the structures gives the following results. The pattern of formed hydrogen bonds in TL1 is qualitatively similar to those obtained for TL1 in x-ray and NMR structures. In TL2, the distances vary much between the x-ray structures, and consequently, the formation of the hydrogen bond observed in the MD trajectories is always in agreement with a least one x-ray structure. The hydrogen bond between U1(O2′) and G4(O6) is partially formed in TL1 but not in TL2 for s15-rna-mg. In all x-ray structures, this hydrogen bond is formed in TL1, but in the NMR structure, it is partially formed.
The hydrogen bond between U1(O2′) and U2(O5′) is stronger than the one between U1(O2′) and G4(O6) for TL1 and TL2 in both s15-rna trajectories. This observation is in agreement with that of Miller and Kollman 48 but contradicts the observations made on the x-ray structures, as already pointed out by Ennifar et al. 4. Nevertheless, the analysis of the NMR structure 1hlx shows that the hydrogen bonds between U1(O2′) and G4(O6) and U1(O2′) and U2(O5′) are both partially formed in the NMR structure, with formation percentages of the same order of value. In all x-ray structures, the UUCG loop is included in a larger RNA sequence. On the other hand, the only structure of an isolated UUCG loop is the NMR structure 1hlx. The nonavailability of a crystallized structure for this molecule may be the sign of a large internal flexibility, which is also observed in MD trajectories.
The previous analysis of the hydrogen bonds inside loops UUCG showed that the removal of Mg2+ ions induces a reorganization of the hydrogen bond network, but according to the limitation of the modeling of ion electrostatics, it is not possible to state that this behavior is specific of Mg2+ ions. One should also notice that 1), the Mg2+ ions observed in the S15-RNA crystal may not all have to be present, 2), in the NMR sample of the P1 helix from group I self-splicing introns, no Mg2+ ions were present, but other ions were present (G. Varani, University of Washington, personal communication, 2006). The comparison among the MD trajectories and the x-ray and NMR structures reveals that floppier hydrogen bonds are generally observed in the NMR structure, which agrees with the observations made during the MD trajectories. Furthermore, the large heterogeneity of distances in TL2 agrees with the larger number of disrupted hydrogen bonds in TL2 during the MD trajectories. The overall analysis performed on the phosphodiester backbone angles and on the hydrogen bonds in TL1 and TL2 shows that the overall topology of the simulated TLs did not change much with respect to the crystallographic structures.
The relative positions of the three RNA helices forming the three-way junction were monitored (Fig. 12) by calculating the two angles formed by the atom triplets (107-C6, 114-N2, 142-O2) and (107-C6, 114-N2, 98-N1). The first angle describes the relative orientation of helices H22 and H20, and the second angle describes the relative orientation of helices H22 and H21. Both angles are slightly larger in the absence of Mg2+ ions, which means that the helices have a tendency to go apart, in agreement with the observations of Batey and Williamson 2.
The positions of the nine crystallographic Mg2+ ions (Fig. 13) were analyzed along the full trajectory s15-rna-mg by determining, for each snapshot, the atoms forming the coordination cage of the Mg2+ ions. The six oxygens closest to each Mg2+ were determined for each snapshot. Table 6 describes the coordination partners of each Mg2+ion along with their residence times.
| Table 6 Coordination numbers of the Mg2+ ions and partner names, with residence times (ns) in parentheses |
| Mg2+ ion | Coordination number | Coordination partners | ||
|---|---|---|---|---|
| 144 | 5 | OW-8570(15.0) OW-8572(15.0) OW-8577 (15.0) | ||
| OW-8584 (15.0) OW-8596 (15.0) | ||||
| 145 | 6 | O2P-112(13.3) OW-4913(15) OW-5012 (15.0) | ||
| OW-5076 (15.0) OW-5151 (15.0) OW-5551 (15.0) | ||||
| 146 | 6 | O6-107(13.3) O6-135(13.3) OW-8017 (15.0) | ||
| OW-8035 (15.0) OW-7978 (15.0) OW-8032 (15.0) | ||||
| 147 | 6 | OW-6061(14.9) OW-6063(14.9) OW-6064 (14.9) | ||
| OW-6171 (14.9) OW-6174 (14.9) OW-8595 (14.9) | ||||
| 148 | 5 | OW-5524(14.4) OW-5568(14.5) OW-5986 (14.6) | ||
| OW-5992 (14.4) OW-6020 (14.5) | ||||
| 149 | 6 | OW-3038(14.9) OW-3048(14.8) OW-5541 (14.9) | ||
| OW-5545 (14.8) OW-5564 (14.8) OW-5565 (14.8) | ||||