INTRODUCTION

TOP ABSTRACT INTRODUCTION METHODS RESULTS AND DISCUSSION CONCLUSIONS REFERENCES

The interactions that stabilize complexes between ribonucleic acids (RNA) and proteins are of great interest for understanding splicing, translation, and other RNA-mediated processes. At present, there is relatively little information available concerning RNA/protein interactions, especially for the single-stranded RNA/protein systems (Antson, 2000). The relative lack of information concerning RNA/protein complex structure and stability results, in part, from the technical difficulty in determining the structures of RNA/protein complexes. RNA/protein complexes are often difficult to crystallize, and determining RNA structures by NMR is challenging due to poor spectral dispersion and problems with isotopic enrichment. For these reasons, only a few structures of RNA/protein complexes become available each year (Cusack, 1999). Even for those RNA molecules for which structural information is available, identifying the interactions that stabilize RNA/protein complexes is challenging because of the complexity of interactions that contribute to complex stability (Saenger, 1983; Nolan et al., 1999). In many cases, RNA duplexes are not involved in the interactions that stabilize RNA/protein complexes because the RNA bases in helical regions are largely unavailable for hydrogen bond formation. Thus, many of the specific contacts that stabilize RNA/protein complexes involve ribonucleotides in single-stranded or loop regions of RNA (Oubridge et al., 1994; Price et al., 1998; Honda et al., 1999).

As a component of the spliceosomal U2 small nucleic ribonucleoprotein particle (RNP), the U2B" protein interacts with U2 snRNA hairpin IV (Fig. 1). Recently, an x-ray structure of the U2 snRNA hairpin IV/U2B" complex was reported (Price et al., 1998). The structure of this complex revealed U2B" interacted with U2 hairpin IV in a similar manner to that previously observed to occur for the U1A protein in complex with U1 snRNA hairpin II. The U2A'/U2B"/U2 snRNA forms a sandwich-like ternary complex in which U2A' interacts strongly with U2B", and interacts weakly with U2 snRNA hairpin IV. This strong interaction between U2A'/U2B" alters the conformation of U2B", allowing U2B" to bind U2 snRNA (Scherly et al., 1990). The U2 snRNA hairpin IV includes a stem consisting of five Watson-Crick basepairs and one non-canonical U-U basepair (Gmeiner and Walberer, 2000; Fig. 1). The U-U basepair occurs at the base of an 11-nucleotide loop. Although the U2A' protein is essential for U2B" binding to U2 snRNA hairpin IV, the direct interaction between U2A' and U2 snRNA hairpin IV is very weak, and is limited to the stem region of the RNA in the x-ray structure (Price et al., 1998; Boelens et al., 1990).

View larger version (16K): [in this window] [in a new window]   FIGURE 1   Secondary structure and sequence numbering of the U2 snRNA hairpin IV/U2B" protein complex. The conserved RNP1 and RNP2 sequences for the U2B" protein are shown in red. The filled circles indicate sites of hydrogen bond formation between protein and RNA.

Understanding the physical basis for formation of the U2 snRNA hairpin IV and U2B" complex in vivo requires detailed information about complex structure and stability under fully-hydrated conditions. X-ray structures of RNA/protein complexes provide useful structural information; however, x-ray structures may differ from fully hydrated structures. Furthermore, the crystalline state cannot fully reflect the contribution of mobile water molecules to RNA/protein complex stability. Although evidence for some water molecules can be seen in the experimental electron density map, water molecules that are in fast exchange on the surface of RNA/protein complexes, and those present in cavities of the complex, are not well accounted for in x-ray studies (Schwabe, 1997). Molecular dynamics simulations are one of the best methods available to investigate the dynamic behavior of the RNA/protein interface.

In the present study, the interaction between U2B" protein and U2 snRNA hairpin IV in solution has been investigated using molecular dynamics (MD) simulations to identify the residues of U2B" that provide the greatest contribution to complex stability (Fig. 2). These MD simulations also provide insight into the role of ordered water in complex stability. To perform MD on a system as large as the hydrated U2B"/U2 snRNA hairpin IV complex, the long-range electrostatic energy of the periodic box was calculated using the particle mesh Ewald (PME) method (Darden et al., 1993; Essman et al., 1995). The MD trajectory and the structure in solution for the complex were analyzed, and the mutant free energy was calculated to determine the contribution from each amino acid in the protein to the stability of the RNA/protein complex. Finally, the location and the residence times of waters residing at the interface of the RNA/protein complex were identified, and the stabilizing interaction for those bridging water molecules that play a significant role in the stability of the RNA/protein complex have been characterized using ab initio calculation. The results indicate that the U2B"/U2 snRNA hairpin IV complex has a novel interaction surface stabilized mainly by electrostatic interactions, and that bridging water molecules contribute significantly to complex stability.

View larger version (25K): [in this window] [in a new window]   FIGURE 2   Ribbon diagram of the U2 snRNA hairpin IV/U2B" protein complex. The protein is shown in purple and the RNA in green. Residues that link 2 and 3 of U2B" (e.g., K50 as shown) protrude into the interior of the RNA loop.

	METHODS

TOP ABSTRACT INTRODUCTION METHODS RESULTS AND DISCUSSION CONCLUSIONS REFERENCES

The structure of the U2 snRNA hairpin IV/U2B" protein was obtained from the x-ray crystallographic coordinates at 2.38 Å resolution (Price et al., 1998) in the Brookhaven protein database (identifier 1A9N). Fourteen K⁺ ions were added to the RNA/protein complex by using a Coulombic grid potential using the LEAP module of Amber 6.0 to maintain the system at charge neutrality (Case et al., 2000). The system was then solvated using TIP3P water (Jorgensen et al., 1983) in a periodic box. The RNA/protein complex had at least 10 Å buffer in every direction of the box to permit substantial fluctuations of the conformation during the course of the MD simulation. The MD simulation was performed using the AMBER suite of programs, version 6.0 (Case et al., 2000). Energy minimization (500 steps steepest descent followed by 500 steps conjugate gradient) was undertaken before initiating the MD simulation to remove initial steric clashes. The complete interaction energy was calculated. The constant pressure MD simulation was calculated using anisotropic diagonal position scaling. The time step used was 0.002 ps. The temperature of the system was increased gradually from 100 K to 300 K with 20-ps NPT reassemble. The target pressure was 1 atm. The Berendsen algorithm (Berendsen et al., 1984) was used with a scaling factor time constant of 0.2. The Lennard-Jones cutoff value used was 8 Å. SHAKE constraints were applied to all bonds involving hydrogen atoms. The PME (Darden et al., 1993; Essman et al., 1995) was used in all calculations to consider the long-range electrostatic energy with grid spacing around 1.0 Å. The size of the charge grid was chosen to be a product of powers of 2, 3, and 5 for each dimension so that the fast Fourier transform could be applied to increase the speed of the calculation of the reciprocal sum. Finally, the 2200-ps MD simulation was run under the same conditions as the equilibration procedure. The density of the system was maintained near 1 g/cm³. In all calculations the AMBER94 force field (Cornell et al., 1995, 1996) was used. The center of mass translation was removed periodically (at every restart or 100 ps) during the production dynamics (Cheatham and Kollman, 1997) to avoid energy drains (Harvey et al., 1998; Chiu et al., 2000). The analysis for the trajectory was carried out using Carnal, Ptraj, and in-house scripts (shell, awk, perl) similar to our previous work (Guo and Gmeiner, 2000; Guo et al., 2000).

Mutations were introduced by using hydrogen with a C-H distance of 1.09 Å to replace the side chain of other residues (XxxGly). A detailed description of the method can be found on studies of protein/peptide systems (Massova and Kollman, 1999). The mutant binding free energy was calculated using these same methods. The stable and well-equilibrated final 1800 ps of the 2200-ps trajectory was used for the mutant analysis with snapshots each 5 ps resulting in highly converged solutions with low standard deviations. A total of 360-snapshot conformations were used for every mutant residue system. The internal energy was calculated using the ANAL module in AMBER6. The solvent-accessible surface area was obtained from the program MSMS with a probe radius of 1.4 Å (Sanner et al., 1996). The solvation energy was calculated by using the generalized Born method. (Still et al., 1990; Srinivasan et al., 1998). Previous studies have found that the relative entropic contribution was negligible, and the results are in good agreement with experimental measurement for the mutant free energy between wild-type and mutant complexes in other systems (Reyes and Kollman, 2000a; Massova and Kollman, 2000). In addition, the computational cost for normal mode analysis (NMA) is expensive, considering all snapshots for large systems. Therefore, we did not consider the entropic contribution in the mutant free energy calculations.

The ab initio calculations were carried out using the Gaussian98 program (Frisch et al., 1998). One specific bridging water molecule (water 624), which interacts with cytosine 10 and threonine 89, was calculated as a model system. The starting geometry was obtained from the minimization of the MD averaged structure over the residual periods in the trajectory. The fully geometry optimizations were finished on the Hartree-Fock level with the 6-31G* basis set for the complex, and for each molecule separately. The Becke's three parameter hybrid function using the LYP correlation function (B3LYP) was evaluated for the density functional theory (DFT) calculation (Becke, 1993; Lee et al., 1988). The interaction energy was corrected for basis set superposition errors (BSSE) using the full counterpoise correction (Boys and Bernardi, 1970). The vibrational zero-point energies (ZPE) correction was considered from the frequency calculations without scaling. In the calculation of interaction energy the deformation energy was also taken into account.

	RESULTS AND DISCUSSION

TOP ABSTRACT INTRODUCTION METHODS RESULTS AND DISCUSSION CONCLUSIONS REFERENCES

Stability of the trajectory

Root-mean-square deviations (RMSD) from the starting structure coordinates during the MD simulation were used to assess the stability and the reliability of the simulation. RMSD from the starting structure over the course of the MD simulation is shown in Fig. 3 for the U2 snRNA hairpin IV/U2B" complex, and also separately for the RNA and protein components of the complex. Both the RNA and protein components of the complex underwent structural fluctuations during the first 500 ps of the trajectory, although the amplitude of the fluctuations was larger for the RNA component. Equilibration of the structure for the complex was achieved after the initial 750 ps of the trajectory. The all-atom RMSD of the U2 snRNA hairpin IV averaged over the remainder of the trajectory was 2.0 ± 0.3 Å, while that for the U2B" protein was 2.0 ± 0.5 Å. The RMSD of the whole system was 2.4 ± 0.60 Å. The RMSD of the complex during the MD simulation indicates that the x-ray crystal structure did not undergo a large change in aqueous solution.

View larger version (16K): [in this window] [in a new window]   FIGURE 3   RMSD, relative to the initial coordinates, during the 2200-ps molecular dynamics trajectory. RMSD for the U2B" protein is shown in the top panel, while those for the U2 snRNA hairpin IV and whole system are shown in the middle and bottom panels, respectively. The trajectory is stable after the initial 750 ps of MD.

U2 snRNA hairpin IV

The average values for the pseudorotation phase angle of the ribose sugars (P), the glycosidic bond torsional angle , and the torsional angles in the phosphodiester backbone for U2 snRNA hairpin IV over the course of the trajectory were, with several notable exceptions, similar to those characteristic of A-form helical geometry. The pseudorotation phase angles P adopted by the ribonucleotides in the stem region of hairpin IV (C0-U5 and U17-G22), indicated that these ribose sugars uniformly adopted C3'-endo conformations throughout the trajectory. The ribonucleotides in the 5' portion of the loop region of hairpin IV alternated between C3'-endo and C2'-endo sugar puckers with U7, G9, and A11 adopting C3'-endo sugar puckers and U8, C10, and G12 adopting C2'-endo sugar puckers. U13 also adopted a C2'-endo sugar pucker, breaking the strict alternating pattern, and A14 also adopted a C2'-endo sugar pucker. The UGCA sequence is common to both U1 snRNA hairpin II and U2 snRNA hairpin IV, and these residues were found in the x-ray structures (Price et al., 1998) of both complexes to be engaged in nearly identical interactions with residues conserved between the U1A and U2B" proteins.

The loop region of U2 snRNA hairpin IV was stabilized by a number of intramolecular hydrogen-bonding and base-stacking interactions. The 5' portion of hairpin IV (A6-G12) was stabilized by three pairs of base-stacking interactions. A6 stacked efficiently on U7 and U8 stacked efficiently on G9, but little base overlap occurred between U7 and U8. A single intramolecular hydrogen bond, between U7 O4 and G9 NH₂, stabilized this region of the loop structure. The torsion angles for A6-G9 were near canonical A-form values except for U8, which adopted a C2'-endo sugar pucker and anomalously high values of , , and . A6 in this region also adopted an anomalously low value of . The two purines, A11 and G12, also stacked efficiently on one another, and this region of the loop was stabilized by a bifurcated hydrogen bond from the 2'-OH of A11 with both O4' and O5' of G12. This region of the loop was also stabilized by an intraresidue hydrogen bond between G12 NH₂ and O1P oxygen of this residue. C10 stacked neither with the A11/G12 purine stack nor the U8/G9 stack, and made no intramolecular hydrogen bonds. Rather, C10 was involved in a number of contacts with the U2B" protein that stabilized the complex (see below).

The pairwise stacking arrangements in the 5' portion of the loop region (A6-G12) initiated the reverse in directionality of the phosphodiester backbone of the RNA, and the A11/G12 purine stack was roughly perpendicular to the axis of the helical stem region. The four nucleotides in the 3'-portion of the loop formed a "step-ladder structure" that was stabilized by intramolecular base-stacking between A14, C15, and C16. The RNA conformation in this region was also stabilized by intramolecular hydrogen bond formation between O2P of A14 and the N4 amino group of C16, and between A15 O1P and the 2'-OH of U13. Base-stacking did not continue from C16 into the stem region, rather U17, the first residue in the 3' portion of the stem, region stacked with A6/U7 in the 5' portion of the loop region. U13 and C16 adopted C2'-endo sugar puckers and anomalously high values of  and . A14 adopted a syn orientation about  and the rotation about C4'-C5' bond for A14 was in the ap range, while other nucleotides in the loop were in the +sc range. Obviously, O5' and the phosphate of A14 were rotated away from the ribose ring. This opened up the nucleotide exposing its sugar, base, and phosphate groups to the surface of the protein to facilitate recognition and binding contacts (Saenger, 1983). Other torsion angles in this region were similar to A-form values.

The non-Watson-Crick U5-U17 basepair (Gmeiner and Walberer, 2000) in the U2B"/U2 snRNA hairpin IV complex closes the RNA loop. This U-U basepair is one of the points of dissimilarity with the U1A/RNA complex, which includes a C-G basepair to close the loop. The average hydrogen bond distances were 2.87 ± 0.11 Å (U5:N3-U17:O4) and 2.96 ± 0.11 Å (U5:O₂-U17:N3) with an ideal planer geometry. The angles are 169.4 ± 6.2° and 167.9 ± 6.6°, respectively. The non-Watson-Crick U-U basepair displayed backbone torsional angles in the same range as those for other ribonucleotides in the stem, except the  torsion angle that is slightly larger than average for U17. The remainder of the stem, except for U2, was nearly canonical A-form duplex with C3'-endo sugar puckers. The rotation about C4'-C5' bond for U2 is in the ap range, extruding the U2 nucleobase, and allowing for U2 to contact the protein.

RNA/protein interaction and mutagenesis

The U2B" protein recognizes hairpin IV of U2 snRNA using an RNP motif (Uhlenbeck et al., 1997; Nagai, 1996; Burd and Dreyfuss, 1994; Mattaj, 1993). The conserved RNP1 and RNP2 sequences that occur in the two middle strands of a four-stranded anti-parallel -sheet each contribute one highly stabilizing base-stacking interaction with the RNA (Y13 in RNP2 and F56 in RNP1). The loop region of hairpin IV of U2 snRNA was involved in most of the contacts that stabilized the complex with U2B". U2B" mainly interacted with U2 snRNA in three regions. Region 1 involved 1 and helix A (Asn-15-Lys-23). Region 2 involved residues near the C-terminus of 2 and the N-terminus of 3, and those residues that link these -strands (V44-F56). Region 3 involved residues near the C-terminus of 4 and those residues that link this strand to helix C (Q85-D92). These regions were selected by comparison of the contact area of the U2B" protein and the U2 snRNA hairpin IV based on the x-ray structure and the structure derived from the MD simulations.

Mutagenesis is an important experimental method used to determine the residues important for macromolecular association and has been used to explore protein/protein (Bogan and Thorn, 1998; Jensen et al., 2000), DNA/protein (Lundquist et al., 1997; O'Neill et al., 1998), and RNA/protein (Gordon et al., 1999; Harada et al., 1998) interfaces. Computational mutagenesis is conceptually similar to experimental mutagenesis; however, it does not require expression of mutant proteins. Rather, the difference in binding free energies between the native complex and a mutant complex of identical structure, except lacking one amino acid side chain, are compared. The mutant free energy, G, is indicative of the contribution provided by the one amino acid side chain to the stability of the native complex. The mutant free energy consists of the solvation energy, G_solv, and the molecular mechanics energy in vacuo, E_gas

(1)

The solvation free energy G_solv was calculated using the general Born approach based on the continuum solvent model (Still et al., 1990), while the nonpolar contribution was calculated from the solvent-accessible surface area (Sanner et al., 1996)

(2)

The molecular mechanics energy in vacuo, E_gas, can be obtained from the internal energy E_int (including bonds, angle, and torsional terms), the electrostatic energy, E_elec, and the van der Waals energy E_vdw,

(3)

Computational glycine mutagenesis was used to identify the amino acids that provided the greatest stabilizing contributions to the RNA/protein complex. The binding free energy (difference between complex and free RNA and protein) for the wild-type and each mutant protein was calculated in terms of the components to the free energy listed in Eqs. 1-3. Here, we assume the unbound RNA and protein have the same conformation as in the complex. In the free energy calculations for the mutated protein, the side chain of the mutated amino acid was replaced by a hydrogen atom (glycine), so that the mutant exclusively reflected the interaction between the side chain of the mutated residue with all other components of the complex. Substantial values for the mutant free energy indicate the side chain plays an important role in complex stability, but does not assure that the side chain is involved in direct RNA/protein interactions. Negligible values of the mutant free energy, however, indicate that the mutated amino acid does not contribute significantly to complex stability, and is unlikely to be important in complex formation. The mutant free energy must be calculated over a substantial region of the MD trajectory during which the structure is stable. In the present case, the final 1.8 ns of the 2.2-ns MD trajectory were used for the mutant free energy analysis.

The mutant free energies for amino acids in the three regions implicated in RNA/protein complex formation are listed in Tables 1-3. Negative values correspond to unfavorable changes in the free energy of the complex as a consequence of the mutation. In region 1, the amino acids with large unfavorable mutant free energies were N15, N16, D19, K20, K22, and K23. N15 had an unfavorable mutant free energy of 8.69 kcal/mol, which derived mainly from loss of a hydrogen bond between N15, R83, and U9 of the RNA (N15:OD1-U9:O2', N15:OD1-R83:NH2). Similarly, the unfavorable mutant free energies for N16, D19, K20, and K22 resulted from loss of hydrogen bonds. Unfavorable mutant free energies in region 2 resulted mainly from lost van der Waals interactions. For example, V44 mainly interacted with A12 and U14 of the RNA by van der Waals interactions that were absent in the G44 mutant. K47 did not form a hydrogen bond with any nucleotide or residue in the native complex, and the unfavorable mutant free energy for this residue resulted from the loss of electrostatic interactions with C17 and U18 of the RNA. F56 interacted with C11 and stacked with A12 in the native complex; however, these interactions were absent in the G56 mutant. The loss of hydrogen bonds resulted in unfavorable mutant free energies for amino acids in region 3. K88 formed several hydrogen bonds with C11, while S91 formed a hydrogen bond with A12. These hydrogen bonds were one of the most important sources of the large mutant free energy.

View this table: [in this window] [in a new window]   TABLE 1   The mutant free energy (kcal/mol) by glycine scanning (1-helix A)

View this table: [in this window] [in a new window]   TABLE 2   The mutant free energy (kcal/mol) by glycine scanning (2-3)

View this table: [in this window] [in a new window]   TABLE 3   The mutant free energy (kcal/mol) by glycine scanning (4-helix C)

Many of the residues with larger mutant free energies were charged residues (e.g., K20, K22, K23, K47, K50, K88, R52), or had a polar side chain (e.g., N15, N16). These results are consistent with the unfavorable mutant free energy arising from loss of an ion pair or hydrogen bond present in the native complex. Similar results have been obtained in protein/protein complexes (Sheinerman et al., 2000). Besides the hydrogen bond and ion pair interaction, charged residues such as Lys and Arg can also form cation- interactions with the aromatic side chains of proteins or bases of nucleic acids, and loss of these interactions can also result in an unfavorable mutant free energy (Gallivan and Dougherty, 1999).

Hydration of the RNA/protein complex

Nucleotides in the loop region of hairpin IV formed a series of hydrogen bonds with residues in U2B" that stabilized the complex. Hydrogen bond strength is not readily predictable and can vary between 1 and 5 kcal/mol, depending on the heavy atom pair, internuclear distance, angle, local electrostatic environment, and solvent accessibility (Fersht, 1987; Erion et al., 2000). Direct formation of hydrogen bonds between RNA and protein are summarized in Table 4. The occupied time is the total time that the hydrogen bond existed during the 2200-ps MD simulation, and is indicative of hydrogen bond stability. Only hydrogen bonds with residence times longer than 1100 ps are included in the table. In Table 4 there are some hydrogen bonds listed that involve backbone atoms; for example, the N4 amino group of C10 formed hydrogen bonds with the backbone amide oxygen of Y86. O2 and N3 of C10 formed a bifurcated hydrogen bond with the backbone amide proton of K88, and this hydrogen bond persisted during the course of the entire MD trajectory of 2200 ps.

Formation of hydrogen bonds and salt bridges between protein residues and nucleotides in the 3' portion of the loop also showed the role played by the ion pair and polar residues in the RNA/protein complex. The number of hydrogen bonds and salt bridges formed between U2B" and this region of RNA was limited, although several hydrophobic contacts between RNA and protein stabilized this region of the complex. The salt bridge formed between the  amino group of K23 and O2P of C16 was, however, one of the RNA/protein interactions sustained for >2000 ps during the MD simulation (Table 4). The R52 guanidinium group also formed a salt bridge with U17 O1P stabilizing the complex in this region. Formation of these salt bridges reduced the occupancy of water molecules associated with the phosphodiester backbone at these positions. One other salt bridge, involving the phosphodiester at C1 in the stem region and K22, also stabilized the complex. A number of residues in the 1-helix A of B2", including K22, also form hydrogen bonds with nucleotides in the stem region of hairpin IV, stabilizing the complex.

Water at the RNA/protein interface

The stability of the RNA/protein complex is derived not only from direct interactions between protein and RNA, but also from interfacial water molecules. Water in the U2B"/U2 snRNA hairpin IV complex can be divided into three types: 1) rapidly exchanging bulk water; 2) surface water forming the first solvation shell; and 3) bridged water molecules that form hydrogen bonds simultaneously to both RNA and protein. Bridging water molecules frequently reside in the cavities of macromolecular complexes for relatively long durations. Several studies have described the occurrence of bridging water molecules at the DNA/protein interface (Schwabe, 1997). Relatively few studies on the importance of bridging water molecules in stabilization of RNA/protein complexes have been described, however, due in part to the lack of structural data for RNA/protein complexes.

Water molecules that form hydrogen bonds simultaneously with both U2 snRNA hairpin IV and U2B" in the RNA/protein complex (bridging water molecules) contribute significantly to complex stability. The criteria used for hydrogen bonding of water molecule was X ··· Y < 4.0 Å between hydrogen bond donating atoms in RNA (X), and protein (Y) and an X ··· H ··· Y angle >120°. To identify those water molecules that were truly involved in bridging hydrogen bond formation, only those water molecules that simultaneously bonded with both RNA and protein were considered as bridging water molecules. During the course of the 2200-ps MD trajectory, there were 1153 water molecules with residence times >1 ps within the inner solvation shells of both the RNA and the protein. A total of 45 bridging water molecules were identified in the trajectory for the U2B":U2 snRNA hairpin IV complex, all with residence times longer than 200 ps (Table 5). Additional bridging water molecules that underwent fast exchange were also identified at the RNA/protein interface. The bridging water molecule having the longest residence time in the trajectory formed hydrogen bonds with O2P of C16 on U2 snRNA hairpin IV (C16:O2P) and the side chain of E24 in U2B" (E24:OE2). Distance fluctuations for two of the most persistent interfaced water molecules are shown in Fig. 4. Water 624 had large fluctuations of position during the initial several hundred picoseconds of the trajectory, but was stable during the remainder of the trajectory. In contrast, water 2115 was stable nearly from the beginning of the trajectory.

View larger version (26K): [in this window] [in a new window]   FIGURE 4   The distance between the bridging water molecule and U7:O4 and D19:OD1 during the course of the 2200-ps trajectory. The structure of this bridging water molecule was remarkably stable during the trajectory. (A) The residence time is 1800-2000 ps. (B) The residence time is 2000-2200 ps.

Due to the high mobility of most waters solvating the complex, it is not too surprising to note that different water molecules can reside in the same position of the RNA/protein interface at different times during the trajectory. Conversely, a single water molecule can reside at different sites in the complex during the trajectory. For example, WAT258 formed a hydrogen bond with the side chain of D42 and the phosphate oxygen of U13 of duration 66 ps. Later in the trajectory, this same water molecule formed a hydrogen bond with the side chain of D92 and the sugar ring of G12 of 14-ps duration. Both WAT277 and WAT623 were found to hydrogen-bond to the base of A11 and the side chain of T89 at different instants during the trajectory. Thus, hydration of the RNA/protein interface is dynamic, with multiple water molecules contributing to the stability of the complex throughout the trajectory.

Defining the number of water molecules that reside at the RNA/protein interface is difficult because this number is time-dependent. A simple estimate is to use the time average over the trajectory. i.e.,

(4)

where t_i is the residence time of each bridging water and T is the trajectory time, which is 2200 ps in the present simulation. Using this approach, the average number of bridging water molecules residing at the RNA/protein interface at any instant is 20 (considering only the water with residence time >200 ps). Because in this estimate only water molecules that resided for at least 200 ps were considered, the estimate is the lower limit of the number of bridging water molecules that occurred during this simulation. The average value of the maximum observed residence time of water molecules (Qian et al., 1993; Sunnerhagen et al., 1998) residing longer than 200 ps during the MD simulation was 979 ps.

It is also interesting to note that most residues with a large mutant free energy can hydrogen-bond with RNA, either directly or indirectly, through bridged water, or both (Tables 4 and 5). For example, OD1 of N15 hydrogen-bond with O₂' of U8 on RNA for 1576 ps, and also indirectly interact with O2P of G9 through bridging water WAT895 for 592 ps. NZ of K20 does not form hydrogen bonds with RNA that persist for long times, but three bridging water molecules (WAT2416, WAT3092, and WAT3522) do hydrogen-bond with O4 of U5, O1P of U2, and O4 of U5. Mutation at this residue would cause a reduction in the number of hydrogen bonds between RNA/protein and protein/bridging water, as well as the direct interaction.

Interaction energy associated with water

Based on the above analysis of the MD trajectory, water molecules reside on the interface of the RNA/protein complex. However, it is difficult to calculate the free energy associated with the bridging water molecules directly from the MD calculation in this case. In this paper, the interaction energy associated with those residual water molecules was estimated using ab initio calculations for a model system. Water molecule 624, associated with Cyt-10 and Thr-89, as identified from the analysis above, were used as a model system for the ab initio calculations (Table 6). The geometry of the starting complex was taken from the minimized structure for the complex averaged over the residence time during the trajectory. The electronic structure was fully optimized at the HF/6-31G* level. The optimized geometries of the Cyt-10 and Thr-89 were taken directly for the single-point calculations.

The interaction energy E can be understood as the sum of pairwise dimer contributions and a three-body term, E³, which accounts for cooperative effects (Sponer et al., 1997; Brandl et al., 1999, 2000) in the water 624/C10/T89 system.

(5)

The direct interaction between C10 and T89 without water is straightforward with consideration of the BSSE correction. The interaction energy was further corrected by consideration of the deformation energy and the zero-point energy. The deformation energy in the C10/T89 dimer is 0.7 kcal/mol for both HF and DFT calculations. The zero-point energy correction is 0.4 kcal/mol for HF and 0.3 kcal/mol for the DFT calculation. The calculated interaction energy of 1.3 kcal/mol (HF) and 1.6 kcal/mol (DFT) is much smaller than the value in the water-mediated trimer (Table 6). The total deformation energy from each component in the trimer is 0.8 kcal/mol (HF) and 1.0 kcal/mol (DFT). The zero-point energy is increased from 0.4 kcal/mol (dimer) to 4.1 kcal/mol (trimer) in the HF level calculation, and from 0.3 kcal/mol (dimer) to 2.6 kcal/mol (trimer) in the DFT level calculation. The E³ value of 0.3 kcal/mol (HF) and 2.3 kcal/mol (DFT) indicated that DFT indicated a stronger cooperativity than HF calculation in the water-mediated trimer. The difference of the 7.2 kcal/mol (HF) and the 10.5 kcal/mol (DFT) between the dimer and trimer is the stabilization energy associated with the water-mediated trimer compared with the direct interaction in the dimer (Table 6). The stabilization energy is derived from the additional hydrogen bonds due to the bridged water and orbital interactions among those three components in the trimer. The results of the ab initio calculation indicated that water 624 greatly stabilized the direct interaction between C10 and T89. To explore the overall effects from all bridged waters at the protein/RNA interface, it would be necessary to do a statistical average for the total stabilizing energy including all bridging waters. However, considering the large number of the bridged water in the MD trajectory, it would be time-consuming to complete calculations at the ab initio level for this system. Therefore, the ab initio calculations reported here do not provide a definite value for the entire stabilization energy from all of the bridged waters in the RNA/protein complex, rather they only provide a possible range for the stabilization interaction energy for each water-associated trimer at the RNA/protein interface. We expect future studies to provide additional understanding of this kind of interaction.

Comparison with U1A'/U1 snRNA hairpin II

The molecular dynamics of the U1A'/U1 snRNA hairpin II complex (Reyes and Kollman, 1999) were investigated by Reyes and Kollman using both alanine scanning mutagenesis (Reyes and Kollman, 2000a), and binding thermodynamics (Reyes and Kollman, 2000b). The computational results from these studies were in excellent agreement with the available experimental results. The loop structure of U1 snRNA hairpin II is very similar to U2 snRNA hairpin IV, the subject of the present investigation. The major differences between these stem loops are that C12 in U1 snRNA hairpin II was replaced by the G12 in U2 snRNA hairpin IV, and there is one more nucleotide (A14) in U2 snRNA hairpin IV. Although U2 snRNA hairpin IV binds U2B" in an orientation similar to the way U1 snRNA binds to the U1A' protein, the details of the interaction differ in several ways (Price et al., 1998). In the U1A'/RNA complex, the most important residues for binding were in the 2-3 loop and RNP regions, which is similar for U2B"/U2 snRNA hairpin IV. The other two interaction areas (1-helix A and 4-helix C) in U2B" were not, however, observed in the U1A' protein complex.

MD simulations of the U1/U1A' complex RNA/protein complex also analyzed hydration of the complex and the role of bridging water molecules in complex stability (Tang and Nilsson, 1999). In this RNA/protein complex, which is highly homologous to the U2B"/U2 snRNA hairpin IV complex considered in the present study, the lifetime of the bridging water molecules identified was typically <20.0 ps (Tang and Nilsson, 1999). It is surprising to note such a large difference in the stabilities of bridging waters between these two complexes. One possible reason for the shorter residence time for bridging water molecules in the previous study compared to the present results may derive from differences in computational methods. In the present study, PME has been used to properly account for long-range electrostatic interactions in the RNA/protein complex, especially hydrogen bond formation. Appropriate accounting of long-range electrostatic interactions is very important for deriving meaningful results for MD simulations involving nucleic acids, which should routinely be considered using the PME sum (Essman et al., 1995). Other studies also observed residence time as long as 5.1 ns in zipper-like DNA (Spackova et al., 2000).

	CONCLUSIONS

TOP ABSTRACT INTRODUCTION METHODS RESULTS AND DISCUSSION CONCLUSIONS REFERENCES

In summary, a 2200-ps MD simulation of the solvated U2 snRNA hairpin IV/U2B" protein complex was completed and the trajectory has been analyzed. The MD simulation in water provided several insights into the factors that provided important contributions to complex stability. Analysis of the mutant free energies indicated the most important regions contributing to the stability of the RNA/protein complex. Charged and polar residues played an important role due to favorable electrostatic interactions. For example, several lysine residues showed a large mutant free energy in the RNA/protein interface that arose due to interaction with the phosphodiester backbone of the RNA. The formation of direct hydrogen bonds and salt bridges between the RNA and protein contributed to complex stability. The overall interaction between U2 snRNA hairpin IV/U2B protein differed from the U1 snRNA/U1A RNA/protein complex, although the two complexes had similar secondary structures; 1153 waters were observed in bridging between RNA and protein in this 2200-ps trajectory. Only 45 of them had a residence time >200 ps (~9% of the trajectory). Other bridge waters exchanged rapidly at the interface. On average, 20 bridging water molecules occurred at the interface of the complex, and the average residence time of those bridged waters was 979 ps. Calculation for the interaction energy resulting from bridging water molecules was determined for the C10/T89 interaction using ab initio calculations. These results indicate that bridged water play an important role in stabilization of the RNA/protein interface.