|Home | About | Journals | Submit | Contact Us | Français|
We develop a partial charge-based Tightly Bound Ion (PCTBI) model for the ion effects in RNA folding. Based on the Monte Carlo Tightly Bound Ion (MCTBI) approach, the model can account for ion fluctuation and correlation effects, and can predict the ion distribution around the RNA. Furthermore, unlike the previous coarse-grained RNA charge models, where negative charges are placed on the phosphates only, the current new model considers the detailed all-atom partial charge distribution on the RNA. Thus, the model not only keeps the advantage of the MCTBI model but also has the potential to provide important detailed information unattainable by the previous MCTBI models. For example, the model predicts the reduction in ion binding upon protein binding and ion-induced conformational switches. For Hepatitis C virus (HCV) genomic RNA, the model predicts a Mg2+-induced stabilization of a kissing motif for a cis-acting regulatory element in the genomic RNA. Extensive theory-experiment comparisons support the reliability of the theoretical predictions. Therefore, the model may serve as a robust starting point for further development of an accurate method for ion effects in RNA conformational equilibrium and RNA-cofactor interactions.
Metal ions play an important role in the stabilization of three-dimensional (3D) structures of nucleic acids (DNAs and RNAs), due to their ability to neutralize/screen the negative charges on nucleic acids backbones and specific ion-nucleic acid interactions.1, 2 Quantitative prediction for ion binding properties is essential for the understanding of nucleic acids structural formation and folding stabilities.3–20 However, despite tremendous theoretical and experimental efforts, accurate prediction of the ion effects at various solution conditions for different RNA structures remains one of main challenges in computational modeling of RNA structure and stability,21–25 especially for solutions involving multivalent ions such as Mg2+.
According to the different charges and sizes and the different ways of binding to the RNA, we can distinguish three types of ions.26–28 First, in one extreme, ions can be dehydrated and trapped at specific sites in RNA29 and ion-RNA binding is stabilized by specific interactions such as chelation with the electronegtive atomic groups of RNA. These ions are site-bound (or chelated). Second, in the other extreme, the majority of the bound ions can remain hydrated and bind to RNA diffusely. The large number of diffusely bound ions may dominate the ion electrostatic energy of the system.30 These ions “freely” diffuse at a distance from the RNA surface without any specific binding sites. These ions are called diffuse ions or diffusely bound (DB) ions. Third, there exist ions which are clustered around the RNA surface due to ion-RNA Coulomb attraction, resulting in high local concentration of ion. These (hydrated) ions are neither trapped at specific binding sites like chelated ions nor “free” to move around like diffuse ions. The high local concentration of the ion can cause strong coupling/correlation between the ions due to volume exclusion and Coulomb interaction. As a result, the local electric field for an ion is not only related to the position of the ion itself but also to the positions of the nearby ions. Such a correlation effect is more pronounced for multivalent ions such as Mg2+ due to the high ionic charge than for monovalent ions.13, 31–33 We call such ions, which involve the strong correlation, as the tightly bound (TB) ions.
A number of different methods have been used to investigate the ion effects in RNA folding and stability. Among these methods, molecular dynamics simulation34–37 has the advantage to provide insights into atomic details of the ion binding process, such as the site-specific ion binding in RNA kissing loop formation,38 the 5S rRNA loop E motif,39 SAM-I riboswitch,40 and many other systems.41–44 However, the accuracy of MD simulations is limited by the accuracy of the ion-dependent force field as well as the completeness of the sampling for the ions and water molecules,34 As a result, accurate calculations require significant computer resources and exceedingly long computer time even for small RNAs. The classical Counterion Condensation (CC) theory45 and the (nonlinear) Poisson-Boltzmann (NLPB) theory46–51 are simple yet useful theories for the prediction of the ion effects. However, both theories ignore ion correlation and ion fluctuations, two potentially important effects for ions around nucleic acids.31, 52 The CC theory treats ion binding with average (partially) neutralized charges while NLPB considers ions with an average continuous distribution. Recently, several new theories have been developed to account for the ion correlation effect for biomolecular systems.53–55 For example, the three-dimensional reference interaction site model (3D-RISM) was developed based on the Ornstein and Zernike integral equation theory,54 and a generalized counterion condensation theory was proposed to explicitly treat Mg2+ ions with ion-ion correlations.55
The tightly bound ion (TBI) model is another highly promising model that can accurately and efficiently predict the ion binding effects for RNAs.28, 52, 56–61 The key idea of the TBI model is to explicitly enumerate the many-body, correlated ion distributions for the TB ions while using a mean-field description (NLPB) for the DB ions. In this way, the model takes into account both the correlation and the fluctuation effects. Extensive comparisons between TBI predictions and experimental results support the validity and accuracy of the TBI model for predicting ion binding properties and ion-mediated nucleic acids stability for simple helices, pseudoknots, kissing complexes, and other more complex tertiary folds.28, 52, 56–61 However, because the original model is based on the explicit enumeration of ion distributions, which is time-consuming, the model is computationally inefficient. The low computational efficiency renders applications to medium (100–200 nts) and large RNA structures (> 200 nts) impractical.
Recently, we developed a new version of the TBI model, Monte Carlo tightly bound ion (MCTBI) model,62 which employs Monte Carlo Insertion Deletion (MCID) algorithm to sample the ion distributions of TB ions. In comparison with the original TBI model, the MCTBI model has two main advantages: more efficient computational speed and a higher resolution (3D coordinate-based) description for the ion distribution. The MCTBI model enables prediction of the highest probability regions for ion binding. However, similar to the original TBI model, the MCTBI model is based on a coarse-grained (CG) charge model for RNA, where the charges on a nucleotide are assumed to “collapse” into a point charge located at the center of the phosphate. In such a CG charge model, each phosphate carries an electronic charge −e and the rest atoms are assumed to be neutral.52, 60, 61 The coarse-grained charge distribution has the advantage to reduce the computational time. However, the over-simplification of the charge distribution prevents the model from applications to many biologically important problems which require information about realistic and detailed charge distributions on RNA. These problems include predictions of ions binding sites at higher accuracy and resolution63 and modeling of ligand-RNA and protein-RNA docking.64 These and many other ion-dependent RNA biology problems demand a new model that goes beyond the CG charge-based theory.
In this paper, we report a new model, PCTBI model, which uses all-atom partial charges for RNA atoms and can treat ion correlation and fluctuation effects. Tests against the previous CG charge-based model using the experimental data supports the validity and accuracy of this new model. Furthermore, compared with the CG charge model, the partial charge-based model reveals important details about the ion effects, including the influence of cofactor (protein) binding on ion distribution around the RNA. As an application of the model, we investigate the ion-dependent conformational equilibrium and conformational switches in Hepatitis C virus (HCV) genomic RNA.
Belonging to the Flaviviridae family, HCV is a Hepacivirus that causes the most common chronic blood-borne viral infection in the United States.65 Replication of HCV in human liver cells triggers and invades host immune responses, which lead to extensive damage to liver.66
Cell-culture based viral replication assays suggested a potential long-range RNA:RNA interaction in the 3′-end of the genomic RNA between 3′X and an upstream cis-acting regulatory element termed 5BSL3.2.67 The base pair complementarity between 5BSL3.2 and 3′X are phylogenetically conserved, and has been shown to be essential for HCV viral protein translation and RNA synthesis.68–70 Previous gel electrophoresis and surface plasmon resonance (SPR) studies have shown that in vitro transcribed 5BSL3.2 and 3′X can form an RNA:RNA complex and Mg2+ is required in the gel and running buffer in order to detect the RNA complex band.71, 72 We apply the PCTBI model developed in the present study to investigate the physical mechanism for the Mg2+-dependence of the folding stability and the shift of conformational equilibrium for the 5BSL3.2:3′X kissing complex. Our results indicate that Mg2+ ions can indeed stabilize the 5BSL3.2:3′X docked state over the undocked state and for 3′X, adding Mg2+ to the solution shifts the conformational equilibrium toward the kissing conformation.
The development and test of the PCTBI model involve a number of RNA structures. The 80-nucleotide (nt) A-form RNA duplex and 48-nt B-form DNA helix are generated using the latest version (version 2.1) of the program X3DNA.73 Other RNA structures, such as the 76-nt yeast tRNAPhe (PDB ID: 1TRA74), 58-nt fragment of rRNA and rRNA-protein complex (PDB ID: 1HC875), and 72-nt Adenine riboswitch (A-riboswitch; PDB ID: 1Y2676), are generated from the Protein Data Bank coordinate files.77 For the 5BSL3.2:3X complex in the HCV genomic RNA, cell-based viral replication assays68–70 have suggested that the kissing complex is formed between 5BSL3.2 and the SL2 hairpin loop in the 3′X.67 We therefore focus on the two RNA segments 5BSL3.2 and SL2. The conformations are generated by the Vfold model combined with molecular dynamics (MD) simulations.78 In the MD simulations, we use the structure predicted by the Vfold model78 as the initial structure, and then apply a harmonic potential U(r) = k(r − r0)2 to constrain the base-pairs, where k = 1.0 kcal/mol and r0 = 3.5 Å. For each 2D structure (defined by the set of base pairs contained in the structure), totally more than 5000 different 3D conformations are generated in the MD simulation.
The partial charges of the atoms in RNA are assigned using the “Dock Prep” module in Chimera.79 Before implementing the “Dock Prep” module, all the ions in the PDB are removed at first. Then, for a given structure, we run the “Dock Prep” module runs multiple times for solvent deletion, alternate locations deletion (keeping the highest occupancy), hydrogen addition, partial charges addition, and output with Mol2 format. Charges for standard residues in RNAs are taken from AMBER ff14SB80, 81 and charges for nonstandard residues, including ligand for the A-riboswitch, are calculated using the ANTECHAMBER module with the AM1-BCC charges.82, 83 In order to test the sensitivity of the different partial charge assignment methods, we also use PDB2PQR method84 to calculate the partial charges for the 58-nt fragment of rRNA with AMBER99,85 CHARMM22,86 or PARSE87 force fields.
For a given RNA with all-atom structure and partial charges assignment, we first run NLPB calculation to roughly compute the ion concentrations for each ion species around the RNA. In the NLPB computation, RNA is located in a large solution box, whose size is larger than six times of the Debye length in order to reduce the boundary effect.28 The mixed salt solution in the box can contain divalent cations Mg2+, monovalent cations K+ (or Na+), and monovalent anions Cl−. Their bulk concentrations , and satisfy the charge neutrality condition: .
To accelerate the NLPB calculation, we apply the three-step focusing process,88 where grid sizes 1.7 Å, 0.85 Å, and 0.425 Å for the three steps, respectively. From the ion concentration, we determine the average ion-ion distance and the ion correlation strength, which is defined as the ratio between the average ion-ion Coulomb energy and the thermal energy.28 Because the ion concentration is location-dependent, the ion correlation strength is coordinate-dependent. According to the spatial distribution of the correlation strength, we classify the region around RNA into TB region and DB region, for the TB ions (of strong correlation) and the DB ions (of weak correlation), respectively. For monovalent ions, which involve weaker Coulomb interactions than multivalent ions, the TB region (of strong correlation) is negligible and therefore, all the monovalent ions can be treated as DB ions.
Unlike monovalent ions, multivalent ions (such as Mg2+ ions), usually involve a significant TB region, which forms a thin layer of several Ås plus the ion radius from the RNA surface. Here we assume that ions are hydrated with radii rMg2+ = 4.5 Å, rK+ = 4.0 Å, rNa+ = 3.5 Å, and rCl− = 4.0 Å, respectively.58, 59 According to the above analysis, for a monovalent/divalent ion mixture solution, only for divalent ions (e.g., Mg2+) we need to distinguish TB ions from DB ions. The DB ions can be treated with the mean-field NLPB theory while the (strongly correlated) TB ions require sampling of many-body, discrete ion distributions.
For an RNA with N phosphates, we allow Nb, the number of the tightly bound Mg2+ ions (in the TB region), to vary from 0 to N. The total partition function of the many-body system can be calculated as the sum over Nb:
where Z(Nb) is the partition function for a given Nb (the number of Mg2+ ions in TB region):
Here Zid represents the ideal partition function without the RNA in the system (uniform ion solution). is the bulk concentration of Mg2+ ions and the factor Nb! accounts for the nondistinguishability of the Mg2+ ions. kB and T denote the Boltzmann constant and temperature. The statistical weight W (Nb) is equal to the configurational integral for all the Nb ions:
Here Ri denotes the coordinate of ion i in the TB region. ΔGb in Eq. 3 and ΔGd in Eq. 2 represent the interaction energy of all the charged particles in TB region (TB ions and atoms of RNA) and the free energy of DB ions, respectively. ΔGb includes volume exclusion, Coulombic interaction and dielectric polarization energies. In our calculation, the excluded volume effect is modeled by a Lennard-Jones (LJ) potential and the polarization energy is computed from the Generalized Born (GB) model:
The first term in the equation above is the Coulombic energy for the interaction between charges in the TB region. Zi(or j)e is the charge of particle i(or j), εR (= 20 in our calculation) is the dielectric constant of RNA, and rij is the distance between particles i and j. The second term is the LJ potential with uo (= 0.35) as the LJ constant and σij as the equilibrium distance between particles i and j. Here we set σij as the addition of the radius of the two particles. The third term above is the mutual polarization energy induced by other charges, εW (= 78) is the dielectric constant of water and Bi (or j) is the Born radius for particle i(or j). The forth and fifth terms represent the self-polarization energies of RNA (subscript R) and ions (subscript I). The free energy ΔGd of DB ions can be calculated from the effective single-particle ion distribution solved from NLPB:89, 90
Here α denotes the ion species. zαe is the charge of ion species α, ψ(r) and ψ′(r) are the electrical potentials at position r with and without ions in the solution, and cα(r) and represent the local (at r) and bulk concentrations, respectively.
We use a simple cubic lattice in TB region to configure the ion distributions for the TB ions. Assuming that the TB region contains Ns lattice sites with lattice grid length lb. Then the statistical weight in Eq. 3 can be computed from the following formula
where mi represents the number of all the available (empty) sites for TB ions. ΔGb is the interaction energy of all the charged particles in the TB region (see Eq. 4). The lattice size lb is adaptive in order to reach an optimal balance between the resolution of ion distribution and the efficiency of computation. For a low bulk concentration of Mg2+ ions, lb = 0.425 Å is determined from the grid size of the third focusing step in the NLPB calculation (for the DB ions). For a high ion bulk concentration, there is a larger TB region (sampling space) so that the sampling of the ion distribution could be time-consuming. Therefore, we increase the value of lb such that the number of lattice sites Ns inside the TB region is ≈ 200 × N. To circumvent the problem of the exceedingly long computational time for the exact enumeration of the ion distributions, we use the Rosenbluth-Rosenbluth (RR)-based Monte Carlo method91 and the Monte Carlo Insertion Deletion (MCID) algorithm,62 to sample the ion distributions in the TB region. The central idea of MCID algorithm is to sample ion distribution by inserting particles one by one according to the energy followed by a re-sampling for the ions by randomly removing high-energy ions. The purpose of the resampling/relaxation is to further enhance the sampling quality for the important (low-energy) ion distributions, we allow the multiple-ion system to relax.92 Specifically, we perform a re-sampling for the ions by randomly removing high-energy ions. The details of MCID algorithm has described in the Ref. 62.
HCV-3′X and 5BSL3.2 stock solutions were created at a concentration of 25 μM in 10 mM Tris-HCL (pH 7.48). These starting stocks were boiled for 3 minutes, and then snap cooled at 4 °C. The 3′X and 5BSL3.2 samples were then mixed at 1:1 ratio in the presence of 1x incubation buffer (final concentration: 140 mM KCl, 10 mM NaCl, 2 mM MgCl2, 10 mM Tris-HCl, pH 7.5), and incubated at 37 °C for 60 minutes. The mixture was analyzed by electrophoresis in acrylamide gels containing 0, 1 and 2 mM MgCl2. The gels were run at 4 °C, and the RNA bands were stained by stains-all and imaged using a Kodak camera.
Purified 5BSL3.2 RNA samples were dissolved in 90% H2O, 10% D2O with 20 mM Tris-HCl, pH 6.5 at different MgCl2 concentrations. NMR spectra were acquired on a Bruker Avance III 800 MHz spectrometer equipped with TCI cryoprobe at the NMR core facility, University of Missouri, Columbia. For the assignment of imino protons, two-dimensional NOESY spectra were recorded for 5BSL3.2 with 0 and 5 mM MgCl2 at 15 °C. Data were processed with NMRPipe and NMR Draw,93 and analyzed with NMRViewJ (One Moon Scientific, NJ).
To test the accuracy of the new partial charge-based PCTBI model, we calculate the fraction of excess bound ions fα per nucleotide for ion species α, and compare the theoretical predictions with the experimental data. fα is given by the following formula:
In experiments, the excess number of bound ions can be measured by “Ion Counting” methods,94 such as fluorescence titration95–99 and buffer equilibration and atomic emission spectroscopy (BE-AES).31, 100 The fluorescence titration measurements take advantage of a fluorescent dye for an ion species such as Mg2+ to count the bound ions. Therefore, the method focuses on single ion species in each measurement. The BE-AES method that uses a filter to create a sample solution (with nucleic acids) and a reference solution (without nucleic acids) allows counting of the excess bound ions in the solution by atomic emission spectroscopy.31 The method has a lower detection limit100 and can also measure anion depletion. In the following calculations, unless otherwise stated, the partial charges of RNA atoms are calculated using the “Dock Prep” module in Chimera79 with AM1-BCC (or AM1 for short) charges.82, 83
Fig. 1 shows the binding fractions fMg2+ and fNa+ (or K+) (the average number of bound ions per nucleotide) of divalent and monovalent ions for the different RNAs as predicted by the original CG charge-based (MCTBI) model, where each phosphate group carries −e charge and other atoms are neutral, and the partial charge-based (PCTBI) model. The results for the binding fractions lead to three conclusions. First, our theoretical predictions are in accordance with the experimental data. Second, the PCTBI model gives slightly better predictions than the previous MCTBI model, such as the aRNA case at [Mg2+] = 0.1 mM. In order to make direct comparisons between the two models, we also show the prediction data in Table 1 for a 24-bp B-DNA, for which the experimental data is available.100 Indeed, the relative error ΔE (see table 1) shows that the PCTBI is (slightly) better than MCTBI. Third, fK+ of RNA-protein complex predicted by the PCTBI model is lower than that predicted by the original MCTBI model [see Fig. 1F]. The computational time increases with the RNA length (number of nucleotides) and the salt concentration. For example, for the 58-nt rRNA in a 20 mM K+ background solution shown in Fig. 1C, the computer times are about 12 and 30 minutes for 10−6 M and 10−2 M Mg2+, respectively. The computations were performed on a PC with an Intel i7-4790 processor and 16 GB RAM. On the same computer, for the 80-nt RNA shown in Fig. 1A, the computer times are about 30 and 60 minutes for 10−6 M and 10−2 M Mg2+, respectively.
The original MCTBI model, unlike the current PCTBI model, cannot treat the partial charges on the protein, which can be critical for the RNA-protein docking. Specifically, the MCTBI model neglects the protein-induced change of Coulomb interactions as the residues of protein are assumed to be electrically neutral. In contrast, the PCTBI model can account for the detailed partial charge distributions not only on the RNA but also on the protein. Of particular importance is the charge distribution on the RNA-protein interface. The protein surface on the solvent side has positive partial charges, which tend to reduce cation binding. Therefore, the coarse-grained charge model, which cannot reflect such detailed charge distribution on the protein, would over-estimate the ion binding (see Fig. 1F).
To test the sensitivity of the predictions to the different partial charge models, we compute the ion binding fractions for the 58-nt rRNA fragment based on the different partial charge assignments, including the default method (“Dock Prep” with AM1),80–83 pdb2pqr84 with AMBER99 (AMB),85 CHARMM22 (CHA),86 and PARSE (PAR)87 force fields. Although the net charges on the 58-nt rRNA predicted by the above four methods are the same (−56e), the charge distributions are different. Using the partial charge assignment for a phosphate group in a standard nucleotide as an example (see Table 2), AM1 and AMB have the same charge assignment since they both use Parm9981 as the force field for a standard nucleotide. The charge assignments by PAR and CHA are different from the above two methods and the difference for CHA is more significant as shown by the net charge of the phosphate group (Table 2). As shown in Fig. 2, the partial charge models of AM1 (black line), AMB (red line), and PAR (green line) give nearly the same predictions, which are better than the predictions by CHA (blue line). Here, we note that for the case of the 58-nt rRNA with background salt of 20 mM K+, Mg2+ may be dehydrated and trapped at a specific binding site at high [Mg2+].75 Such a dehydration effect, however, is neglected in the current model. This lack of ion dehydration effect may contribute to the theory-experiment difference for [Mg2+] > 0.1 mM.
One of the important applications of the PCTBI model is to predict the electrostatic free energies as a function of the ion concentrations and the ion-induced shift of the populational shift of the different conformations and corresponding conformational switches. As shown in Fig. 3 for six RNA systems, increasing the Mg2+ concentration results in a decreasing in the electrostatic free energy ΔG, which would contribute to the stabilization of the RNA structure. As shown in the result for the 58-nt rRNA fragment with 20 (blue line) and 60 mM K+ (cyan line), high concentration of monovalent cation can also help stabilize the RNA. However, at high concentration of divalent cation, the stabilization effect from the monovalent cation becomes negligible. Furthermore, the free energies of the rRNA (cyan line) and the rRNA-protein complex (magenta line) with 60 mM K+ indicate that the rRNA-protein complex is favored by the electrostatic interactions. It is important to note that such effects cannot be predicted by the original MCTBI model.
With the use of the detailed partial charge distribution for RNA (and protein) and the sampling of ion distributions based on the 3D coordinates of the ions, the current new model enables more reliable predictions for the probable bound ion distributions and ion binding sites. The probability of finding a bound ion at position (grid point) k is given by
where n(Nb, k) is the number of Nb-ion distributions (out of the totally Nsample sampled distributions) with site k occupied by a TB ion.
For the tRNA (PDBID: 1TRA74) with 32 mM Na+ and 2 mM Mg2+ and the rRNA-protein complex (PDBID: 1HC875) with 60 mM K+ and 2 mM Mg2+ as examples, Fig. 4 shows the TB ion distributions of probability larger than 5% as predicted by the PCTBI model [Figs. 4A and C] and the original coarse-grained charge-based TBI model [see Figs. 4B and D]. The 3D structures of RNA used here were determined by X-ray diffraction, and the specific bound Mg2+ ions (cyan balls in Fig. 4) were determined from the electron density maps.74, 75 Comparison between the predicted high-probability ion distributions and the experimentally observed ion binding sites leads to two conclusions. First, the experimentally observed sites are in close proximity to our predicted “high probability” binding sites. We note that the predicted ion binding regions are usually close to the experimentally observed ion binding sites, however, the predicted binding sites may not agree exactly with the those shown in the crystallographic RNA structures. This is due to two reasons. First, the (MCTBI and PCTBI) models assume fully hydrated ions while experimentally observed specific site-bound ions often involve dehydration.29 Therefore, the (fully hydrated) ions in the models cannot enter the solvent-inaccessible pockets in RNAs. To best estimate the likelihood of ion binding, we use 5% as the cut-off probability for the probable ion binding sites. Second, ions bound to an RNA structure can be dynamic. Therefore, it is possible that the bound ions can distribute in a broad region instead of become fixed at a specific site shown in the crystallographic structure. Nevertheless, the results suggest that the theory may offer a low-resolution prediction for ion binding sites, which can be used as starting points for further refinements by including the ion dehydration effect. Second, the “high probability” locations predicted by the PCTBI model are more focused and accurate than those predicted by the original coarse grained charge-based TBI model (see the locations labeled by black and red cycles in Fig. 4). Because a Mg2+ ion has the same number of electrons as water molecules, it is possible that certain ion binding sites are not shown in the experimental structures.63, 76, 101, 102 Therefore, it might be possible that certain predicted highly probable binding sites are not explicitly shown in the crystallographic structures. Furthermore, as reported in a previous study,40 the predicted ion binding regions may also correspond to regions where ions have low diffusion.
In this section, by combining the experimental data and the theoretical predictions, we investigate the Mg2+ ions effect on the formation of the 5BSL3.2:3′X complex in HCV genomic RNA. Because the SL2 kissing loop is the main part in 3′X to form the kissing complex with 5BSL3.2, we focus on the 5BSL3.2 and the SL2 domains in our calculation.
In the experiment, the 5BSL3.2:3X complex was loaded to native polyacrylamide gels containing 0, 1, and 2 mM Mg2+. In agreement with the previous findings,71, 72 the 5BSL3.2:3′X complex band can be detected only in gels containing Mg2+ [see Fig. 5]. When Mg2+ was deprived in electrophoresis, the complex dissociated into free 5BSL3.2 and 3′X in the no Mg2+ gel, suggesting that Mg2+ is needed to maintain the complex. It is worth noting that the free 5BSL3.2 RNA migrated much faster in the Mg2+-free gel as its position was significantly lower than free 3′X, whereas in the 1 mM Mg2+ gel it was much closer to 3′X, and in the 2 mM Mg2+ gel it migrated to almost the same position as 3′X. These data suggest that Mg2+ is necessary to prevent dissociation of the 5BSL3.2:3′X complex, and Mg2+ might be involved in the folding of 5BSL3.2 due to the fact that 5BSL3.2 exhibited a pronounced decrease in migration rate in gels containing higher concentrations of Mg2+, whereas the migration rate of 3′X was relatively insensitive.
To confirm the sensitivity of 5BSL3.2 folding to Mg2+, one-dimensional imino proton spectra were collected for 5BSL3.2 with different Mg2+ concentrations. New signals were detected as the Mg2+ concentration increased, whereas the signals detected in the absence of Mg2+ remained unchanged. Assignment of these exchangeable imino protons on the two-dimensional NOESY spectra confirmed that formation of the kissing-loop conformation of the 5BSL3.2 relies on Mg2+. When the Mg2+ concentration was less than 1 mM, the imino proton signals observed were assigned to residues G9268–G9273 and G9310–G9313, demonstrating the formation of the bottom stem of 5BSL3.2. No signals indicative of a well-formed structure on the top stem were detected. The new imino signals detected at higher Mg2+ concentrations were assigned to residues of the 5BSL3.2 upper stem, demonstrating that the kissing conformation of 5BSL3.2 is stabilized by Mg2+ while the bottom stem remained unchanged [see Fig. 6a]. It is conceivable that the 5BSL3.2 upper stem residues are rapidly exchanging among various conformations in the absence of Mg2+ and cannot give rise to detectable imino signals, whereas the kissing conformation is favored in the presence of Mg2+ and thus afford new signals.
To quantitatively analyze the Mg2+ effects on the folding of 5BSL3.2, we use the Vfold method78 to predict 2D structures of 5BSL3.2 from the sequence. The 5BSL3.2 has two alternative 2D structures (shown in the insets of 6b): a stem-loop structure (kissing conformation), which can form the 5BSL3.2:3′X complex, and a three-way junction (non-kissing conformation), which cannot form the 5BSL3.2:3′X complex. For each 2D structure, we the use the Vfold3D method combined with molecular dynamics simulations78 to generate ensemble of 3D conformations for 5BSL3.2. We then apply the new partial charge-based model developed here to evaluate the electrostatics free energy for each conformation. For a given 2D structure, we predict its electrostatic free energy as the Boltzmann-weighted average electrostatic free energy over all the conformations. The predicted free energies for the two alternative 2D structures show that an increasing [Mg2+] would cause the free energy of the kissing conformation (of 5BSL3.2) to become lower than the non-kissing conformation [see Fig. 6b], indicating that a high [Mg2+] would favor the kissing structure over the non-kissing structure. Fig. 7 shows the “high probability” (> 5%) locations of Mg2+ ions for the kissing conformation (a representative 3D structure) of 5BSL3.2 with [K+] = 150 mM and [Mg2+] = 5 mM. The bulge pocket region in the kissing conformation (yellow) of 5BSL3.2 involves significant charge build-up of the RNA, resulting in a high probability for Mg2+ binding around the pocket region, especially in a solution of high [Mg2+]. Therefore, high [Mg2+] favors the kissing conformation of 5BSL3.2. As a caveat, we note that the stability of the 5BSL3.2 is determined by multiple factors, such as canonical and non-canonical hydrogen bonding, base-staking, and even the interactions with 3′X domain. Therefore, our results here only show the Mg2+-induced shift of the relative stability between given RNA structures.
In contrast to the sensitive Mg2+-dependence of the 5BSL3.2 folding, consistent with the gel electrophoresis results, no significant NMR signal changes were observed in the absence and presence of Mg2+, suggesting a weaker [Mg2+]-dependence of the 3′X folding. The quantitative effect of Mg2+ ions on the docking part of 3′X (SL2) is computed. Our Vfold2D model predicts that SL2 folds into a single 2D structure and the Vfold3D model gives the 3D structure of SL2. SL2 and the 5BSL3.2 kissing (stem-loop) conformation have different 3D structures and thus different [Mg2+]-dependence. The 5BSL3.2 kissing conformation involves a pocket structure in the bulge region (from 9298A to 9304G) and thus significant RNA backbone (local) charge build-up, which induces strong ion-binding. Therefore, 5BSL3.2 stability is highly sensitive to [Mg2+]. In contrast, the SL2 structure lacks such a ion-binding pocket and shows less [Mg2+] sensitivity. The Boltzmann-weighted average electrostatic free energy predicted by the PCTBI model for the Vfold3D-generated 3D conformational ensemble indeed shows a much weaker [Mg2+]-dependence (see Fig 6c). The result indicates that the predicted SL2 structure is stable over a wide range of [Mg2+] and the stability is not sensitive to Mg2+. The conclusion is in agreement with the experimental finding.
To understand the [Mg2+]-dependence of the 5BSL3.2-SL2 docking, we generate a docked 5BSL3.2:SL2 complex state [see the inset of Fig 6d] and an undocked state with separated 5BSL3.2 and SL2 structures. We find that the electrostatics free energy of the docked state is much lower than the undocked state, and the free energy gap increases with an increasing [Mg2+] (see 6d). Compared with the structure of undocked state, the docked state has a more compact structure and, meanwhile, create some extra TB regions for Mg2+ ions between the 5BSL3.2 and SL2. The Mg2+ ions in these TB region would subject to the attraction from the 5BSL3.2 as well as SL2. Therefore, the docked state could attract more TB Mg2+ ions, and these Mg2+ ions can indeed help stabilizing the 5BSL3.2-SL2 docked complex (over the undocked state).
In this paper, we describe a new partial charge-based model, PCTBI, to predict ion effects for RNA. Although both PCTBI and MCTBI models take into account the ion correlation and fluctuation effects, the current new model (PCTBI) has the advantage of accounting for atomistic charge distributions for RNAs. Unlike the previous coarse grained charge models, this new model enables calculations for the ion-mediated electrostatic interactions for RNAs at atomistic charge resolution. For a great variety of biologically significant problems, such as RNA-protein interactions, a detailed account of the charge distribution at the atomistic level is essential. Therefore, the new model would open up the possibility for many important applications such as RNA-protein docking and the kissing complex formation in HCV genomic RNA. The study described here leads to several main conclusions, as shown below.
From a computational point of view, the development of the PCTBI model involves two key ingredients. First, the consideration of ion correlation and fluctuation effects demands sampling of an ensemble of many-ion distributions. The use of atomistic partial charge model renders an increase in computational time. Second, the use of Monte Carlo sampling method significantly improves the computational efficiency for the sampling of ion distribution. The combination of the two strategies enables the PCTBI model to predict electrostatic free energy and ion binding properties with higher resolution and accuracy. Future development of the model involves consideration of ion dehydration effects and further improvement of the sampling method in order to treat larger RNAs.
This research was supported by NIH grants R01-GM117059 and R01-GM063732 (to SJC) and Zhejiang Provincial Natural Science Foundation China under grant LQ14A040004 (to LZS).