Concentration-dependent and configuration-dependent interactions of monovalent ions with an RNA tetraloop

Naturally occurring polyelectrolytes, such as nucleic acids, accumulate high concentrations of monovalent ion species as a means of screening the strong electronegative charge of their sugar-phosphate backbones and allowing the molecule to adopt secondary structures or form backbone-backbone interactions.^1–3 Even the simplest ribonucleic acid (RNA) molecules have been observed to adopt an ensemble of configurations⁴ and both experimental and computational studies have shown that a complex solution of cosolutes can have a profound effect on the relative populations of that ensemble.^5,6 Characterization of RNA-solution interactions is thus essential to provide a comprehensive description of stable RNA configurations, especially when considering the complex mixture of ions, cosolutes, and crowding agents present in intracellular environments where RNA molecules fold.⁷ 

Small RNA molecules have important biological functions in gene expression and small-molecule-sensing;^8,9 therefore, understanding the interactions of RNA oligomers with ions and solvent is important. However, the distribution of ions around small RNA oligomers is not well known, and the Manning ion condensation theory,¹⁰ which works well for large polyelectrolytes, does not apply to smaller, more spherically shaped, oligomers. Molecular dynamics (MD) becomes an ideal tool for studying these systems by describing the ions and solvent explicitly, and the thermodynamic formalism used to describe large polyelectrolytes can be employed to interpret the results of MD simulations and reveal the effects of salt concentration and oligomer configuration.

The distributions of concentrated salt solutions around nucleic acid structures have been an active field of study for over forty years, and the effect of salt concentration and polyelectrolyte structure on these distributions has been closely investigated.¹⁰ The Gibbs-Duhem equation [Eq. (1)] forms the theoretical basis of these phenomena and is exemplified by dialysis. In this experiment, a chamber containing a charged molecule, such as a nucleic acid with its neutralizing cation (i.e., K⁺:RNA⁻), is separated from a concentrated salt solution (KCL) by a semipermeable membrane. The concentrated salt solution (referred to as the “bulk solution”) will distribute its constituent species (K⁺, CL⁻) into the solvation layers around the nucleic acid via the electrostatic potential between the ions and the polyanionic nucleic acid backbone. The distribution of ions between the two chambers will, at equilibrium, reflect a balance between the electrostatic energy and the osmotic pressure of the salt,

By this relation, the free energy of a nucleic-acid:salt system is equal to the sum of the products of the population and chemical potential of each salt species (N_i and μ_i, respectively) for M different salt species (i), given a system entropy (S), temperature (T), volume (V), and pressure (P). The populations (N_i) for different salt species are important parameters in this equation since in equilibrium they will distribute across the membrane and establish the Gibbs-Donnan equilibrium, balancing osmotic pressure and electrostatic charge.

Multiple studies by Anderson and Record have addressed the thermodynamics of how these cosolute species distribute around polymers.^11–15 The number of ions for a given salt species that accumulate near a nucleic acid relative to the bulk solution in dialysis is termed preferential interaction coefficients (Γ^PIC). By accumulation (or exclusion) of a certain number of a given salt species (⁠$|ΓiPIC|$⁠), the bulk chemical potential is kept constant and the Donnan equilibrium is maintained. The value of Γ^PIC can be determined experimentally for three-component solutions [(m₁) water, (m₂) nucleic acid, and (m₃) salt] by measuring the change in the chemical potential of the salt as a function of the presence of the nucleic acid [μ₂₃, Eq. (2)].¹⁶ The chemical potentials for the nucleic acid and ions are constant when this is achieved,

Based on a two-domain perspective (local/bulk) and assuming a negligible concentration of nucleic acid, the preferential interaction coefficient for each salt species can be solved by taking the observed number of ions (N_i) at a given distance (r) from the nucleic acid and subtracting the number of ions that would be expected in a bulk sample (⁠$NiNW$⁠)_bulk given the number of water molecules (N_W) within the same distance (r). Both of these numbers can be determined by averaging over the course of a canonical (NVT) MD simulation and taking pairwise distances between the nucleic acid and the surrounding ions and water [Eq. (3)],

This ion and water counting method has been used to a great effect in describing the preferential interactions of denaturing cosolutes around RNA oligomers and protein systems based on all-atom molecular dynamics simulations^5,17 and has been used to calculate the distributions of different-sized monovalent cations for large RNA kissing loops.¹⁸ 

While the calculation of preferential interactions is useful as an ion-counting technique, it must be acknowledged that $ΓiPIC$ is agnostic to the charge state of either the nucleic acid or the ion. If the system is electrically neutral, the local distribution of the ionic cloud around a nucleic acid can be described by a similar charged-particle-counting technique that identifies populations of excess ionic charge around the nucleic acid.^19,20 These populations are referred to as Donnan coefficients (⁠$Γ±,RNAD$⁠) and are often compared to an “ideal” Donnan limit that is half the charge of the nucleic acid (0.5⋅$$Z$$⁠). These values are based on equations for polyelectrolyte-ion systems determined by Record and Anderson.¹³ Provided that the counterion and co-ion populations are nearly equivalent, and the nucleic acid concentration is vanishingly small, Eqs. (4) and (5) are sufficiently accurate for measuring the Donnan coefficients with the same two-domain perspective used for Γ^PIC,

Recent atomic emission spectroscopy (AES) experiments have indicated that cationic preferential interaction coefficients (⁠$Γ+PIC$⁠) for large DNA and RNA molecules are linearly dependent on salt concentration up to 1m NaCL.²¹ The authors of this work hypothesized that additional studies using all-atom molecular dynamics were necessary in order to verify the same trends at high salt concentrations.²² For this purpose, we chose to investigate the concentration and configuration-dependent distributions of a monovalent salt (KCL) around a well-characterized nucleic acid system: an eight-nucleotide RNA oligomer (sequence: GCGCAAGC). This sequence has been observed to form Watson-Crick base-pairs between the 5′ and 3′ terminal ends and generate a tetraloop (GCAA) even in the absence of divalent cations [Fig. 1(a)].^23,24 Previous all-atom molecular dynamics simulations of this oligomer have shown spontaneous formation of both canonical right-handed (A-RNA) base-pairs and non-canonical left-handed (Z-RNA) base-pairs.⁴ Both structures form intramolecular hydrogen bonds with half of all Watson-Crick edges and consequently limit the exposure of these RNA moieties to the solution. In some cases, these base-pairs are not formed, the RNA molecular surface is more exposed to the solution, and the molecule is described as “unfolded.” All three RNA conformational states can form and interconvert in unbiased replica-exchange simulations and were observed to have the most even distribution between different configurations in intermediate simulation temperatures (e.g., 350 K).

In this study, we characterize the radial distributions, preferential interaction coefficients, and Donnan coefficients of monovalent salt (0.1 to 3m KCL) for this model nucleic acid using all-atom molecular dynamics. We highlight two major facets of solution distributions around RNA: (I) the concentration-dependence of preferential interactions for water and salt species and (II) the effect of RNA configuration on ion accumulation and exclusion. The results of our study show that molecular dynamics can reproduce the experimental values and trends observed in concentration-dependent ion distribution experiments and that these trends extend from sub-molal concentrations akin to cells to high salt concentrations similar to solutions for crystallographic preparation.

Our study employs molecular dynamics simulations of three RNA configurations with different base-pairing interactions (Fig. 1) and nine KCL salt concentrations from 0.1 to 3 $molkg$ (m).

The three selected RNA configurations, and the range of KCL concentrations employed, are intended to provide a comprehensive picture of the behavior of concentrated KCL solutions around some of the most basic RNA structures at a computationally tractable scale. The primary structural difference between the three chosen configurations is the geometry of the 5′ and 3′ terminal bases that form CG Watson-Crick base-pairs and constitute the hairpin stem region. In the native A-RNA configuration [Fig. 1(a)], these bases form a canonical right-handed helix, with all stem nucleobases in an anti-χ configuration and all ribose sugars forming a C3′-endo pucker. The non-canonical Z-RNA configuration [Fig. 1(b)] forms a left-handed helix, the guanosine nucleobases adopt a syn-χ configuration, while the ribose sugars of the smaller cytosine nucleobases adopt a C2′-endo pucker. The only criteria for the unfolded configurations [Fig. 1(c)] are that they do not form an ordered helical structure.

In Table I, concentrations are given in molality (m) in order to satisfy the equations for preferential interaction coefficients (Γ^PIC) and Donnan coefficients (Γ^D) as first described by Record.¹³ 

Salt concentrations in this work are selected in order to span a range of salt concentrations from mammalian cell cytoplasm (0.1m) to the limits of KCL solubility (∼3m). The sizes of our simulation boxes are chosen by taking the maximum end-to-end distance of this RNA oligomer in unfolded configurations from previous simulations (L_U ∼ 40 Å) and making cubic boxes with edges greater than the sum of this length and twice the predicted Debye screening length for each monovalent salt concentration: L > (L_U + 2·δ). This way, any ions at the RNA oligomer surface are screened from the ions at the RNA oligomer surface in the periodic image. Debye screening lengths are determined for monovalent 1-1 electrolyte solutions based on the equation $δ=ϵkbT2⋅e2Na[M]12$⁠, where ϵ = 78 for water, k_b is Boltzmann’s constant, T = 350 K, e is the monovalent charge, N_a is Avogadro’s number, and [M] is the concentration in molarity. To provide appropriate distances, and to include at least 30 ions of each ion species in the solution, sub-molal salt concentrations are described with three times as many water molecules as systems with ≥1m KCL.

Our computational models use modified AMBER99 force field parameters for RNA²⁵ and a fixed-charge model for monovalent ions.²⁶ Both force fields were originally modified to fit numerous experimental measurements for RNA and monovalent ions in the context of the TIP3P water model, which we employ for the solvent.²⁷ Each system is equilibrated in a 10-ns NPT simulation at 350 K, 0.1 MPa using a velocity-rescaling thermostat²⁸ and the Berendsen barostat²⁹ with coupling constants of 1.0 ps and 0.2 ps, respectively. Seed frames for production runs are drawn from snapshots where the temperature is 350 ± 2 K and the pressure is 1 ± 2 bar. Production simulations are run using NVT ensembles at 350 K for 100 ns each, and snapshots of these trajectories are saved every 2 ps. The RNA configurations are not restrained during the simulation, and any time blocks involving transitions between different configurational states are discarded from the analysis. This was necessary, as several folding events occur in unfolded-initialized simulations at all KCL concentrations. Blocks including 100-ns of simulation are used from the unfolded trajectories, while in folded-initialized simulations no transitions occurred and all configurations are included in the analysis. Each system is run in quadruplicate to insure consistency of the results, making a total of 400 ns for each system and thus a total of 10.8 μs of the total simulation time. All simulations are run with the GROMACS 4.5.5 simulation package.³⁰ An important note about the kinetics of these simulations is the lack of long time scale ion coordination. Individual ions are rarely coordinated to particular atoms on the RNA surface for any length of time exceeding 10 ps, and thus all calculations shown in this work represent general distributions of each ion species around different RNA molecules rather than any tightly coordinating interactions.

Radial distribution functions are analyzed using the g_rdf program in GROMACS 4.5.5. Preferential interaction coefficients (⁠$Γ±,RNAPIC$⁠) and Donnan coefficients (⁠$Γ±,RNAD$⁠) are determined based on Eqs. (3) and (4) with an in-house code written using the MOSCITO package.³¹ Linear fits for the Donnan coefficients are determined as a function of concentration over two ranges (0.1–1m KCL) and (1–3m KCL) using GNUPLOT version 5.0 and are plotted with the same program.³² All concentration-dependent and configuration-dependent effects on ion distributions are observed in the first few solvation layers around the RNA oligomers, while at increased distances these distributions become bulk-like [g_±(r) → 1, $Γ±,RNAPIC→$0, $Γ+,RNAD→Z$⁠, and $Γ−,RNAD→$0]. For this reason, only the first solvation layers are plotted and described in the results.

Radial distributions around polar (O,N) and non-polar (C) RNA atoms provide a general description of the relative distributions of each species [$gCL−(r)$ , $gK+(r)$ , and $gOW(r)$ , “RNA” omitted from the subscript for brevity] and permit analysis of the effect of folded (A-RNA, Z-RNA) versus unfolded RNA geometries. A representative spread of the concentration-dependence is given by 0.1m, 1m, 2m, and 3m KCL for all three configurations and all three solution species [Fig. 2 (CL⁻); Fig. 3 (K⁺); Fig. 4 (O_W)]. All other concentrations are plotted in the supplementary material. The same relative displacement is observed for all solution species: accumulations around polar atoms occur as near as 2.5 Å, and accumulations around non-polar atoms occur as near as 3 Å.

Accumulation of CL⁻ around polar RNA atoms is considerably greater than around non-polar atoms [Fig. 2(a)]; however, these values are well below unity at all concentrations, indicating overall exclusion of CL⁻. Increasing KCL concentration 30-fold only increases $gCL−$(r) roughly two-fold for both polar and non-polar RNA atoms. At all concentrations, the highest $gCL−$(r) values (at r = 3.25 Å) are observed in unfolded configurations (0.2–0.35), while Z-RNA values are the lowest [$gCL−$(r) ≤ 0.1], while A-RNA values consistently intermediate, regardless of KCL concentration. It should be noted that for greater distances (r > 3.25 Å), A-RNA and Z-RNA $gCL−$(r) values are quite similar, and both are weaker than unfolded configurations. These observations indicate that both folded configurations show greater CL⁻ exclusion than unfolded configurations.

In contrast to the interactions of CL⁻ with RNA, K⁺ interactions show tighter interactions near both polar and non-polar RNA atom types in all configurations and all KCL concentrations. Polar atoms show non-zero $gK+$(r) as close as 2.5 Å and peaks near 2.75 Å, while non-polar atoms begin to accumulate K⁺ at 3 Å with a broader, flatter peak than the comparable non-polar peak for $gCL−$(r). At concentrations below 0.5m KCL (Fig. 3 and the supplementary material), Z-RNA shows the highest $gK+$(r) values for both polar and non-polar RNA atoms. At concentrations >0.5m KCL, the $gK+$(r) value at the 2.75 Å peak for the unfolded configurations begins to exceed that of Z-RNA, though $gK+$(r) values at distances greater than 2.75 Å around Z-RNA are nearly equal to the values of unfolded configurations up until 2m KCL, while A-RNA is weaker than both. The greater accumulation of K⁺ around Z RNA, relative to A RNA, is consistent with observations of Z RNA being more stable at high salt concentrations.³³ 

These observations indicate that the ion distributions in the first solvation shells are strongly affected by the RNA configuration and that unfolded structures accumulate at least trace amounts of all ion species regardless of the KCL concentration. Even at the highest concentrations, in this study, the unfolded configurations show $gK+$(r) > 1 at the RNA surface, indicating that the salt is more concentrated around the unfolded RNA than in the bulk solution even when the salt is near the concentration at which it might begin to precipitate in the experiment (∼3m).

An important consideration in these simulations is the behavior of water as a function of changing the KCL concentration and box size. As stated in the methods, the sub-molal (0.1–0.75m KCL) simulations have nearly three times as many water molecules as the systems with ≥1m KCL concentrations in order to effectively screen between the periodic boundaries as the Debye lengths increase.

The $gOW$(r) distributions are minimally affected by the KCL concentration despite an order of the magnitude change in the number ions. The primary differences among the different simulations come from the configuration of the RNA molecule (Fig. 4 and the supplementary material). For A-RNA and unfolded configurations, the polar atoms show nearly equivalent values for $gOW$ at the two peaks located at 2.75 Å and 5 Å, while Z-RNA is consistently weaker at both positions. There is a minor exclusion of water from polar RNA atoms in the region between these two peaks and in this space; A-RNA is more similar to Z-RNA, indicating a slightly greater average density of water around unfolded configurations.

One other notable feature in $gK+$(r) (Fig. 3) and $gOW$(r) (Fig. 4) distributions is the relative strength of the secondary polar RNA atom peak near r ∼ 5 Å. For K⁺; this peak maximum is at 4.75 Å and is weaker than the more proximate peak at r ∼ 2.75 Å, while for water the first and second peaks are approximately equivalent.

From these observations, we see that K⁺ is accumulating within the first solvation layers of the RNA surface, regardless of the configuration, and predominantly within the first layer. Within the same region, CL⁻ ions also accumulate but to a much lower degree. These two ions show opposite dependencies on concentration: K⁺ accumulation weakens with increasing KCL concentration while CL⁻ accumulation increases. Thus the ratio of counterion to co-ion density (⁠$gK+$ : $gCL−$⁠) is greatest at 0.1m KCL and weakest at 3m KCL.

The g(r) values can address the relative densities of different solution species in the solvation layers around RNA, but in order to account for the numerical excess (or deficit) of different ion species around each RNA configuration, an ion-counting measure is necessary. This is where preferential interaction coefficients (⁠$Γ±,RNAPIC$⁠) become necessary to enumerate the ions distributed at the RNA surface.

To quantitate how K⁺ and CL⁻ interact with RNA, and how these ions distribute into the solvation layers around RNA, the accumulated excess of each ion within 10 Å (approximately two solvation layers) is described using preferential interaction coefficients (⁠$Γ±,RNAPIC$⁠), as described in Eq. (3).

The number of CL⁻ ions excluded from RNA (⁠$Γ−,RNAPIC$⁠, Fig. 5 and the supplementary material) increases with concentration. By 3m KCL, the number of excluded CL⁻ ions is the same as the charge of the RNA oligomer. Both folded RNA configurations consistently exclude most CL⁻ ions, though all three configurations show notably similar trends with different offsets beyond the first 3 Å from the RNA surface. For concentrations >1m KCL, the value of $Γ−,RNAPIC$ converges by approximately 6 Å from the RNA surface, and approximately 2 CL⁻ ions are excluded for every additional 1m increase in the KCL concentration. A small accumulation of CL⁻ occurs at r ∼ 2.5 Å around unfolded configurations and shifts the overall trend of $Γ−,RNAPIC$ for the unfolded configurations relative to the folded configurations. In sub-molal KCL concentrations, the offset is barely perceptible, and similar CL⁻ exclusion occurs for all configurations. However, by 3m KCL, this offset accounts for one less CL⁻ ion excluded from the unfolded configurations relative to the number of CL⁻ ions excluded from both folded configurations.

To more thoroughly characterize preferential interactions of the counterions around different RNA configurations, we determine the preferential interaction coefficients for K⁺ (⁠$Γ+,RNAPIC$(r), Fig. 6, and the supplementary material). At sub-molal KCL concentrations, all RNA configurations show a greater degree of preferential accumulation of K⁺ than CL⁻ exclusion within 10 Å of the RNA surface. However, for concentrations ≥1m KCL, the trend reverses (⁠$|Γ−,RNAPIC|>|Γ+,RNAPIC|$⁠), and K⁺ accumulation is the weaker of the two preferential interaction coefficients at the RNA surface. Thus it is more the exclusion of co-ions than the accumulation of counterions that defines the most proximate RNA solvation layers in increasing monovalent salt concentrations.

Interestingly, $Γ+,RNAPIC$ reaches a maximum (2.75–3) at 10 Å for all three configurations at sub-molal concentrations, but this maximum begins to decrease for ≥1m KCL concentrations (Fig. 6 and the supplementary material). For all KCL concentrations, the value of $Γ+,RNAPIC$(r) increases by approximately 1.5 over the distances 4.0 ≤ r ≤ 5.5 Å, indicating that on average nearly 2 K⁺ ions accumulate at the outer edge of the first solvation layer around all configurations of this RNA molecule across the entire range of KCL concentrations.

Atomic emission spectrometry by Jacobson et al. determined that the preferential cation coefficients (⁠$Γ+PIC$⁠) for large linear RNA and DNA molecules will decrease linearly as a function of salt concentration in sub-molar NaCL concentrations (0.069–1.02M) and are predicted to reach zero near 1.4M.²¹ We see an agreement with our own results: a monotonic decrease in $Γ+,RNAPIC$ for all configurations with increasing salt concentration and $Γ+,RNAPIC≤0$ within the first solvation layer for all salt concentrations ≥1.5m KCL.

As with the anionic preferential interaction coefficients (⁠$Γ−,RNAPIC$⁠), the most critical ranges that define structural-specificity of the cationic preferential interaction coefficients (⁠$Γ+,RNAPIC$⁠) are within 3 Å of the RNA surface. A-RNA configurations only accumulate K⁺ within the first 5 Å up to 0.75m KCL and begin approaching the preferential exclusion around 1m KCL [$Γ+,RNAPIC$(<4 Å) ∼ 0]. This indicates that the surface of A-RNA attracts roughly the same number of K⁺ ions as would be found in an equivalent volume of 1m KCL bulk solution. At higher concentrations, the A-RNA shows increasing K⁺ exclusion.

The Z-RNA configurations, in agreement with $gK+$(r), have the highest $Γ+,RNAPIC$ values below 0.5m KCL, and at increasing KCL concentrations the values are intermediate between unfolded and A-RNA configurations, though within 5 Å of the RNA surface they are most frequently within the error of unfolded configurations. Both the Z-RNA and unfolded configurations consistently show a small increase in $Γ+,RNAPIC$ at 3 Å, indicating that on average at least one K⁺ ion preferentially interacts right at the RNA surface of these configurations.

These observations contrast with the exclusion of CL⁻ ions, which show Z-RNA to be entirely within the same error and following the same trend as A-RNA. From these findings, we determine that the non-canonical Z-RNA accumulates more K⁺ ions than the native A-RNA, even though both folded configurations exclude the same number of CL⁻ ions.

The great utility of $Γ±,RNAPIC$ is thus its ability to discern the raw counts of ions right at the RNA surface and reveal how differences in RNA structures affect the accumulation or exclusion of different ion species. The most telling interactions are occurring right at the RNA surface but also extend into the first few RNA solvation layers. In this work, all folded configurations exclude the most co-ions, and non-canonical Z-RNA and unfolded configurations accumulate the most counterions.

The last important piece of this puzzle is the distribution of the ionic charge and how the concentration of KCL and the configurations of RNA affect this distribution via the Donnan coefficients (⁠$Γ±,RNAD$⁠).

In order to account for the charge distributions of monovalent salt around these RNA models, we plot the distant-dependent Donnan coefficients for K⁺ [$Γ+,RNAD$(r)] in Fig. 7. Based on the equations for Donnan coefficients [Eqs. (4) and (5)], and the electroneutrality principle employed in determining these equations, the Donnan coefficients for CL⁻ (⁠$Γ−,RNAD$⁠) show the same trends as those for K⁺ (⁠$Γ+,RNAD$⁠) but with an offset of −7 based on the charge of the RNA molecule. These plots are shown in the supplementary material.

Both ion species in these simulations are monovalent, so the scale of K⁺ accumulation and CL⁻ exclusion—as calculated by the preferential interaction coefficients—is directly related to the charge distribution calculated by the Donnan coefficients. In sub-molal KCL concentrations, the absolute values for K⁺ preferential interactions exceed those of CL⁻ (⁠$|Γ+,RNAPIC|>|Γ−,RNAPIC|$⁠), and thus the number of accumulated K⁺ ions has a greater effect on $Γ+,RNAD$⁠. This fact is reflected by both the non-canonical Z-RNA and unfolded configurations showing the highest Donnan coefficients.

As KCL concentration increases, the number of excluded CL⁻ ions at the RNA surface outstrips the number of accumulated K⁺ ions, and in molal KCL concentrations (1–3m) CL⁻ exclusion becomes the dominant effect on charge distribution around these small RNA oligomers. The greatest CL⁻ exclusion is observed in A-RNA and Z-RNA configurations, and as a result the $Γ+,RNAD$ values for both folded configurations is less than the $Γ+,RNAD$ values of unfolded configurations over the same distance (r) from the RNA surface.

As with the preferential interaction coefficient, the Donnan coefficient shows its greatest configuration-dependence within 3 Å of the RNA surface—the first solvation layer. Only a small increase in $Γ+,RNAD$ is observed for both folded configurations, and this value weakens considerably with increasing KCL concentration. For A-RNA, the small increase disappears by 0.75m KCL. In all subsequent KCL concentrations, Z-RNA configurations remain slightly higher than A-RNA due to a greater accumulation of K⁺ ions within this first solvation layer. For KCL concentrations ≥0.5m KCL, unfolded configurations maintain the greatest ion accumulation in the first solvation layer and consequently show the highest Donnan coefficients.

To describe the concentration-dependence of $Γ+,RNAD$ near the RNA surface, we take the average value of $Γ+,RNAD$ from 4.0 to 6.5 Å for each configuration and plot as a function of the KCL concentration (Fig. 8). We see that the $Γ+,RNAD$ values at 1m KCL for A-RNA (3.001), Z-RNA (3.116), and unfolded configuration (3.323) all show a close agreement with the extrapolated per-nucleotide Donnan coefficient for calf-thymus DNA at a similar concentration (0.98M NaBr, $Γ+,7−ntD$ = 3.234).

Increasing the KCL concentration from 0.1 to 1m produces a near-linear reduction in the Donnan coefficient for all RNA configurations. An additional increase from 1m to 3m KCL produces an even stronger reduction that still approximates a linear concentration dependence. Based on these trends, line fits can be generated that describe the concentration-dependence of Donnan coefficients for each RNA configuration (Fig. 8). Fitting each configuration to a single trendline for all concentrations produces much larger least-square errors than fitting to separate lines for sub-molal (0.1–1m) and molal (1–3m) KCL concentrations, so two lines are described for each configuration. The data are interpreted in Table II.

From these data, we see that the solvation layers around folded tetraloop configurations become net-neutral (⁠$Γ+,RNAD=0$⁠) by approximately 2.5 or 3m KCL. By contrast, the solvation layers around unfolded configurations retain positive values at these same high KCL concentrations. Interestingly, our calculated Γ^D values for this small RNA oligomer at 1m KCL compare well with the Γ^D values for large DNA at 0.98M NaBr.¹ 

We conclude from these results that our simulations can describe a concentration-dependence in the Donnan coefficients for different configurations of this small RNA oligomer and that our calculated values agree with values derived from experiments on a much larger nucleic acid.

After parsing out the concentration and configuration dependencies in our results, we find that even small RNA configurations describe multiple equilibrium phenomena for proximate water and ions. From basic radial distribution calculations, we observed that polar RNA atoms (which represent backbone phosphate oxygens, sugar hydroxyls, and nucleobase carbonyls and amines) form the strongest interactions at the closest distances for all solution species. Non-polar RNA atoms (represented by all carbon species in sugars and nucleobases) accumulate solution species farther from the RNA surface and do not form interactions as sharp as those around polar atoms. Accumulations of CL⁻ increase with concentration and favor unfolded configurations, while K⁺ accumulations decrease with the KCL concentration and are predominantly observed around Z-RNA in sub-molal concentrations and unfolded configurations in higher (≥1m KCL) concentrations.

Using the preferential interaction coefficients (⁠$Γ±,RNAPIC$⁠) as ion-counting measures and the Donnan coefficients (⁠$Γ±,RNAD$⁠) to assess charge-distributions, we have further elaborated on the concentration-dependent and configuration-dependent distributions of these RNA simulations and arrived at several important conclusions regarding this small RNA system in the context of high monovalent salt concentrations.

First, the configuration-dependence for preferential interactions depends on the most proximate RNA-ion interactions within the first solvation layer (r ≤ 3.0 Å). Within this narrow range, the unfolded configurations show both K⁺ and CL⁻ accumulations at all KCL concentrations, non-canonical Z-RNA configurations also show K⁺ accumulations but minimal changes in CL⁻ distributions, and canonical A-RNA only weakly accumulates either K⁺ or CL⁻. At distances greater than 3.0 Å from the RNA oligomer surface, the preferential interaction coefficients are almost wholly dependent on the KCL concentration and show similar behavior for all three configurations. The two folded configurations, which form intra-molecular hydrogen-bonds that limit the number of polar moieties in solution that would attract or repel salt species, show similar levels of CL⁻ exclusion, but disparate levels of K⁺ accumulation: non-canonical Z-RNA accumulating more than A-RNA in all KCL concentrations, and more than unfolded configurations for KCL concentrations <0.5m. Due to the ion accumulations in the first solvation layer, unfolded configurations consistently maintain the highest preferential interaction coefficients for both ion species at increasing distances from the RNA surface in all KCL concentrations ≥0.5m.

Second, the relative distributions of salt species within the solvation layers around RNA oligomers are correlated with the concentration of KCL. In sub-molal KCL concentrations, the accumulation of K⁺ exceeds the exclusion of CL⁻ (⁠$|Γ−,RNAPIC|<|Γ−,RNAPIC|$⁠), while for concentrations ≥1m KCL, the exclusion of CL⁻ outstrips the accumulation of K⁺ (⁠$|Γ−,RNAPIC|>|Γ+,RNAPIC|$⁠). Increased CL⁻ exclusion has the added effect of reducing the Donnan coefficients to values less than the ideal Donnan limit (⁠$Γ+,RNAD<0.5⋅Z$⁠). Discontinuities are observed in the concentration-dependence of Donnan coefficients at 1m KCL, allowing two separate trends to be described for $Γ+,RNAD$ as a function of KCL concentration: a sub-molal (0.1–1m) trend where most values remain near the ideal Donnan coefficient and a molal (1–3m) trend where the coefficients drop rapidly with increasing KCL concentration.

Last, our simulations match experiments both qualitatively and quantitatively in terms of cation preferential interaction coefficients (⁠$Γ+,RNAPIC$⁠) and Donnan coefficients (⁠$Γ+,RNAD$⁠) and allow us to describe the relative charge density of our model RNA system. In agreement with the work of Jacobson and Saleh, we see a monotonic decrease in $Γ+,RNAPIC$ as a function of salt concentration and predict near zero excess K⁺ accumulation in the first solvation layer of all configurations at ∼1.5m KCL.²¹ The extrapolated Donnan coefficients of calf-thymus DNA, as measured by Strauss et al. in 0.98M NaBr, are in close agreement with the Donnan coefficients calculated for our RNA oligomer near similar concentrations of KCL.¹ Based on the linear reduction of $Γ+,RNAD$ in increasing KCL concentrations, the folded configurations show $Γ+,RNAD∼$ 0 at 3m KCL, while unfolded configurations are ∼0.5 at the same concentration. From these observations, we conclude that this small RNA oligomer, when folded into a canonical or non-canonical tetraloop, does not disrupt the charge distribution of a 3m KCL solution. Further analyses will be needed to confirm these findings.

Ribonucleic acids are deeply affected by the nature of their surrounding solutions. The change in salt concentration has been observed to affect ion condensation on large polymers^10,34 and induce transitions in helical states from native to non-native configurations in high salt concentrations similar to those used in this work.^35,36 The precipitation of RNA contaminants through increasing the salt concentration is a technique that has seen extensive use in the purification of larger nucleic acid assemblies including plasmids.³⁷ 

A great deal of work in recent years has attempted to provide more accurate descriptions of the ion atmospheres around RNA through experiments and calculations.^{19–21,38–40} Significant developments in computer modeling have aimed at describing the ion atmosphere and the Donnan equilibrium coefficients around large DNA and RNA molecules using mean-field Poisson-Boltzmann and Debye-Hückel theories.^39,40 Nevertheless, ion correlation effects are not included in mean-field approximations, and at high salt concentrations, these effects become relevant for obtaining an accurate model of the equilibrium distributions of the ion atmosphere. The correlation effects are included in explicit simulations and have been used to model monovalent ion distributions around RNA kissing loops.¹⁸ This and many mean-field studies considered monovalent salt concentrations below 1m and did not venture into higher concentrations where experiments suggest that the charge density of the solution would be on par with that of the nucleic acid itself.⁴¹ 

Multiple efforts have employed the use of MD simulations to quantify the behavior of water and ions around RNA with the goal of recapitulating the observations of experiments and to explain the behavior of nucleic acids at the all-atom level.^42–45 For small RNA oligomers, which are not accurately described by theories developed for large, linear polyelectrolytes, this is especially important.^34,46 Even our model, which details many important phenomena for small RNA oligomers and high salt concentrations, is itself limited by the water model (TIP3P), which does not include ion polarization effects nor does it account for other hydration effects that might be captured in more advanced water models.^47,48 Additional improvements to these models may include the use of biomolecular force fields that also employ polarization and more sophisticated water models.⁴⁹ 

From the results of this work, we have shown that concentrated monovalent salt solutions and varying RNA configurations have distinct effects on ion and water distributions at the surface of a small RNA oligomer, and these effects can be modeled with high performance molecular dynamics and all-atom force fields. These findings lend new insights into the utility of molecular dynamics for describing equilibrium solution behavior around small nucleic acids and the ability to compare and contrast these results with experimental data for larger nucleic acid molecules.

See supplementary material for plots for the radial distribution functions, preferential interaction coefficients, and Donnan coefficients for all three RNA configurations at several intermediate KCL concentrations (0.25m, 0.5m, 0.75m, 1.5m, and 2.5m).

The authors would like to thank Dr. Alan A. Chen for the development of the force fields used in this study. The Stampede Supercomputer at Texas Advanced Computing Center (TACC) and the computer clusters at the Center for Nonlinear Studies provided the data for this study. Funding for this work was provided by the NSF through Grant No. MCB-1050966 and by US DOE LDRD funds from Los Alamos National Laboratory.