the folding and functioning 
 of rna molecules usually depends upon 
 their interactions with monovalent and divalent cations .
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 . however , in some cases , more complex tertiary structures 
 are stabilized through the interaction of ions with specific rna atoms , 
 including bases and phosphates . 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 .
, 
 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 .
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 .
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 . melting 
 experiments have suggested that the tertiary structure formation of 
 gac rna requires the presence of one bound monovalent ion ( nh4 and k preferred ) and one or two bound divalent ions ( with a preference for mg ) .
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 .
the latter structure 
 is the result of re - refinement of a previous structure from the same 
 group ( pdb code : 1qa6 ) . 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 1 ) .
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 , c1092g1099 . 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 2 and table 1 ) : an mg ion ( mg167 according to the nomenclature by misra and draper ) that directly bridges the phosphate oxygen 
 of a1073 and o4 of the u1094 ; an mg ( 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 2 ) .
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 .
as the rna can be neutralized more effectively by cations with higher 
 charges , it might be expected that a bound mg ion will 
 be a better candidate to stabilize the tight ion - binding pocket in 
 the gac rna , consistent with the 1mms structure .
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 ) , they found a bound thallium in the same position .
however , simultaneous binding of three cations 
 ( mg167 , k / mg58 , and mg163 ) to rna in a relatively confined area ( figure 2 ) , 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 ) .
binding of mg167 is supported 
 by hydroxyl radical footprinting experiments and nonlinear poisson  boltzmann calculations , 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 .
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 .
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 u1065a1069 ( blue atoms and bonds ) 
 and hairpin c1092g1099 ( 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 .
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 ) .
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 .
ion interactions have led to quantitative measurement 
 of the ion - dependent stability of many rna conformations .
these experimental approaches include monitoring 
 the uv absorbance during thermal unfolding of tertiary rna conformations , measurement of the force required 
 to mechanically unfold certain rna structures using optical tweezers , calorimetric analysis of the ion - binding events 
 inferred to occur during the rna unfolding / refolding procedure , and fluorescence titrations of proposed ion - binding 
 events at different temperatures .
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 .
a further confounding difficulty in both approaches 
 is that the ion - binding events are usually accompanied with major 
 conformational changes in rna .
in other 
 words , it is difficult to separate thermodynamics of ion binding from 
 thermodynamic changes resulted from folding .
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 mg present in the binding pocket using the 
 quantum theory of atoms in molecules .
this analysis provides a quantum mechanics ( qm)-based study of the 
 intermolecular bonds of our systems and a solid analysis of the interaction 
 properties .
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 with ff12sb force field , 
 which includes ff99 with updated x and / parmbsc0 modifications . a total of 2500 cycles of the steepest descent 
 minimization with the hawkins  cramer 
truhlar generalized 
 born implicit solvent model 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 box such that 
 the minimum distance of the residues to the edge of the box was 12.0 
  .
then , 28 mg ions with allnr  nilsson  villa 
 parameters and 28 k and cl ions with joung  cheatham parameters were added to the box using tleap . the cpptraj 
 program 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 mg 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  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 .
shake was used to constrain bonds to hydrogen , and particle mesh ewald 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 ) .
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 . 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  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 mg 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 .
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:(1)mg0 : rna + 11 crystallographic 
 mg + neutralizing k + 100 mm excess kcl(2)mg15 : rna 
 + 11 crystallographic 
 mg + neutralizing k + 100 mm excess kcl + 
 15 mm excess mgcl2(3)mg80 : rna + 11 crystallographic 
 mg + neutralizing k + 100 mm excess kcl + 
 80 mm excess mgcl2 mg0 : rna + 11 crystallographic 
 mg + neutralizing k + 100 mm excess kcl mg15 : rna 
 + 11 crystallographic 
 mg + neutralizing k + 100 mm excess kcl + 
 15 mm excess mgcl2 mg80 : rna + 11 crystallographic 
 mg + 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 .
crystallographic k and mg ions as well as water residues were retained 
 but the terminal gtp residue was removed ( as discussed previously ) . 

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 mg 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 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 .
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 .
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 ) .
initial structures for the umbrella sampling 
 windows were chosen differently in three sections along the ion - leaving 
 path as follows:(1)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.(2)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.(3)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  restraints , and the reaction coordinate 
 value was increased from 6.9 to 20.0  with a 100 kcal / mol  distance restraint .
all other rna heavy atoms were restrained 
 in position with 20 kcal / mol  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  .
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  restraints , and the reaction coordinate 
 value was increased from 6.9 to 20.0  with a 100 kcal / mol  distance restraint .
all other rna heavy atoms were restrained 
 in position with 20 kcal / mol  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  .
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 was used with collision frequency of 5 ps to hold the temperature at 298.15 k ; the shake algorithm was used to constrain the bonds involving hydrogen 
 atoms ; the direct space cutoff was set to 8.0  ; and particle 
 mesh ewald 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  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 mg 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 mg 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 mg , a mg located in bulk solvent
this way , the reaction coordinate 
 between 0.5 and 20.0  was sampled at 0.1  intervals for 
 a leaving mg at 0 , 15 , and 80 mm excess mg concentrations with 2 ns simulation time per umbrella ( us_mg0_mg , 
 us_mg15_mg , and us_mg80_mg simulation sets ) . for no excess mg 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 .
two additional umbrella sampling simulations were conducted 
 for 
 excess mg 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  mg mutated 
 system in the first scenario were changed to the parameters developed 
 by qvist as well as hydration free 
 energy ( hfe ) and iod parameters developed by li et al .
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 , 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 ) . 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 level and a 6 - 31g(d ) basis set .
the 
 calculated wave function was used to extract the electron density 
 and its topological properties with the aimall suite of programs . to calculate the energy difference between the 
 two ions ,
the same solvated system was calculated switching k for mg and the properties of the electron density 
 were extracted using the same methodology .
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 mg ( 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 mg ions , it is not possible to converge the mg 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  positional restraints 
 on all rna atoms ( competition_1 set ) .
figure 3 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 3 , 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 
 mg 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 
 mg , but the difference between the free energies of binding 
 of two ions is not high enough to completely prevent mg from binding .
distance of the closest k ( top ) and mg ( 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 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 3 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.73.5 
  from six oxygen atoms in the crystal structure , as illustrated in figure 2 .
figure 4 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 4 , 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 1 , 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 .
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 
 mg 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 mg , 
 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 can not 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 ,
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 5 ) , 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 . 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 1 ) the water molecules can leak into the pocket 
 and directly interact with k ( figure 5 ) .
however , the comparison between the two 
 histograms in figure 5 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 ) .
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 can not 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 mg , 15 mm mg , and 80 mm mg were added to independent md simulation 
 experiments , of which the 80 mm mg is similar to the 
 ion environment used in the crystallization of the 1qa6 and 1hc8 structures .
a series of md simulations were performed 
 during which the ion ( na , k , or mg ) 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 ion was restrained using a harmonic distance restraint 
 with 100 kcal / mol  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 mg with 100 mm k ; 5 mm mg with 100 mm na ; or 5 mm mg with 100 mm 
 nh4 . here
, we 
 replace k in 
 the crystallographic pocket with mg 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 mg ( in all scenarios ) and na ( in 
 the scenario with no excess mg ) . 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 .
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 1 ) .
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 mg , in which mg163 had already dissociated before the k ion - exchange event . in the resulting pmf plots ( figure 6 ) ,
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 .
pmf profiles 
 for pulling out k , mg , and 
 na from the gac rna monovalent ion - binding site .
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 mg to 
 the monovalent 
 binding pocket is always energetically less favorable than binding 
 of k. binding of mg is energetically less 
 favorable than binding of na in the case with no excess 
 mg .
these data support the hypothesis that the pocket 
 is specific for monovalent ions and agree 
 with the experimental observations that k stabilizes the 
 gac rna tertiary structure better than a bound na . in the presence of mg in mg163 position , it is 
 unfavorable 
 for mg 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 mg 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 mg in mg163 position and with no 
 excess mg 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 
 mg in the absence of mg in the mg163 position 
 decreases with increasing bulk mg concentration .
therefore , 
 mixed occupancy of k and mg in this pocket 
 in the whole ensemble seems possible especially at higher mg 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 . to test how the magnesium binding to this pocket might be affected 
 by different magnesium models , a magnesium - pulling experiment was 
 repeated .
mg parameters 
 included the traditional qvist parameters and hfe and iod parameters developed by li et al .
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 .
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 allnr  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 mg models 
 vary in magnitude but none exceed the binding free energy of k ions .
pmf profiles for pulling out mg from the 
 gac rna monovalent 
 ion - binding site with different magnesium models .
pmf plots are calculated 
 under the same conditions described in figure 6 . as a complementary approach to study the energetic preference of 
 either k or mg 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 ; 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 2 ) .
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 2 .
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 mg .
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 i d of residue ) .
electron 
 density redistribution of the binding pocket with either 
 k or mg 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 mg 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 .
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 mg , 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 
the 
 surplus charge is extracted from the three backbone phosphate oxygen 
 atoms present in the pocket .
a different picture occurs with mg , 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 .
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 3 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 mg show very few instances of mg binding .
however , mg 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 
 mg is not huge in the absence of mg163 .
this supports 
 the identification of potassium in this pocket as described by conn 
 et al . rather than magnesium as suggested 
 by wimberly et al .
also , preference of 
 k over na in this pocket in low mg concentration ion pull - out simulations agrees with the trend in 
 melting experiments observed by lu and draper and wang et al .
, where 
 k / mg stabilized the gac tertiary structure 
 more efficiently than other monovalent and divalent ions .
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 mg 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 .
further support for such conclusion 
 comes from the fact that the closest coordinating oxygen atom in 1mms ( op1 of a1073 ) is located 2.6  away from the chelated 
 mg ( table 1 ) , whereas the mean experimental value for mg 
o 
 distances in rna structures is 2.08  according to the cambridge 
 structural database and distances longer than 2.4 
  are suspected to be misidentified .
however , according to the results of the ion pull - out simulations 
 and the free - energy difference between k and mg for their binding to the pocket ( figure 6 ) , 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 .
the md simulations suggest that occupancies of the mg163 and 
 k58 
 positions by related ions are mutually exclusive at lower mg 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
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 mg ion
mobile water molecules could relax the pocket geometry , allowing higher 
 mobility of the nucleotides and also exchange of the bound k.
 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
Methods
Results
Conclusions