Hofmeister Series for Metal-Cation–RNA Interactions: The Interplay of Binding Affinity and Exchange Kinetics

A large variety of physicochemical properties involving RNA depends on the type of metal cation present in solution. In order to gain microscopic insight into the origin of these ion specific effects, we apply molecular dynamics simulations to describe the interactions of metal cations and RNA. For the three most common ion binding sites on RNA, we calculate the binding affinities and exchange rates of eight different mono- and divalent metal cations. Our results reveal that binding sites involving phosphate groups preferentially bind metal cations with high charge density (such as Mg2+) in inner-sphere conformations while binding sites involving N7 or O6 atoms preferentially bind cations with low charge density (such as K+). The binding affinity therefore follows a direct Hofmeister series at the backbone but is reversed at the nucleobases leading to a high selectivity of ion binding sites on RNA. In addition, the exchange rates for cation binding cover almost 5 orders of magnitude, leading to a vastly different time scale for the lifetimes of contact pairs. Taken together, the site-specific binding affinities and the specific lifetime of contact pairs provide the microscopic explanation of ion specific effects observed in a wide variety of macroscopic RNA properties. Finally, combining the results from atomistic simulations with extended Poisson–Boltzmann theory allows us to predict the distribution of metal cations around double-stranded RNA at finite concentrations and to reproduce the results of ion counting experiments with good accuracy.


■ INTRODUCTION
Metal cations are essential for the folding and function of nucleic acids. Given the negative charge of the sugar− phosphate backbone, positively charged ions are required to screen the electrostatic interactions. In addition, metal cation binding to specific sites stabilizes the three-dimensional structure further and assists or even enables catalytic reactions, for instance, in ribozymes.
Until recently, most scientific work has focused on Mg 2+ , as it represents the dominant cellular divalent cation. However, an increasing number of experimental results indicate that a large variety of nucleic acids properties depend not only on the salt concentration and valence of the ions but also on the ion type. 1−25 Examples include the stability and folding of large RNA ribozymes 1−4 or the melting temperature of canonical and noncanonical nucleic acid structures [17][18][19]25 (see Table 1 for further examples).
Three examples are highlighted in the following. The folded structure of large ribozymes is stabilized most efficiently by Mg 2+ ions while the stabilization efficiency decreases with decreasing charge density of the ions. 1−4 On the other hand, the folding kinetics is 20−40-fold faster if Mg 2+ ions are replaced by Ba 2+ or Na + ions, 14 indicating competing effects in thermodynamic stabilization and folding kinetics. Finally, the catalytic activity can be reduced significantly by replacing small amounts of Mg 2+ by Ca 2+ ions. 20,21 What is the origin of these ion specific effects? In general, such ion specific effects are ubiquitous and known for a large variety of physical properties such as osmotic coefficients, solubility of gases and colloids, protein precipitation, or the catalysis of chemical reactions. 26 In most cases, anions and cations can be ranked reproducibly in a Hofmeister series according to their influence on these macroscopic properties. 27 At the same time, the widespread applicability of the Hofmeister series suggests that ion specific phenomena have a common origin. By now we know that water structuring effects and direct ion−macromolecule interactions are equally important. 28−30 Therefore, in order to gain microscopic insights into the origin of ion specific effects in RNA systems, a quantitative description of cation− RNA interactions is required.
Over the past decade, many high-resolution structures of RNAs have been resolved, in some cases revealing bound cations in the crystallographic electron density map. Experimental techniques, especially high resolution nuclear magnetic resonance and single molecule Forster resonance energy transfer spectroscopy, provide increasing insight into the folding landscape and the dynamics of conformational changes. 31−34 However, resolving conformational changes and monitoring how cations are involved is challenging experimentally, since most cations are spectroscopically silent or difficult to detect. Moreover, ion exchange rates can be much faster than the experimental resolution, prohibiting a quantitative characterization. 35 As an alternative to spectroscopic methods, ion counting (IC) experiments can be used to gain insight into the ionic atmosphere around nucleic acids by measuring how many ions of each species are present in the atmosphere. 11,12,36−38 However, IC experiments do not allow one to resolve the exact interactions, since they cannot resolve the affinity of the metal cations toward different binding sites on RNA.
Here, molecular dynamics simulations (MD) can contribute important insights and provide a unique description of the dynamics including the behavior of the solvent and the ions. 33,39−41 Over the past years, MD simulations have been used to determine the binding affinities of metal cations to different binding sites. 42−49 Additionally, theoretical models and simulations have been used to investigate the effect of different cations on helix−helix interactions, 47 the distribution of metal ions around RNA motifs, 48 or competitive binding to tertiary structures. 49 However, a general theory that quantifies cation−RNA interactions and encompasses different binding sites on RNA is still missing due to the complex interactions between RNA, ions, and water. Theoretical predictions from Manning or Poisson−Boltzmann (PB) theory 50−52 neglect all nonelectrostatic contributions to ion−RNA interactions and are therefore inapplicable to predict ion specific behavior or to reproduce the results from ion counting experiments. Modifications to PB theory accounting for finite size effects or the tight bound ion model 52−54 provide improvement. However, they neglect the subtle interplay of ion hydration and partial ion dehydration upon binding which has to be taken into account by explicit water MD simulations if ion specificity is to be understood beyond the heuristic fitting of experimental data. 36,55 In this work, we fill this gap by using all-atom molecular dynamics simulations and enhanced sampling techniques to obtain atomistic insights into the interactions of metal cations and RNA. For the three most common ion binding sites on RNA, we calculate the binding affinity and exchange kinetics of eight different mono-and divalent metal cations. The two complementary observables display different ion specific features. The binding affinity of metal cations at the backbone follows a direct Hofmeister series but is reversed at the nucleobases. Direct and reversed affinities of different binding sites on RNA are the origin of the observed high specificity and selectivity of metal ion binding to RNA. 15,17−19 In addition, the exchange rates and consequently the lifetime of contact pairs between different metal cations and RNA cover almost 5 orders of magnitude, ranging from tens of picoseconds to tenths of microseconds. The interplay between the considerably different lifetimes and ion binding affinities, are the fundamental reason for the observed ion specificity in RNA stability and folding (Table 1). Finally, combining the results from our simulations of single ions with extended Poisson−

■ METHODS
Atomistic Model and Simulation Setup. In this work, we investigate the interactions between different mono-and divalent metal cations and an RNA dinucleotide. The RNA dinucleotide consists of two guanine nucleobases (G) and contains three different ion binding sites: 15 The backbone binding site consisting of two nonbridging oxygen atoms of the phosphodiester linkage (ions can bind either to atom O1P or O2P), the guanine N7 position and the O6 carbonyl group of guanine (see Figure 1). All force field parameters of the RNA dinucleotide were taken from Amber99sbildn* 56 with parmbsc0 57 and χ 0L3 58 corrections, and the TIP3P water model was used. 59 For the metal cations, Li + , Na + , K + , Cs + , Mg 2+ , Ca 2+ , Sr 2+ , and Ba 2+ , we used recently optimized force field parameters. 60 This choice was motivated by the fact that the optimized parameters yield accurate ion pairing properties as judged by comparison to experimental activity coefficients and accurate exchange kinetics as judged by experimental water exchange rates. In turn, these parameters allow us to capture the fine differences between the distinct metal cations and to quantify their binding affinity and exchange kinetics.
All simulations were performed using the Gromacs simulation package 61 version 5.1.4. Periodic boundary conditions were applied and the particle-mesh Ewald method 62 was used with tin foil boundary conditions, a Fourier spacing of 0.12 nm and a grid interpolation up to order 4 to handle long-range electrostatic forces and the electroneutrality condition. Bonds to hydrogen atoms were constrained by the LINCS algorithm 63 and a 2 fs time step was used. Close Coulomb real space interactions were cut off at 1.2 nm and Lennard-Jones (LJ) interactions were cut off after 1.2 nm, respectively. Long-range dispersion corrections for energy and pressure were applied to account for errors stemming from truncated LJ interactions. Initially, the RNA was placed in a cubic box (4.0 × 4.0 × 4.0 nm 3 ) filled with TIP3P water molecules. Harmonic restraints (1000 kJ mol −1 nm −2 ) were applied on the heavy atoms of the dinucleotide to prevent large conformational changes during the temperature replica exchange molecular dynamics (TREMD) simulations. Therefore, our results do not include the conformational adjustment of the RNA upon ion binding. Subsequently, we performed an energy minimization, equilibration, and production run for each simulation setup. After energy minimization using the steepest descent algorithm, we employed NVT and NPT simulations as pre-equilibration using the Berendsen scheme. For each 1 ns simulation, the temperature was kept constant with the Berendsen thermostat 64 at a value of 300 K and pressure was kept at 1 bar using the Berendsen barostat. 64 Production runs were performed in the NPT ensemble using the isotropic Parrinello−Rahman barostat, 65 and the temperature was maintained using the velocity rescaling thermostat with a stochastic term 66 with coupling constants 5.0 and 0.1 ps for pressure and temperature, respectively.
Umbrella Sampling. We used umbrella sampling simulations to calculate the potential of mean force (PMF) as a function of the distance between a single cation and the ion binding sites on the RNA. For the backbone, we use the distance between the cation and the O2P atom. Two umbrella force constants k 1 and k 2 were chosen. For all ions except Mg 2+ , k 1 was used for r < 0.4 nm and k 2 for larger distances with a window spacing of 0.01 nm (60 windows in total). For Mg 2+ , k 1 with window spacing 0.005 was used for 0.25 nm < r < 0.34 nm and k 2 with window spacing 0.01 otherwise (70 windows in total). The values for k 1 and k 2 for all cations and binding sites are shown in Table S3. Each umbrella simulation was performed for 10 ns discarding 1 ns for equilibration. The PMF was calculated using the weighted histogram analysis method (WHAM), 67 and a Jacobian correction was applied since distances are used rather than Cartesian coordinates where P(r) is the probability distribution obtained from WHAM. Note that P(r) is directly related to the free energy profile F(r) via F(r) = −k B T ln P(r). Errors were calculated from block averaging (see the Supporting Information).
Kinetic Rate Coefficients. The binding kinetics was determined from TREMD simulations in combination with a kinetic rate model. 68 The kinetic rate model extracts kinetic rate coefficients from replica exchange simulations by using the unbiased trajectory segments within a maximum likelihood framework. 68 We used 26 replicas of the dinucleotide system with a single cation placed randomly in the simulation box. A temperature range between T = 300 and 400 K was used with a series of temperatures obtained using the temperature generator. 69 To obtain multiple binding events, the simulations were performed for 100 ns for Na + , K + , and Cs + , 240 ns for, Li + , Sr 2+ , and Ba 2+ , and 400 ns for Ca 2+ . Exchanges between replicas were attempted every 0.5 ps.
To gain insight into the binding kinetics, we define an indicator function h(r) that defines whether the ion is bound or unbound where κ = 5.0 and γ = 1.05. b 1 and b 2 are the cutoffs for the bound and unbound state, respectively. They were defined for each cation such that only transitions between inner-sphere and outer-sphere conformation were counted (see Table S4). According to our definition, the calculated rate coefficients correspond to the exchange between the inner-sphere and outer-sphere conformation. Based on these rates, one can also derive an expression for the exchange rates between inner-sphere conformation and bulk (see the Supporting Information). All rates are listed in Table S6. The dynamics of the populations p i of the bound (b) and unbound state (ub) follows from the phenomenological rate equations for a two state system The kinetic rate coefficients, k + and k − , can be calculated from 68 where N is the total number of transition between the states. t ub/b is the total time that the ion spent in the unbound or bound state. Errors for the rate coefficients are calculated from transition counts. 68 To validate the approach, we performed additional straightforward MD simulations at 1 M salt concentration, 100 ns for LiCl, NaCl, KCl, and CsCl, and 500 ns for BaCl 2 , SrCl 2 , and CaCl 2 . The results from TREMD and MD are in excellent agreement (see Figure S2) in the Supporting Information. Binding Free Energies. We consider the binding process D + M ↔ DM of a dinucleotide (D) and a metal cation (M) in dilute solution. The equilibrium state obeys where [i] denotes the concentration of species i. Let p 0 and p 1 be the fraction of RNA without and with a cation bound, respectively. Expressing p 0 and p 1 in terms of integrals over the PMF, we can rewrite eq 6 where [D] tot is the total concentration of RNA in the system, [M] = 0.025 is the concentration of cations in the simulation box, r L is the radius of the sphere that contains the same number of water molecules as our simulations box (N water = 2140), and r † is the position of the maximum in the PMF (see Table S5 in the Supporting Information). Alternatively, the equilibrium constant can be calculated directly from the kinetic rate coefficients: To facilitate comparison with experiments, we compute a standard free energy of binding with the standard concentration c 0 = 1 M. Free Energy Perturbation. Since the exchange kinetics for Mg 2+ is very slow, eq 8 cannot be used to get a reliable estimate of the binding affinity. Therefore, we performed free energy perturbation simulations to validate our results using the Bennett acceptance ratio method. 70 The simulations were also performed for Ca 2+ and compared to our other methods. We calculated the binding free energy as the difference of decoupling the cation in bulk and in the ion binding site. In each setup, the pathway of transferring the ion from the gas phase to bulk water or to the binding site was split into two separate processes: First, a neutral Lennard-Jones particle was created, then the charge was increased in a second step. Along the solvation pathway, a transformation parameter λ was used with 40 evenly distributed replicas. Each replica was simulated for 1 ns discarding the first 200 ps for equilibration. To improve convergence, soft-core potentials are applied for LJ potentials with α = 0.3, linear λ scaling, and a radius power of 6. To ensure that the cation remained at the correct position along the solvation pathway even when the interaction between the cation and the rest of the system was decreased to zero, a harmonic position restraint was applied with force constant k = 10 6 kJ/(mol nm 2 ) for Mg 2+ and k = 10 4 kJ/(mol nm 2 ) for Ca 2+ . In a separate simulation, the contribution of releasing the position restraint was calculated using 27 uneven λ-values. 71 Note that the force constant for Mg 2+ is much higher than that for Ca 2+ . A larger value was required for Mg 2+ since the equilibrium position in the binding pocket is in the repulsive part of the Lennard-Jones interaction potential. Consequently, a lower value implicates that the cation leaves the correct position in the binding pocket upon decreasing the interaction leading to a defective result for the binding free energy. Further details can be found in the Supporting Information.
Extended Poisson−Boltzmann Theory. In the following, we model a dsRNA as a uniformly charged cylinder of radius R = 1.3 nm with a uniform charge density σ RNA . R was taken as the averaged RNA radius for dsRNA in the A-form. 72 We assume that the dsRNA consists of G−C pairs which contain only three different ion binding sites (one on the backbone and two on the guanine base). Since the ion binding sites are segregated, we assume that the cylinder consists of three types of surface patches corresponding to the three different binding sites (see Figure 2B). In the model, the ion can interact locally either with the backbone or site N7 or O6 through the respective PMF, V i bb , V i N7 , and V i O6 . The PB equation including the ion specific PMFs in cylindrical coordinates reads 28,73 where Φ(r) is the electrostatic potential, ϵ 0 is the dielectric constants of vacuum, ϵ is the relative dielectric constants of TIP3P water (ϵ = 82), 74 where c i0 is the bulk salt concentration of species i, ξ bb = 0.5, and ξ N7 = ξ O6 =0.25 since we have two backbone sites per guanine. V + bb , V + N7 , and V + O6 are the PMFs for the backbone and N7 and O6 on the nucleobase, respectively. In this work, we do not consider the effect of Langmuir pubs.acs.org/Langmuir Article co-ions, i.e., we set their PMFs to zero (V − bb = V − N7 = V − O6 = 0). However, anions may play an additional role. 75 Note that the PMF for the backbone results after subtracting the long-ranged Coulomb interaction in order to avoid double counting within PB theory (see the Supporting Information). To ensure that the ionic density does not exceed its physical limit, the proper ionic densities c i (r) are calculated from the unrestricted densities c̃i(r) using a Fermionic distribution 76,77 Here, we chose a i = 0.7 nm, corresponding to the effective distance of phosphate groups in a dsRNA. The PB equation was solved numerically on a one-dimensional grid with pm resolution yielding the electrostatic potential Φ(r) and the radial ionic densities c i (r). The electrical potential satisfies the bulk boundary condition Φ(r → ∞) = 0. In addition, we use the constant charge boundary condition, dΦ(r)/dr| r=R = −σ RNA /ϵϵ 0 at the surface located at r = R where σ RNA = −0.831e/nm 2 corresponds to the surface charge density (e is the fundamental charge) of negatively charged phosphate groups. The number excess of cations (or depleted anions) per nucleotide, commonly expressed as Γ i , is calculated from the ionic distributions where h R = 6.78 nm is the height of the RNA with a net charge of −46e to resemble the experimental setup, 11 N A is Avogadro's number, and V f = 10 24 nm 3 /L is a conversion factor. In experiments, the competition among ions is measured from titration of a second cation species in the presence of a constant buffer concentration. 11 where c is the concentration of the titrated salt, Γ 0 and Γ 1 are the number of associated background cations at the starting and end states, n is the Hill coefficient, and c 1/2 is the competition constant defined as the concentration of the competitive cation at which half of the number of the background cations are replaced. 11 n and c 1/2 are fitted simultaneously.

■ RESULTS AND DISCUSSION
The present study focuses on the interactions of mono-and divalent metal cations with the three most common ion binding sites on RNA. In the following sections, we first present the free energy profiles of metal cation binding and the ordering of cations according to their binding affinity. Subsequently, we discuss the kinetics of ion binding. Finally, we calculate the ionic distributions around dsRNA from extended Poisson−Boltzmann theory and compare to ioncounting experiments.
Interaction between Cation and RNA Binding Sites. Figure 3 shows the free energy profiles underlying the RNA− Langmuir pubs.acs.org/Langmuir Article cation interaction as a function of the distance between the ion and the three ion binding sites investigated in this work: The backbone binding site ( Figure 3A), the N7 binding site ( Figure  3B), and the O6 binding site ( Figure 3C). All free energy profiles show two metastable states that are separated by a free energy barrier. The first minimum corresponds to inner-sphere conformations, in which the cations form a direct contact and have a reduced water coordination number (Figure 3, snapshots (i)). The second minimum corresponds to outersphere conformations, in which the contact is mediated by an additional water molecule (Figure 3, snapshots (o)). The free energy profiles for the backbone show two opposing trends: The binding affinity of the cations, which is dominated by the deepest minimum in the free energy profile, increases with increasing charge density of the cations. Therefore, metal cations with high charge density such as Mg 2+ or Li + bind preferentially. At the same time, the free energy barrier increases with increasing ionic charge density since the ions become more strongly hydrated. The removal of one water molecule from the first hydration shell in order to facilitate a direct contact becomes therefore increasingly difficult. Consequently, the ion binding kinetics gets slower with increasing charge density of the cations.
A similar behavior is also evident at the N7 and the O6 binding site. Again, the barrier height increases with increasing charge density of the cations. However, the binding affinity at the N7 binding site is exactly opposite compared to the backbone. Here, metal cations with high charge density such as Mg 2+ or Li + are most strongly repelled and cations with low charge density such as K + or Cs + bind preferentially. In addition, the outer-sphere conformation becomes more stable than the inner-sphere conformation for the divalent metal cations and for Li + . The situation at the O6 binding site is similar for the monovalent ions. However, inner-sphere and outer-sphere conformation have similar stability for the divalent cations except for Mg 2+ .
The binding free energies ΔG b 0 of the different cations in inner-sphere and outersphere conformation are compared in Table 2. For inner-sphere binding at the backbone, the binding free energy follows a direct Hofmeister series: Cs + ≈ K + ≈ Na + > Ba 2+ > Li + > Sr 2+ > Ca 2+ > Mg 2+ (Figure 4A). At the N7 binding site, the ordering according to the binding free energy is exactly opposite and follows a reversed Hofmeister series: Mg 2+ > Ca 2+ > Sr 2+ > Li + > Ba 2+ > Na + > K + ≈ Cs + ( Figure  4B). At the O6 position, the situation is less clear, and the ordering follows a partially reversed Hofmeister series: Mg 2+ > Ca 2+ ≈ Sr 2+ ≈ Li + ≈ Ba 2+ ≈ Na + > K + ≈Cs + ( Figure 4C). Note that the reversal of binding affinities when going from the backbone binding site to the N7 binding site is in accordance with Collins' empirical like-seeks-like rule: 78 The oxygen atoms of the phosphate group have a high charge density (see Table  S1 in the Supporting Information). Consequently, cations with high charge density such as Mg 2+ and Li + have a higher binding affinity. By contrast, the binding sites on the nucleobases have a lower charge density and cations with a similarly low charge density show a higher binding affinity. However, the situation is more complex in reality due to the competition of inner-and outer-sphere conformation ( Table 2).
Our results are in qualitative agreement with experimental results for a mononucleotide: 15 In experiments, Mg 2+ ions show a higher binding affinity to the phosphate group compared to Ca 2+ ions. Moreover, at the nucleobase the simulations predict outer-sphere conformations for both cations with similar binding affinities in agreement with the experimental results. However, the simulations slightly overestimate the measured binding affinities by about 1−2k B T. Our results also explain the unusual behavior of G-quadruplexes. In G-quadruplexes, the ordering of the ions according to their stabilization efficiency is exactly opposite compared to canonical RNA structures (Table 1). Since the most important binding sites in G-quadruplexes are the guanine bases along the central pore, the stability is largely determined by the interactions of the cations with the O6 binding site. The reversed ordering of the cations according to their stabilization efficiency is a direct consequence of the reversed binding affinity at the nucleobases.
Taken together, our results show that the different binding sites are highly selective: Ion binding sites on the backbone preferentially bind metal cations with high charge density in inner-sphere conformations. By contrast, ion binding sites on the nucleobases preferentially bind cations with low charge density. Here, the probability of finding divalent metal cations in inner-sphere conformations is considerably lower. Moreover, at the N7 position, the inner-sphere conformations for divalent cations with high charge density is only transient. Divalent cations are therefore expected to be found solely in outer-sphere conformation in equilibrium. This result is in agreement with recent results by Leonarski et al. 79,80 and sheds further light on the controversial discussion regarding the  81 Cation Exchange Kinetics. In the following, we turn to the kinetics of metal cation binding to the RNA dinucleotide. The indicator function h, which defines whether an ion is bound (eq 2), allows us to gain insight into the binding kinetics. For Na + and Li + at the backbone binding site, h is shown in Figure 5A,B. For Na + , we observe frequent exchanges between the bound state and unbound state, while for Li + exchanges are much rarer. At the same time, ion binding coincides with the removal of one water molecule from the first hydration shell ( Figure 5C). Figure 5D-G shows the kinetic rate coefficient for the different metal cations at the backbone and N7 binding site. (The results for the O6 binding site are shown in Supporting Information Figure S2 and Table S6). Note that, with our current approach, we are not able to determine the exchange rate for Mg 2+ . For Mg 2+ , the exchange kinetics is on the hundred millisecond time scale 82 which is out of reach of our TREMD simulations.
At the backbone binding site, the association and dissociation rates decrease with increasing charge density of the ions (Figure 5D,E). Since ion association requires the removal of one water molecule from the first hydration shell, the binding kinetics is determined by the time scale of water exchange which shows the same trend. 60,83 On the other hand, ion dissociation requires the replacement of one phosphate oxygen by a water molecule which is associated with an even higher free energy barrier (see free energy profiles in Figure  3A). Since the binding affinity increases with increasing charge density of the ions, the energy barrier increases with increasing charge density. Consequently, cation exchange is fastest for ions with a low charge density (Cs + ) and slowest for ions with a high charge density (Mg 2+ ). Overall, the ordering of the ions according to their dissociation kinetics follows a direct Hofmeister series at the backbone binding site.
At the N7 binding site, the association rates decrease with increasing charge density of the ions ( Figure 5G) and are similar in magnitude as at the backbone site. This behavior is expected, since it reflects that cation binding is largely determined by the kinetics of water exchange and therefore  Langmuir pubs.acs.org/Langmuir Article independent of the chemical details of the binding site. By contrast, the dissociation rates increase with increasing charge density. The ordering of the ions according to their exchange kinetics is exactly opposite as at the backbone and the ordering corresponds to a reversed Hofmeister series. This reversal results from the reversed binding affinities and is reflected in the free energy profiles ( Figure 3B). Here, the energy barrier increases with decreasing charge density. In summary, our results show that ion dissociation rates cover almost 5 orders of magnitude. Consequently, the lifetime of contact pairs formed between metal cations and the RNA backbone strongly depends on the charge density of the cation. For example, a K + −RNA contact pair is stable on the order of 0.1 ns while a Ca 2+ −RNA contact pair is stable for about 0.1 μs. This considerably longer time scale explains why Ca 2+ ions are much more efficient than K + in stabilizing folded RNA structures. 1,5 On the other hand, the longer lifetime of contact pairs of metal cations with high charge density can hinder RNA folding since folding intermediates or misfolded states are stabilized much longer as reflected in folding kinetics experiments of RNA (folding time: Mg 2+ > Ca 2+ > Sr 2+ > Ba 2+ ). 13,14 Ion Atmosphere and Comparison to Ion Counting Experiments. In the following, we provide a quantitative description of the ion atmosphere around dsRNA by combining simulation results and PB theory. In general, the ion atmosphere consists of site-specific and diffusive ions. In our theoretical model, specifically bound ion are taken into account via the free energy profiles obtained from the all-atom molecular dynamics simulations while the diffusive ions are treated on a mean field level (see Methods). Therefore, the model allows us to calculate ion excess without any free parameters and to compare directly to the results from ion counting experiments on dsDNA. 11,12 The extended PB theory reproduces the experimental titration curves for mono-and divalent with good accuracy ( Figure 6A−F). Moreover, the competition constant for monovalent ions agrees quantitively with the results from experiments ( Figure 6G). The competition constants reveal preferential binding of Li + compared to Na + , K + , and Cs + . Based on the simulation results, the preferential binding of Li + results from the higher affinity at the phosphate group compared to the other ions. The binding affinities at the guanine base are much smaller and therefore not relevant ( Figure 4). We find that including the results from MD simulations via the free energy profiles provides an improved description of the ion atmosphere compared to approaches that rely solely on steric effects. For example, the recent threedimensional reference interaction site model (3D-RISM) predicts that Na + outcompetes Cs + by a factor of 1.8 84 while no significant preference is predicted from our model in agreement with experimental results. 12 For the divalent metal cations, the agreement of the competition constants from simulations and experiments is not as satisfactory as that for monovalent cations ( Figure 6H). The extended Poisson−Boltzmann model captures the overall trend of the competition constants indicating that the specific binding affinity of the cations decreases according to Mg 2+ > Ca 2+ > Sr 2+ > Ba 2+ as expected from the simulations. However, the model overrates the competition constants since the calculated binding affinity with the current Mg 2+ force field is too high compared to experimental results 15 as discussed previously. Additionally, the subtle differences between dsDNA and dsRNA in ion counting, recently reported, 38 cannot be resolved by our model (data not shown).
The model presented here takes the diffusive ions and the site-specific ions into account while the conformational changes of the RNA are neglected. In larger and complex RNA structures, the ionic distribution is affected by conformational changes of the RNA. For instance, the addition of cations decreases the overall extension due to increased electrostatic screening. At the same time, the structural changes of the RNA affect the ionic atmosphere and can alter the binding of the site-specific ions. Here, predictive models 52,85 and all-atom simulations 39,46,48 could be used to gain further insights.

■ CONCLUSION
Predicting how metal cations affect the physicochemical properties of RNA is challenging due to the entangled contributions of backbone, nucleobases, ions, and water. In order to quantitatively describe ion−RNA interactions and to gain microscopic insight, we performed MD simulations for the three most common binding sites on RNA. The results clearly reveal which metal ion preferentially binds to the phosphate or nucleobase binding site and whether the interaction involves an inner-or outer-sphere contact ( Table  2). In addition, the calculated free energy profiles and the exchange rates include the contributions of ion hydration and direct ion−RNA interactions and are the key to understand ion specific behavior. In particular, the higher binding affinity of metal cations with high charge density to binding sites involving phosphate groups and the much longer lifetime of contact pairs rationalize the higher stabilization efficiency and longer folding times of large ribozymes in the presence of those ions compared to ions with low charge density. Moreover, the reversed behavior of RNA quadruplexes (highest stabilization efficiency of ions with low charge density) is a direct consequence of the reversed binding affinity at the nucleobases predicted from our simulations. So far, we considered only the interactions of metal cations with three individual binding sites. These are idealized cases while RNA macromolecules contain complex binding sites in which the metal cations are coordinated by several functional groups. However, experimental results have shown that the binding affinity in RNA systems is largely additive, 15 opening up the possibility for the modeling of complex RNA binding sites in the future. Furthermore, ion binding may affect RNA flexibility and provoke local structural changes. This implies that a complete quantitative description would require to include both the RNA degrees of freedom and the ion specific interactions between the metal cations and the RNA binding sites.
In all cases, Mg 2+ was found to play a distinct role due to its high binding affinity and slow exchange kinetics. In fact, the exchange kinetics for Mg 2+ is so rare that the rate could not be determined even with the acceleration of TREMD. Here, more work is required for a reliable calculation.
Finally, we have combined the results from the MD simulations with extended Poisson−Boltzmann theory in order to predict the distribution of metal cations around dsRNA at finite salt concentrations. The advantage of this approach is that it treats the specifically bound ions realistically including the changes in the hydration shell upon binding via the free energy profiles. This allows us to quantitatively describe the ion atmosphere and to calculate ion competition constants for cations in close agreement with experimental results. In the future, the results for the divalent metal cations could be improved by optimizing the force field parameters of Mg 2+ based on experimental binding affinities 42 or by including ion−ion correlations that might become important at high salt concentrations.
Force field parameters, umbrella force constants, convergence of the free energy profiles, definition of bound and unbound states, inner-sphere and outersphere free energy, kinetic rate coefficients, free energy perturbation and removal of the Coulomb part from the potential of mean force (PDF)