Computational Assessment of Potassium and Magnesium Ion Binding to a Buried Pocket in GTPase-Associating Center RNA.

An experimentally well-studied model of RNA tertiary structures is a 58mer rRNA fragment, known as GTPase-associating center (GAC) RNA, in which a highly negative pocket walled by phosphate oxygen atoms is stabilized by a chelated cation. Although such deep pockets with more than one direct phosphate to ion chelation site normally include magnesium, as shown in one GAC crystal structure, another GAC crystal structure and solution experiments suggest potassium at this site. Both crystal structures also depict two magnesium ions directly bound to the phosphate groups comprising this controversial pocket. Here, we used classical molecular dynamics simulations as well as umbrella sampling to investigate the possibility of binding of potassium versus magnesium inside the pocket and to better characterize the chelation of one of the binding magnesium ions outside the pocket. The results support the preference of the pocket to accommodate potassium rather than magnesium and suggest that one of the closely binding magnesium ions can only bind at high magnesium concentrations, such as might be present during crystallization. This work illustrates the complementary utility of molecular modeling approaches with atomic-level detail in resolving discrepancies between conflicting experimental results.

Introduction

The folding and functioning of RNA molecules usually depends upon their interactions with monovalent and divalent cations.[1−3] Diffusing cations are crucial, and in general, they are the major driving force in the stabilization of RNA structures via neutralizing the negative charges on RNA phosphate moieties.[3] However, in some cases, more complex tertiary structures are stabilized through the interaction of ions with specific RNA atoms, including bases and phosphates.[3−6] Depending on the nature of the local RNA structure and the ion type, the interacting ions either “associate” with RNA via hydrogen bonding to their first shell of waters (i.e., water-mediated interactions) or directly bind to one or more RNA atoms.[3] The latter case, which is referred to as “chelation”, is possible only after one or more water residues between the ion and the RNA are removed from the first hydration shell surrounding both the ion and the RNA, a process that involves a high energetic barrier. Consequently, considering the large number of structurally characterized RNA molecules, there are fewer cases of identified direct ion–RNA interactions, as compared to those of associated ions.[7]

Some RNA molecules have been studied extensively by NMR and crystallography as model systems, whose conformational and folding characteristics depend on specific ion binding as well as diffuse ion concentrations.[8−12] A model system that has been studied extensively in terms of its ion interactions is the GTPase-associating center (GAC) RNA (also known as 58-nucleotide RNA and here referred to as GAC RNA), which is a highly conserved fragment of the 23S ribosomal RNA that binds to the L11 ribosomal protein.[13] Melting experiments have suggested that the tertiary structure formation of GAC RNA requires the presence of one bound monovalent ion (NH4+ and K+ preferred)[14,15] and one or two bound divalent ions (with a preference for Mg2+).[16] There are two co-crystal structures of GAC RNA and L11 protein: a 2.6 Å resolution structure of GAC from eubacterium Thermotoga maritima and a 2.8 Å resolution structure of GAC RNA from E. coli in complex with the C terminal domain of L11 protein from Bacillus stearothermophilus, obtained by Wimberly et al. (PDB code: 1MMS)[13] and Conn et al. (PDB code: 1HC8), respectively.[17] The latter structure is the result of re-refinement of a previous structure from the same group (PDB code: 1QA6).[18] In spite of residue differences at few positions, the RNA conformations are very similar in both structures (backbone root-mean-square deviation (RMSD) of 0.6 Å, Figure ). An important feature of the RNA tertiary structure conformation is that a bulged region containing nucleotides A1070, G1071, and C1072 is clamped between a U-turn on its 5′ side, its 3′ tail nucleotides, and a hairpin, C1092–G1099. Clamping this bulged region brings some phosphate and base oxygen atoms in relatively close proximity and results in the formation of a negatively charged area that is expected to interact with cations to maintain its stability. In both structures, three cations directly bind to oxygen atoms in this region, with ion–oxygen distances (IODs) matching direct chelation (Figure and Table ): an Mg2+ ion (Mg167 according to the nomenclature by Misra and Draper[19]) that directly bridges the phosphate oxygen of A1073 and O4 of the U1094; an Mg2+ (Mg163) that binds to phosphate oxygen atoms of A1069 and A1070; and a cation in the buried region between those magnesium ions that directly interact with the bulge phosphate oxygen atoms (Figure ). This third cation is identified as magnesium in the 1MMS structure (Mg58) but as potassium in the 1HC8 structure (K58). According to the 1HC8 structure, the third cation is hypothesized to be the monovalent ion that is essential for final tertiary folding, as characterized in the melting experiments.[14,17] As the RNA can be neutralized more effectively by cations with higher charges, it might be expected that a bound Mg2+ ion will be a better candidate to stabilize the tight ion-binding pocket in the GAC RNA, consistent with the 1MMS structure.[13] However, Conn et al. suggested that this stabilizing and bound ion is a K+ ion based on its distances to the chelating oxygen atoms, and in an experiment of crystallizing the GAC complex with thallium (which has a similar nominal ionic radius to potassium[12]), they found a bound thallium in the same position.[17] However, simultaneous binding of three cations (Mg167, K/Mg58, and Mg163) to RNA in a relatively confined area (Figure ), as depicted in the crystal structures, at full occupancy could be considered unexpected due to electrostatic repulsion, except under very high salt concentrations (i.e., under crystallization conditions). Unfortunately, both crystal structures do not comment on the occupancy. Binding of Mg167 is supported by hydroxyl radical footprinting experiments[20] and nonlinear Poisson–Boltzmann calculations,[19] whereas binding of Mg163 is not supported by noncrystallographic experimental data. Instead, according to nonlinear Poisson–Boltzmann calculations, Misra and Draper concluded that although positions 163 and 167 have larger electrostatic potential compared to that in other magnesium sites in the 1HC8 crystal structure, binding to position 163 is energetically unfavorable due to the ion dehydration cost and repulsive force from the monovalent ion.[19] Therefore, it can be hypothesized that the occurrence of Mg163 near the K58 in 1HC8 or near the Mg58 in 1MMS may be due to crystallization conditions (i.e., very high salt concentrations), crystal packing, partial occupancy, or a combination of all three factors.

Figure 1

Overlay of the GAC RNA crystal structures, 1HC8 (cyan ribbon) and 1MMS (yellow ribbon). The bulged region of A1070, G1071, and C1072 (red atoms and bonds) is clamped between U-turns U1065–A1069 (blue atoms and bonds) and hairpin C1092–G1099 (magenta atoms and bonds). Atoms and bonds of both structures are illustrated. Mg163 and Mg167 (green spheres) have the same identity in two structures, but K/Mg58 (purple sphere) is identified as potassium in 1HC8 and magnesium in 1MMS. Atoms of other residues and hydrogen atoms are not shown for clarity. The molecular graphics were generated using Chimera 1.9.[21]

Figure 2

Potassium-binding pocket according to Conn et al. (PDB entry 1HC8).[17] The potassium is represented as a purple sphere and magnesium ions as green spheres. The nucleic acid backbone is in licorice representation, with phosphorus atoms in yellow, carbon atoms in gray, and oxygen atoms in red. Base atoms are not shown for clarity, except U1094, whose O4 is bound to Mg167. The bound potassium is directly chelated to six phosphate oxygen atoms: A1069@O3′, A1070@OP2, C1072@O5′, OP1, OP2, and A1073@OP1 (@ denotes atom). Mg163 directly bridges A1070@OP2 and C1072@OP1 and is located at a distance of 4.2 Å from the bound potassium. Mg167 is directly bound to A1073@OP2 and U1094@O4 and is located at a distance of 5.7 Å from the bound potassium (OP1 and OP2 are nonbridging oxygen atoms of the nucleotides). The figure was generated using Chimera 1.9.[21]

Table 1

Ion Distances in the Pocketa

pocket oxygen atoms	crystallographic distance according to 1HC8 (Å)	crystallographic distance according to 1MMS (Å)	average distance (Å) in ion-competition simulations	average distance (Å) in ion hydration simulations
A1069@O3′	3.2	3.0	3.7 ± 1.1	3.2 ± 0.2
A1070@OP2	2.7	2.8	3.1 ± 0.9	3.0 ± 0.2
C1072@OP1	2.9	2.9	3.2 ± 0.7	3.0 ± 0.2
C1072@OP2	3.3	3.4	3.2 ± 0.8	3.0 ± 0.2
C1072@O5′	3.5	3.8	4.1 ± 0.6	3.4 ± 0.2
A1073@OP1	2.8	2.6	4.6 ± 1.7	2.7 ± 0.1

Distance between the closest K+ and the oxygen atoms in the monovalent pocket from 20 μs aggregate trajectory of ion-competition molecular dynamics (MD) (COMPETITION_1 set) and 200 ns aggregate trajectory of ion hydration simulations (HYDRATION set) using the 1HC8 crystal structure. Errors are standard deviations.

As reviewed previously, numerous experimental approaches applied to study RNA–ion interactions have led to quantitative measurement of the ion-dependent stability of many RNA conformations.[8] These experimental approaches include monitoring the UV absorbance during thermal unfolding of tertiary RNA conformations,[16,22,23] measurement of the force required to mechanically unfold certain RNA structures using optical tweezers,[24] calorimetric analysis of the ion-binding events inferred to occur during the RNA unfolding/refolding procedure,[25] and fluorescence titrations of proposed ion-binding events at different temperatures.[26] One potential limitation of these experimental methods is the difficulty of separating the binding to a single specific binding site from interactions between the background ion media and the RNA in titration experiments.[20] A further confounding difficulty in both approaches is that the ion-binding events are usually accompanied with major conformational changes in RNA.[20] In other words, it is difficult to separate thermodynamics of ion binding from thermodynamic changes resulted from folding.[20,27] MD simulations can be a complementary approach to help overcome these limitations through the ability to separately investigate each and every ion-binding event at a very detailed atomistic level; however, they have limitations of the applied force field and effective sampling. In silico molecular models enable us to artificially separate the influences by examining the ion-binding event in a limited conformational space close to the RNA native structure and compare the intrinsic free energies of binding for different ions.

In this work, we conducted MD simulations in explicit solvent under different ion conditions, that is, a virtual ion-competition experiment, to analyze the preferential ion binding at the K/Mg58 position depicted in two crystal structures of 1MMS and 1HC8. To further characterize preferential binding at this site, we also employed the umbrella sampling method to estimate the potential of the mean force (essentially the free energy) of pulling specific ions out of the pocket to bulk solvent at various salt concentrations. The umbrella sampling studies were performed in the presence and absence of a magnesium ion in the Mg163 position to analyze the influence this magnesium has on the ion binding in the adjacent K/Mg58 pocket as well as to examine the hypothesis that it can coexist close to the other two cations at both the high magnesium concentrations used in crystallization and the lower magnesium concentrations used in most other solution biochemical assays. Additionally, we perform a characterization of the topology of the electron density of both K+ and Mg2+ present in the binding pocket using the quantum theory of atoms in molecules.[28] This analysis provides a quantum mechanics (QM)-based study of the intermolecular bonds of our systems and a solid analysis of the interaction properties.

Methods

MD Simulations to Investigate Preferential Ion Binding in the Pocket

COMPETITION_1 Simulation Set

The crystal structure of GAC RNA with PDB code 1HC8(17) was used as the initial structure. All of crystallographic ions and water residues, as well as the terminal GTP residue, were removed. The terminal GTP is paired with a uracil in the crystal structure and was removed as it is located at about 30 Å from the monovalent ion-binding site and does not closely interact with this site in the absence of crystal packing. Hydrogen atoms were added to the remaining 57 nucleotides, and the in vacuo topology and coordinates were built using the tLEaP program in Amber 12[29] with ff12SB force field, which includes ff99[30] with updated X[31] and α/γ parmbsc0[32] modifications. A total of 2500 cycles of the steepest descent minimization with the Hawkins–Cramer–Truhlar Generalized Born implicit solvent model[33,34] was performed using SANDER program in Amber 12. The resulting structure (with an RMSD of 0.37 Å to the initial structure) was solvated in a truncated octahedral TIP3P water[35] box such that the minimum distance of the residues to the edge of the box was 12.0 Å. Then, 28 Mg2+ ions with Allnér–Nilsson–Villa parameters[36] and 28 K+ and Cl– ions with Joung–Cheatham parameters[37] were added to the box using tLEaP. The CPPTRAJ program[38] in Amber Tools was used to randomize ion positions in such a way that the distance of each ion to any RNA atom was not less than 6 Å and ions were separated with a minimum distance of 4 Å, causing the initial structures to lack ions in the positions previously occupied in the crystal structures. Ion position randomization was done 20 times in successive steps to generate 20 independent copies of the system. To avoid any initial bias for ion occupancy in the K/Mg58 pocket, the distances of the closest K+ and Mg2+ ions to the pocket were identified in each copy to make sure that almost half of the copies have each ion types as closest to the pocket.

Each copy of the system was equilibrated in four steps: (1) minimization, (2) heating at constant volume (NVT), (3) further heating at constant pressure (NPT), and (4) restrained MD at constant pressure. Details of the equilibration protocol are described in the Supporting Information as EQ1. The production MD simulations were performed using PMEMD.CUDA from Amber 12 for 1 μs per copy with 2 fs time step in the NPT ensemble, whereas RNA atoms were restrained with a 0.5 kcal/mol Å2 force constant to assess the ion binding, with the RNA structure kept close to the crystal conformation. The temperature was held at 298.15 K with weak-coupling algorithm, with a time constant of 10 ps.[39] SHAKE[40] was used to constrain bonds to hydrogen, and particle mesh Ewald[41] with an 8.0 Å cutoff and default parameters was used to treat long-range electrostatics in the production simulations. The coordinates were saved to trajectory files every 5000 steps (10 ps intervals).

COMPETITION_2 Simulation Set

The last frames of the 20 independent COMPETITION_1 simulations were used as initial structures for the COMPETITION_2 simulations, which were run with the same protocol, except with no restraints on the RNA for 500 ns if the ion-binding preference is reproducible with potential subtle changes in RNA conformation.

HYDRATION Simulation Set

To analyze the effect of the ion position on its hydration shell, 20 copies of simulations with the same RNA/ions composition as above were built (named here as ion hydration simulations) except that K58, Mg163, and Mg167 were positioned in their crystallographic sites. These three ions as well as the RNA atoms were fixed with 0.5 kcal/mol Å2 restraints during the course of 10 ns simulations, with the same positional restraint used for the RNA atoms. These simulations were done to assess the ion hydration shell in the pocket exactly in the same pocket conformational situation as that of the crystal structure. The HYDRATION simulation was also repeated with Mg2+ in the K/Mg58 site. A representative structure from both of these HYDRATION simulations was extracted with the 50 closest waters to the K/Mg58 ions to be used in quantum calculations.

Ion Pull-Out Umbrella Sampling Simulations

Umbrella sampling was done to quantitatively compare the free energy of binding (potential of mean force (PMF)) of different ions to that of the GAC monovalent ion-binding site at three different ionic environments. Before the umbrella sampling simulation, we performed unrestrained simulations of GAC RNA in three different ion environments to identify spontaneous ion leaving or unbinding events as follows:

MG0: RNA + 11 crystallographic Mg2+ + neutralizing K+ + 100 mM excess KCl

MG15: RNA + 11 crystallographic Mg2+ + neutralizing K+ + 100 mM excess KCl + 15 mM excess MgCl2

MG80: RNA + 11 crystallographic Mg2+ + neutralizing K+ + 100 mM excess KCl + 80 mM excess MgCl2

To prepare models, the crystal structure with PDB code 1HC8 was used as the initial structure.[17] Crystallographic K+ and Mg2+ ions as well as water residues were retained but the terminal GTP residue was removed (as discussed previously). Hydrogen atoms were added to the remaining 57 nucleotides. Using tLEaP program, the structure was solvated in a truncated octahedral box of TIP3P waters such that the minimum distance of the residues to the edge of the box was 12.8 Å. Then, 33 K+ ions were added to neutralize the whole system (RNA + crystallographic ions) and 26 K+ and 26 Cl– ions were added to provide ∼100 mM excess KCl. To prepare three sets of simulations under different ionic conditions, 0, 4, and 21 Mg2+ ions as well as 0, 8, and 42 Cl– ions were added to the systems to make approximately 0, 15, and 80 mM excess MgCl2 concentrations, respectively (here referred to as MG0, MG15, and MG80 sets). The noncrystallographic ions were put in random positions using CPPTRAJ, as described in the COMPETITION_1 section to prepare three (for the MG0 set) and six (for MG15 and MG80 sets) system copies. The MG0 set was equilibrated with EQ1 protocol discussed in the Supporting Information. The MG15 and MG80 sets were equilibrated in nine steps of minimizations and restrained MD simulations using a newer equilibration protocol, as described in the Supporting Information as EQ2. Hydrogen mass repartitioning was used to facilitate a 4 fs time step in the final step of the equilibration as well as in the production simulations.[42,43] This was done by transferring part of mass of the heavy atoms of the solute to the covalently bonded hydrogen atoms using the parmed.py program from Amber 14.[44] For each copy of each set, the equilibration set was followed by ∼1 μs of production MD using the PMEMD.CUDA program from Amber 14, as described for the production step in the former section but with no restraints.

The reaction coordinate for the umbrella sampling was chosen to be the distance between the ion and the center of mass of nine phosphate groups of A1070, C1072, and A1073 (i.e., the phosphorus and two nonbridging oxygen atoms of each residue). We refer to this position as the reaction coordinate anchor point. Initial structures for the umbrella sampling windows were chosen differently in three sections along the ion-leaving path as follows:

∼2.5 to 6.9 Å: The potassium bound to the monovalent ion-binding site was tracked during the above-mentioned unrestrained simulations (MG0, MG15, and MG80), and all ion-exchange events were analyzed both visually and using CPPTRAJ. For each simulation set, trajectory frames from one ion-exchange event with the least conformational change in RNA near the ion-binding site were chosen for generating the initial structures. The structures for initial umbrella windows with reaction coordinate values between the initial equilibrated distance (∼2.5 Å) and 6.9 Å were extracted from the chosen ion-exchange path trajectory using CPPTRAJ.

0.5 to ∼2.5 Å: The last frame of the umbrella with the reaction coordinate value of the initial distance in the ion-exchange event was used as the initial structure for a set of consecutive 2 ns MD simulations to direct the ion to 0.5 Å on the reaction coordinate. The umbrella sampling simulations in this part of the path with low reaction coordinate values were conducted in a serial fashion, with the last frame of each window simulation used as the initial frame of the next window simulation with a 0.1 Å lower reaction coordinate value. Each of these umbrellas was started with short equilibration simulations described as EQ3 in the Supporting Information.

7.0 to 20.0 Å: The initial frame of the umbrella window with reaction coordinate value of 6.9 Å was used as the initial structure of a short MD simulation, in which the leaving K+ ion was gradually directed to the bulk solvent (reaction coordinate = 20.0 Å) in 130 ps as follows: the pseudoangle involving the ion and the phosphorus atoms of A1071 and C1072, as well as the pseudodihedral involving the ion and the phosphorus atoms of A1073, G1071, and C1072, was kept fixed at the initial value with 200 kcal/mol Å2 restraints, and the reaction coordinate value was increased from 6.9 to 20.0 Å with a 100 kcal/mol Å2 distance restraint. All other RNA heavy atoms were restrained in position with 20 kcal/mol Å2 force constant to avoid potential conformational disruption resulting from the ion displacement near RNA atoms. Then, the initial structures for umbrellas within this distance frame were extracted using CPPTRAJ as well to generate the path from 7.0 to 20.0 Å. The window interval for the entire reaction coordinate was 0.1 Å.

Initial structures generated in sections 2 and 3 of the path were equilibrated using the EQ3 protocol described in the Supporting Information. Production MD was performed for 2 ns at constant pressure and temperature for each umbrella along the whole path, during which a Langevin thermostat[45] was used with collision frequency of 5 ps–1 to hold the temperature at 298.15 K; the SHAKE algorithm[40] was used to constrain the bonds involving hydrogen atoms; the direct space cutoff was set to 8.0 Å; and particle mesh Ewald[41] was used to calculate long-range interactions. A 2 fs time step was used for the production simulations, and the reaction coordinate values were saved for every step. During equilibration and production simulations, the ion was held at the related umbrella value with 100 kcal/mol Å2 distance restraint. Section 2 of the path was sampled in a serial fashion, whereas other sections were sampled in parallel as independently running simulations. This way, the reaction coordinate between 0.5 and 20.0 Å was sampled at 0.1 Å intervals for a leaving K+ at 0, 15, and 80 mM excess Mg2+ concentrations with 2 ns simulation time per umbrella (US_MG0_K, US_MG15_K, and US_MG80_K simulation sets).

To examine the behavior of Mg2+ in the same position, the leaving ion was mutated from K+ using parmed.py program in Amber 14, and the same initial coordinates for each umbrella window were used. When mutating the ion to Mg2+, a Mg2+ located in bulk solvent was simultaneously mutated to K+ to keep the system charge neutral. This way, the reaction coordinate between 0.5 and 20.0 Å was sampled at 0.1 Å intervals for a leaving Mg2+ at 0, 15, and 80 mM excess Mg2+ concentrations with 2 ns simulation time per umbrella (US_MG0_MG, US_MG15_MG, and US_MG80_MG simulation sets). For no excess Mg2+ scenario, the reaction coordinate was also sampled for a leaving Na+ ion by mutating the K+ to Na+ (US_MG0_NA set). The Na+ simulation was performed to compare the affinities of the pocket for Na+ and K+ as a test for experimental agreement.[14,15]

Two additional umbrella sampling simulations were conducted for excess Mg2+ concentrations of 15 and 80 mM for either Mg58 or K58 leaving the binding pocket, where Mg163 was moved to random positions far from the RNA using CPPTRAJ in the initial structure of each umbrella (US_MG15_K_163, US_MG80_K_163, US_MG15_MG_163, and US_MG80_MG_163 sets). To examine the effect of magnesium parameters on its free energy of binding, the Lennard-Jones parameters of all magnesium residues in the K+ → Mg2+ mutated system in the first scenario were changed to the parameters developed by Åqvist[46] as well as hydration free energy (HFE) and IOD parameters developed by Li et al.[47] (US_MG0_MG_AQV, US_MG0_MG_HFE, and US_MG0_MG_IOD sets). All simulation sets are summarized in Table S1.

Histograms of the reaction coordinate values for the umbrellas were analyzed using CPPTRAJ, and extra umbrellas were initiated where the overlap between umbrellas was poor. The WHAM method,[48,49] as implemented by Grossfield (WHAM version 2.0.9), was used to generate the PMF from the umbrella sampling simulations. Other trajectory analyses, such as measurement of RMSD and distances and calculation of the number of waters in hydration shells, were done using CPPTRAJ in Amber 14 (details are available in the Supporting Information).

Quantum Calculations

Using the representative structure of the system with K+ from a cluster analysis of the explicitly solvated MD simulations (HYDRATION simulations), a single point calculation was obtained using density functional theory at the M06-2X[50,51] level and a 6-31G(d) basis set.[52] The calculated wave function was used to extract the electron density and its topological properties with the AIMAll suite of programs.[53] To calculate the energy difference between the two ions, the same solvated system was calculated switching K+ for Mg2+ and the properties of the electron density were extracted using the same methodology.

Results

MD Simulations Exploring Preferential Ion Binding in the Pocket

To assess the occupancy of the monovalent binding site, MD simulations were performed with different ion compositions to probe preferential ion binding. The GAC RNA was solvated in TIP3P water with total added 28 Mg2+ (neutralizing amount) and 28 K+ ions (plus 28 neutralizing Cl–). The ions were initially located at random positions within the simulation box with the minimum distance of 6 Å from any RNA atom and 4 Å from each other. Because of the low exchange rates for Mg2+ ions,[54] it is not possible to converge the Mg2+ exchange events with a single simulation, given the achieved simulation time scales. Therefore, we increased the number of simulation copies and simulation length per copy to ideally improve the sampling within the capabilities of available resources. A total of 20 copies of the simulation were generated with initial placement of the ions at different random positions. A 1 μs production simulation was performed for each copy, as described in the Method section. To better understand the effect of the ion binding on the tertiary structure, as opposed to the folding, the RNA crystal conformation was restrained with 0.5 kcal/mol Å2 positional restraints on all RNA atoms (COMPETITION_1 set).

Figure illustrates the proximity of the ion to the ion-binding pocket by plotting the distance of the closest potassium and magnesium to the center of mass of the phosphate atoms of A1070, C1072, and A1073 in the tight ion-binding pocket. The distance of the bound potassium to this center of mass point is 0.97 Å in the crystal structure (1HC8), which, as can be seen in Figure , is not rigidly maintained in the simulations due to fluctuations in the position of the ion and the residues in the binding pocket. In most cases, an empty binding pocket became occupied by potassium during the equilibrium stages, although due to the ion randomization in the initial structures, 8 copies had closer Mg2+ ion to the pocket and 12 copies had closer K+ ion to the pocket in their initial structures and the distances to the closest ions of different types were very similar. However, among 20 copies, the pocket was occupied by magnesium during only one simulation (shown as cyan, copy 1). This ratio implies that the occupancy of the pocket by K+ is preferred over that by Mg2+, but the difference between the free energies of binding of two ions is not high enough to completely prevent Mg2+ from binding.

Figure 3

Distance of the closest K+ (top) and Mg2+ (bottom) to the center of mass of the phosphorus and two nonbridging oxygen atoms of A1070, C1072, and A1073 (named as the reaction coordinate anchor point) during the ion-competition simulations with restrained RNA (COMPETITION_1). Each line represents running-averaged data over 10 ns windows for one of the 20 copies. The blue dashed line represents the crystallographic distance to the bound potassium. The arrows with copy numbers are markers for points of interest that are further discussed in the text.

Simultaneous occupancy of K58 and Mg163 positions occurred in only two copies (shown as black and indigo, copies 2 and 3, respectively). Interestingly, in one of these copies (black, copy 2), the chelated magnesium stayed closer to the K58 pocket and the potassium pocket remained only partially occupied by potassium that remained close to the pocket but far from the bound magnesium. The buried pocket was also partially occupied in the simulation copy indicated as brown (copy 4). The violet copy (copy 5) in the bottom of Figure indicates occupancy of Mg167 site by magnesium while the buried pocket was occupied by potassium. In summary, the buried ion-binding pocket was occupied by magnesium in only one of the simulation copies and only two other copies showed simultaneous occupancy of K58 and Mg163 (Figure S1).

Similar ion-competition simulations were repeated for 500 ns but without any restraints on the RNA to detect any potential differences in the occupancy of the binding pocket caused by subtle RNA conformational changes (COMPETITION_2 set). The results (Figure S2) show almost similar statistics for ion occupancy in the K58 and Mg163 sites although there is more flexibility in the position of the bound K, which is reasonable due to the greater flexibility of the RNA conformation.

The bound K+ is located within a distance of 2.7–3.5 Å from six oxygen atoms in the crystal structure,[17] as illustrated in Figure . Figure shows the distributions of the ion distances in the pocket to these oxygen atoms in the COMPETITON_1 simulations to check whether the ion retains its exact crystallographic position within the pocket or shifts position toward the sides. According to the short-distance major peaks in Figure , the location of the bound K+ during the RNA-restrained simulations was similar to that of the crystal structure with some flexibility in the binding pocket. Sharper short-distance peaks for OP2 atoms of A1070 and C1072 compared with wider and bimodal distributions for other four atoms imply the tendency of the bound K+ to shift toward the Mg163 binding site within its own pocket while Mg163 is absent in its crystallographic position. This shift is also seen in Table , which compares the crystallographic distance of the K+ in the pocket to that of the pocket oxygen atoms, with average of those distances generated by the simulations. The position of the bound K+ in the unrestrained competitive simulations is also similar to that of restrained simulations (Figure S3) although the peaks for K+ distance from O3′ of A1060 and O5′ of C1072 are less sharp.

Figure 4

Histograms of distances of the bound K+ in the pocket to the six crystallographic interacting oxygen atoms in 20 μs aggregate trajectory of RNA-restrained ion-competition simulations (COMPETITION_1).

The absence of a bound Mg2+ in the position of Mg163 in most simulation copies and the simultaneous shift of the bound potassium toward this position in these simulations imply that the crystallographic occurrence of the Mg163 might be a result of the crystal conditions. When the Mg163 position lacks Mg2+, the bound potassium might tend to shift toward this position rather than being in the middle of the buried binding pocket. However, the possibility of K+ being in different positions within the binding pocket cannot be justified on the basis of the bimodal distribution of ion distances from C1072 and A1073. In fact, such observation can be a force field and/or convergence artifact related to the simulations. Currently, no experimental evidence is available to support it.

In the crystal structure,[17] the K+ binding pocket is isolated from the solvent by RNA. Therefore, it is feasible that the ion is partially or completely dehydrated. No water density is visible around the bound K+. In solution and in simulations with TIP3P waters, K+ has five or six waters in its first hydration shell. During the ion-competition simulations, the bound potassium was often only partially dehydrated with three or four water residues in its first solvation shell (Figure ), implying that the pocket can accommodate a potassium ion together with some of its first-shell waters, allowing the bound potassium to move. The crystal structures do not contain any waters in the pocket, but fully dehydrating a potassium or magnesium upon binding is energetically costly and hence very unlikely.

Figure 5

Distribution of waters in the first solvation shell of the bound K+ in 20 μs aggregate trajectory of ion-competition simulations and 200 ns aggregate trajectory of ion-hydration simulations.

For further investigation of the role of hydration, we set up 20 copies of short MD simulations (10 ns) using the 1HC8 structure to monitor the possible diffusion of water molecules into the pocket (HYDRATION set). The RNA and K58, Mg163, and Mg167 were restrained to their crystallographic positions. The results show that (Table ) the water molecules can leak into the pocket and directly interact with K+ (Figure ). However, the comparison between the two histograms in Figure reveals that the greater chance of having more waters in the pocket and around the ion in ion-competition simulations is probably a result of ion relocation from the exact crystallographic position in the pocket. Yet, the presence of three to four water residues in both simulations supports our hypothesis that the ion in the pocket is partially dehydrated even under its crystal conditions (noting that the lack of ion hydration in the crystal could be simply due to its omission during the refinement).

Ion Pull-Out Umbrella Sampling

Although the observed preference of potassium over magnesium in occupying the binding site in these simulations may be caused by thermodynamic favorability of potassium in the pocket; one cannot rule out the possibility that the higher kinetic barrier for magnesium dehydration has prevented it from populating the pocket within the simulation time scale. To distinguish between these scenarios, we calculated the free energy needed to pull out different ions from the potassium-binding site to the bulk solvent using umbrella sampling (see Method for details). The crystallographic magnesium ions were retained in all experiments, and 0 Mg2+, 15 mM Mg2+, and 80 mM Mg2+ were added to independent MD simulation experiments, of which the 80 mM Mg2+ is similar to the ion environment used in the crystallization of the 1QA6 and 1HC8 structures.[18]

A series of MD simulations were performed during which the ion (Na+, K+, or Mg2+) was restrained to set positions along a reaction coordinate. The reaction coordinate is the distance of the ion to the center of mass of nine phosphate groups of A1070, C1072, and A1073 (phosphorus and two nonbridging oxygen atoms of each residue), referred to hereafter as the “anchor point”. The crystallographic distance on this reaction coordinate for the bound K+ is 0.97 Å in the 1HC8 crystal structure.

To prepare the initial structures for each umbrella, unrestrained MD was performed in three to six copies for each system (MG0, MG15, and MG80 sets) to capture ion-exit events as the inserted ion in the pocket leaves the RNA. One ion exit was observed during 1 μs simulation in two copies of MG0 and four copies of MG15 and MG80.

The trajectory frames in which the leaving K+ was at a distance of less than 6.9 Å from the anchor point were extracted from a given simulation to be used as initial structures for related windows. Beyond 6.9 Å, the ion was free, well exposed to the solvent, and directed to 20 Å (bulk solvent) with a 130 ps simulation. Frames from this simulation were used to obtain the initial umbrella structures from 6.9 to 20 Å. The umbrellas were spaced 0.1 Å from each other. A total of 2 ns of MD simulation was performed for each umbrella. The ion was restrained using a harmonic distance restraint with 100 kcal/mol Å2 force constant for each umbrella. Such a large restraint was required especially at some distance ranges close to the RNA to keep the ion in the desired window. The ion was also directed from the lowest umbrella to 0.5 Å on the reaction coordinate with a series of 2 ns umbrella simulations, in which the last frame of a simulation was used as the initial frame of the next simulation.

This GAC RNA adopts its tertiary structure in 0.1 mM Mg2+ with 100 mM K+; 5 mM Mg2+ with 100 mM Na+; or 5 mM Mg2+ with 100 mM NH4+.[55] Here, we replace K+ in the crystallographic pocket with Mg2+ and Na+ to look directly at the free energy of different ion occupancy in the pocket. The sampling was repeated with mutating the leaving K+ to Mg2+ (in all scenarios) and Na+ (in the scenario with no excess Mg2+). According to the experiments done by Wang et al., not only a specific monovalent ion interaction is needed for folding of a similar construct of the GAC RNA, but K+ stabilizes its tertiary structure better than Na+.[14,15] Although the GAC RNA used in those experiments was the E. coli version with the A1061U mutation compared to the one we simulate here (A1061U), the trend of monovalent ion binding to both should be similar, as the pocket structure and RNA global conformation of the A1061U mutant (1MMS structure) is very similar to that of the GAC in 1HC8 structure (Figure ). Therefore, pulling out Na+ was done in this work to see if the same trend between K+ and Na+ is reproducible for this pocket and support the hypothesis that this pocket is the actual monovalent ion-binding pocket, which is expected from those experiments.

Also, to examine the effect of Mg163 on the binding of each ion to the monovalent ion pocket, the Mg163 was moved to a far random position in the initial structure of each window and the simulations were repeated. The exception was in the system with no excess Mg2+, in which Mg163 had already dissociated before the K+ ion-exchange event.

In the resulting PMF plots (Figure ), the free energy of binding for each ion can be calculated as the difference between the PMF values at 20 Å and in the pocket-related well at <3 Å. All plots use a common and arbitrary zero of free energy of the free ion at 20.0 Å from the ion-binding pocket. Slightly different minima on the reaction coordinate for different ions highlights the fact that they have different binding modes within the pocket. The PMFs for binding in the absence of Mg163 suggest that the K+ ion binds tighter with the increase in the concentration of magnesium. However, this is likely an artifact of the zero of the PMF; instead, likely it is destabilization of the unbound state (rather than tighter binding) due to the increased ionic density.

Figure 6

PMF profiles for pulling out K+, Mg2+, and Na+ from the GAC RNA monovalent ion-binding site. Reaction coordinate represents the distance between the ion and anchor. The PMF plots were generated using the second 1 ns window in each umbrella. All of the PMF plots were arbitrarily shifted to zero relative to their magnitude at 20.0 Å.

PMF calculations show that binding of Mg2+ to the monovalent binding pocket is always energetically less favorable than binding of K+. Binding of Mg2+ is energetically less favorable than binding of Na+ in the case with no excess Mg2+. These data support the hypothesis that the pocket is specific for monovalent ions[17] and agree with the experimental observations that K+ stabilizes the GAC RNA tertiary structure better than a bound Na+.[14,15]

In the presence of Mg2+ in Mg163 position, it is unfavorable for Mg2+ and K+ ions to bind to the pocket. This agrees with the observation in ion-competition simulations showing that crystallographic localization of Mg163 and K58 is mutually exclusive. However, in 80 mM Mg2+ simulations, binding of potassium inside the pocket and in the bulk solvent is almost isoenergetic (when Mg163 occupies the crystallographic binding site).

In the absence of Mg2+ in Mg163 position and with no excess Mg2+ ions, the binding of potassium to the pocket is about −2 kcal/mol, which is favorable; the binding of sodium is almost isoenergetic inside the pocket and in the bulk water; and there is about +3 kcal/mol energetic penalty for the binding of magnesium. However, the free-energy difference for having K+ versus Mg2+ in the absence of Mg2+ in the Mg163 position decreases with increasing bulk Mg2+ concentration. Therefore, mixed occupancy of K+ and Mg2+ in this pocket in the whole ensemble seems possible especially at higher Mg2+ concentrations.

As the binding site is buried inside the RNA, the ions do not go through a smooth path when leaving it and so the PMF plots look rough. The small bumps on the PMF plots represent transient interaction of ions with other RNA atoms along their path. There are several other binding modes near the monovalent ion-binding site at 3 and 4 Å on the reaction coordinate. These binding modes can be described as chelation of the ions to a nonbridging oxygen atom of U1097 on their path to the bulk solvent. However, when comparing the PMF plots of each ion type in different scenarios, it should be considered that the paths were generated from independent simulations, and therefore, the same points on the same ion type plots in different scenarios might not exactly represent the same ion–RNA interactions.

We expect the free-energy values to vary to some extent by using different ion models.[36,46,47] To test how the magnesium binding to this pocket might be affected by different magnesium models, a magnesium-pulling experiment was repeated. We chose conditions of no excess magnesium scenario and using crystallographic initial structures. Mg2+ parameters included the traditional Åqvist parameters[46] and HFE and IOD parameters developed by Li et al.[47] The two latter sets are designed to reproduce correct HFE and correct IOD, respectively. Ion dehydration energy and ion–oxygen interactions are two determining elements in ion binding to the GAC pocket.[17]

PMF plots for these four magnesium models are shown in Figure . The HFE model causes an increase in the dehydration barrier for the ion to enter the binding site without a significant change in the binding free energy. The IOD and Allnér–Nilsson–Villa models give very similar results. Åqvist parameters have a small effect on the dehydration barrier but enhance the binding energy by about −1 kcal/mol. The calculated PMFs for the different Mg2+ models vary in magnitude but none exceed the binding free energy of K+ ions.

Figure 7

PMF profiles for pulling out Mg2+ from the GAC RNA monovalent ion-binding site with different magnesium models. PMF plots are calculated under the same conditions described in Figure .

Electron Density Analysis of the Ion Binding to the Buried Pocket

As a complementary approach to study the energetic preference of either K+ or Mg2+ in the buried binding pocket, we conducted a topological analysis of the electron density using the quantum theory of atoms in molecules (QTAIM) methodology. Because we are using a fixed-charge force field in our MD studies that does not permit the redistribution of charge, the use of QTAIM allows us to quantitatively study the binding properties of each ion on the basis of robust QM calculations. QTAIM analysis shows where the electron density is accumulated between two atoms;[28] in this case, the bond paths are located between the ion, three water molecules (bond paths toward the oxygen atoms obtained from HYDRATION simulations), and three oxygen atoms from the phosphate groups of the RNA backbone (A1073, A1069, and C1073, refer to Figure ). The bond critical point (BCP) is a point along the bond path where the shared electron density reaches a minimum, and it is an approximate measure of the amount of electron density built up in the bonding region. The BCP analysis will aid in the characterization of the interactions with the two different ions inside the binding pocket. The values of the BCP found for the six interacting atoms are shown in Table . Overall, the values are ∼0.05 au higher when the K+ ion is present in the binding pocket. The BP and its associate BCP are an order of magnitude higher in the interaction with the oxygen of the phosphate group of the cytosine C1072. This suggests that the presence of K+ imposes a stronger binding in the binding pocket than that of Mg2+.

Table 2

Bond Critical Point Values between the Six Coordination Bond Paths Extracted Using the Most Representative Structure from the HYDRATION MD Simulationsa

	equatorial coordination				axial coordination
	OP1 from A1073	OP2 from A1069	O from WAT	O from WAT	OP1 from C1072	O from WAT
K+	0.07685	0.071889	0.073829	0.067042	0.10253	0.071924
Mg2+	0.03741	0.033866	0.035134	0.032128	0.05013	0.034016

BPs are present in either oxygen atoms of water molecules (marked as O from WAT) or oxygen atoms from phosphate groups (marked as OP1 or OP2 from ID of residue). Values are in atomic unit.

Electron density redistribution of the binding pocket with either K+ or Mg2+ could affect the energetics of the ion binding. To understand this process, we calculated the electron density population, N(A), for the six points of coordination around both K+ and Mg2+ ions (values are presented in Table S2). Related work using MD and a polarizable force field reported the lack of charge redistribution of monovalent ions versus divalent ions.[56] Atomic net charges of the cations in the binding pocket as calculated under the QTAIM method are 0.84 and 1.789 for K+ and Mg2+, respectively. Overall, the presence of K+ ion in the binding pocket increases the electron density of the coordinated water molecules, although the ion only contributes in 0.848772e– and serves as a charge “bridge”. The surplus charge is extracted from the three backbone phosphate oxygen atoms present in the pocket. A different picture occurs with Mg2+, which yields 1.789e– that gets redistributed toward mainly the water oxygen atoms.

The presented QM analysis only provides an initial exploration of the intricate electronic interactions in the studied systems. More elaborate QM calculations, beyond the scope of this work, could provide a deeper insight into the electron density of the binding pocket.

Conclusions

Both ion-competition MD and ion pull-out umbrella sampling simulations agree that binding of potassium is preferred over magnesium in the binding site (Figures and 6). The electron density calculations also confirm the preference of potassium in this pocket, although they only provide a qualitative insight into the charge-transfer effects of such a large biological system, which likely requires a more thorough QM investigation. The simulations that facilitate competition between K+ and Mg2+ show very few instances of Mg2+ binding. However, Mg2+ is not completely blocked from the monovalent pocket, which is in agreement with the PMF result, as the difference of free energy of binding between K+ and Mg2+ is not huge in the absence of Mg163. This supports the identification of potassium in this pocket as described by Conn et al.[17] rather than magnesium as suggested by Wimberly et al.[13] Also, preference of K+ over Na+ in this pocket in low Mg2+ concentration ion pull-out simulations agrees with the trend in melting experiments observed by Lu and Draper and Wang et al., where K+/Mg2+ stabilized the GAC tertiary structure more efficiently than other monovalent and divalent ions.[14,55] This agreement further supports the hypothesis that the K58 position is the monovalent ion-binding site.

Because Wimberly et al. have not discussed the details for Mg2+ identification in the K/Mg58 position as their crystal structure shows, it is possible that the ion has been misidentified in that position. However, more elaborate experiments have led to accurate identification of K+ in that position later by Conn et al.[17] Further support for such conclusion comes from the fact that the closest coordinating oxygen atom in 1MMS (OP1 of A1073)[13] is located 2.6 Å away from the chelated Mg2+ (Table ), whereas the mean experimental value for Mg2+–O distances in RNA structures is 2.08 Å according to the Cambridge Structural Database[57,58] and distances longer than 2.4 Å are suspected to be misidentified.[58] However, according to the results of the ion pull-out simulations and the free-energy difference between K+ and Mg2+ for their binding to the pocket (Figure ), it is not surprising if magnesium is observed in the pocket at higher magnesium concentrations, which are used under crystallization conditions. However, it is important to note that the free-energy difference would dictate potassium binding at lower concentrations such as those that have been used in melting experiments. Therefore, it is very likely that the chelated monovalent ion described in those experiments have been actually bound to this pocket.[14,15]

The MD simulations suggest that occupancies of the Mg163 and K58 positions by related ions are mutually exclusive at lower Mg2+ concentrations. Occupancy of the Mg163 site in the presence of potassium in the neighboring pocket is more likely to be caused by high magnesium concentration under crystallographic conditions.

In the cell, K+ concentration is typically estimated at ∼150 mM and free Mg2+ at 1–2 mM. Our pulling experiments suggest that under these solution conditions, the Mg163 site will be unoccupied, whereas K58 is stably bound. Regardless of magnesium concentration, Mg163 occupancy, and the chosen magnesium force field, potassium binding is always more favorable than magnesium binding to the monovalent binding site. In the absence of an Mg2+ ion, the pocket will be more accessible to water molecules. Mobile water molecules could relax the pocket geometry, allowing higher mobility of the nucleotides and also exchange of the bound K+.