Predicting Monovalent Ion Correlation Effects in Nucleic Acids.

Ion correlation and fluctuation can play a potentially significant role in metal ion-nucleic acid interactions. Previous studies have focused on the effects for multivalent cations. However, the correlation and fluctuation effects can be important also for monovalent cations around the nucleic acid surface. Here, we report a model, gMCTBI, that can explicitly treat discrete distributions of both monovalent and multivalent cations and can account for the correlation and fluctuation effects for the cations in the solution. The gMCTBI model enables investigation of the global ion binding properties as well as the detailed discrete distributions of the bound ions. Accounting for the ion correlation effect for monovalent ions can lead to more accurate predictions, especially in a mixed monovalent and multivalent salt solution, for the number and location of the bound ions. Furthermore, although the monovalent ion-mediated correlation does not show a significant effect on the number of bound ions, the correlation may enhance the accumulation of monovalent ions near the nucleic acid surface and hence affect the ion distribution. The study further reveals novel ion correlation-induced effects in the competition between the different cations around nucleic acids.

Introduction

Nucleic acids (DNAs and RNAs) are highly charged polyanions. The folding of nucleic acids involves metal ion binding, which causes charge neutralization and screening of Coulomb repulsion between backbone charges.[1] Indeed, electrostatic interactions make significant contributions to nucleic acid folding stability[2−6] and conformational equilibrium and transitions.[7] Moreover, metal ions are essential for many biochemical reactions such as catalytic reactions of ribozymes.[8−13] To understand the nucleic acid structure and function, a computational model for metal ion effects is highly needed.

In general, metal ions around the nucleic acids can be classified into site-specific bound (SSB) ions[14] and nonspecific bound (NSB) ions.[15] SSB ions are partially or fully dehydrated and trapped (chelated) at specific sites or regions such as pockets in three-dimensional (3D) structures,[16] whereas NSB ions are often hydrated and form an “ion cloud” (or “ion atmosphere”)[17−19] around the nucleic acid. SSB ions can interact directly with the nucleic acid[20] and participate in nucleic acid biochemical reactions.[8−12] Physical properties of SSB ions, including ion–nucleic acid[21] and ion–ion interactions,[22,23] can be probed by X-ray and NMR measurements. Theoretically, the binding sites of SSB ions can be predicted by a variety of methods such as the knowledge-based MetalionRNA model.[24]

The interactions between (the large number of) NSB ions and the nucleic acid provide a significant stabilizing force for nucleic acid folding. To understand the ion effects on nucleic acid stability, it is highly desirable to understand the number and the distribution of NSB ions. A particularly important property of NSB ions is the number of excess NSB ions (Γ) in the ion atmosphere. The number of excess NSB ions Γ is the number of excess ions above the number expected based on the bulk concentration. In experiments, Γ can be determined using various “ion count” methods, such as buffer equilibration and atomic emission spectroscopy,[17,19] anomalous small-angle X-ray scattering,[25] and dye indicator.[26−29] However, experimental determination of further details about the NSB ion atmosphere such as the ion distribution and ion-induced electrostatic free energy changes remains a challenge.[30] We need theoretical models to predict ion–nucleic acid interactions and ion binding properties.

Many theories and computational models such as counterion condensation (CC) theory[31−33] and nonlinear Poisson–Boltzmann (NLPB) theory[34−40] have been developed for ion–nucleic acid interactions. These theories and models have led to important insights into ion effects on nucleic acid structures and stabilities. However, most of these theories rely on the mean-field approximation by neglecting the ion coupling (correlation) effect. In an ionic solution, ions can be correlated due to volume exclusion and Coulombic interactions.[30] Such a correlation can be particularly strong for ions distributed in the close vicinity of the nucleic acid surface, where ion concentrations can be high. Previous experimental and theoretical results have suggested that the correlation effect could be important for multivalent ions.[17,30,41,42]

According to the ion correlation strength, the NSB ions can be further classified into two types: the diffusely bound (DB) ions for the weakly correlated ions[15] and the tightly bound (TB) ions for the strongly correlated ions.[42] Although both the DB and the TB ions belong to the NSB ions, they have different binding properties. The DB ions move around the nucleic acid diffusely in the region separated from the nucleic acid surface by a distance,[43] whereas the TB ions are clustered around the nucleic acid surface to form a thin layer of high local concentration.[42] The strong ion correlation lowers the mobility of the TB ions.[30] It is important to note that unlike SSB ions, TB ions are assumed to be fully hydrated and do not form direct chelation with the nucleic acid. Experiment–NLPB comparisons[17,19] showed that the NLPB model may provide reliable predictions for the weakly correlated ions (DB ions) but could also underestimate the ion binding effects for the strongly correlated ions (TB ions).

DB ions can be treated as a continuous fluid with the mean-field theory. However, TB ions, which are strongly correlated, need to be treated at the level of discrete many-particle distributions to account for the coupling between the different ions. Molecular dynamics (MD) simulations are able to fully account for the correlation effect in ion–nucleic acid interactions,[44] such as Mg2+ binding to the SAM-I riboswitch,[45] 5S rRNA Loop E Motif,[46] and ribozyme.[47] Monte Carlo (MC) simulations can also sample the correlated Mg2+ distributions around the nucleic acids.[48−50] The sampling of the discrete ion distributions intrinsically reveals the ion fluctuation effect, which can be important for RNA folding.[30,42]

The previous emphasis on the ion correlation effect has been nearly exclusively focused on the multivalent ions such as Mg2+ ions. Successful attempts have been made to address the correlation for multivalent ions by modifying the NLPB[51−53] or the CC models.[54] However, recent experimental evidence indicated that monovalent ions may also cause ion correlation effects in ion–nucleic acid interactions.[19,55,56] MD simulations can treat correlation effects with discrete monovalent ions;[57−59] however, the applications (for monovalent ions or divalent ions) are limited by the long computational time with the explicit treatment for salt, solvent, and atoms in nucleic acids.[60] To enhance the computational efficiency, the three-dimensional reference interaction site model (3D-RISM) was developed to implicitly treat the solvent in simulations by solving the Ornstein and Zernike integral equation.[61,62] Ions are treated as discrete particles in 3D-RISM. Compared to the NLPB model, the monovalent ion distribution predicted by 3D-RISM is in good agreement with the MD simulation results.[61] The 3D-RISM model partly relies on MD simulations, so the prediction results can depend on the selected force fields.[61]

The tightly bound ion (TBI) model is a hybrid statistical mechanical model for ion–nucleic acid interactions, which combines the explicit enumeration/sampling of the discrete ion distributions for strongly correlated ions and a mean-field approach for the continuous distributions of weakly correlated ions. The model was originally developed to predict the effect of multivalent ions, particularly, divalent ions, by considering explicitly the ion correlation and fluctuation effects.[42] In this model, the discrete many-body distributions for the TB (strongly correlated) ions are enumerated and the (correlated) ion–nucleic acid interactions are calculated for each discrete ion distribution. For the DB (weakly correlated) ions, the mean-field NLPB is applied to compute the electrostatic free energy. Extensive comparisons with experimental data have suggested that the TBI model provides improved predictions for the ion effects in nucleic acid folding stabilities.[63−71] The original TBI model is severely limited by the low computational efficiency due to the time-consuming enumeration method, even for coarse-grained (CG) ion distributions. Recently, to enhance the computation efficiency, MCTBI,[72,73] a new TBI model, was developed using the Monte Carlo method to sample the TB ion distribution. Compared to the original TBI model, the MCTBI model has two notable advantages: a significant (1000-fold) increase in the computational efficiency and the use of full three-dimensional coordinate-based ion distributions.

To account for the correlation and fluctuation effects for the monovalent ions, on the basis of the original MCTBI model,[72] we here develop a generalized MCTBI-like (gMCTBI) model in this paper. In the original MCTBI model, only multivalent cations such as Mg2+ are classified into the TB and DB ions, whereas all of the monovalent cations such as Na+ are regarded as DB ions.[72,73] Therefore, TB ions in the original model include only multivalent cations. In the current gMCTBI model, however, TB ions include not only the (strongly correlated) multivalent ions but also the (strongly correlated) monovalent ions. As a result, the gMCTBI model accounts for three types of ion–ion correlations: monovalent–monovalent, monovalent–multivalent, and multivalent–multivalent. The model has the potential to provide not only improved predictions for the number of bound ions but also more complete and detailed insights into ion binding properties such as the most probable ion distributions. Using gMCTBI, we investigate the ion binding properties around the nucleic acids for Na+ and Na+ mixed with other monovalent cations (such as Li+ and Cs+) or divalent cations (Mg2+) and validate the theoretical predictions with available experimental data. The model predictions for various systems also reveal new features of ion distributions and ion–nucleic acid interactions.

Results and Discussion

To compare our theoretical predictions with the experimental results, we compute the excess number of bound ions Γα for type-α ions; see eqs S4 and S5 in the Supporting Information (SI).[74] To show the results for the detailed distribution of the bound ions, we also compute the linear density ρΓ(x)[45] as the spatial distribution function[61,62] for the excess ions (see eq S6). The physical meaning of the linear density is that ρΓ(x)dx equals the excess number Γ of ions in the region between distances x and x + dx. Here, for DNA helix, x denotes the distance measured from the helix axis. For tRNAPhe, which has a more complicated structure, we define x as the distance to the RNA surface.[45] The explicit treatment for discrete ion distributions enables us to predict the average and the most probable distributions for TB ions (see eq S7). Furthermore, we compare the results with and without the monovalent ion correlation to assess its effect in ion binding. In this section, we mainly study the ion binding properties around the nucleic acids in various salt solutions: NaBr (Section ), mixed monovalent salt solution (Section ), and monovalent–divalent mixed salt solution (Section ). In Section , we will predict the detailed ion binding properties to show the importance of the monovalent ion correlation effects.

Binding Properties of Na+ around B-DNA

We first investigate the ion binding properties for a 24-bp B-DNA helix in a NaBr solution. Figure A shows the results for the excess (depletion) number of Na+ (Br–) ions predicted from the original MCTBI and the current gMCTBI models. Both results agree with the experimental data.[19] The negative charges on the DNA attract the cation Na+ and exclude the anion Br–. According to the linear density profile predicted by the gMCTBI model (Figure B), the cation (Na+) accumulation can be classified into the strong accumulation region (several high peaks in the region x < 20 Å) and the weak accumulation region (a long tail at x > 20 Å). The DNA helix in a high-[Na+] solution can attract more excess Na+ ions (higher peaks) in the strong accumulation region, resulting in a stronger charge neutralization and consequently less excess Na+ ions (a lower long tail) in the weak accumulation region. The decay of the linear density of the excess cations in the weak accumulation region can cause a lower total number ΓNa of the excess Na+ ions. In fact, the long tail of the excess Na+ in the dilute [Na+] can extend to more than 100 Å into solution. As a result, ΓNa decreases with the increase in [Na+] (see Figure A).

Figure 1

Predictions for ion binding for a 24-bp DNA helix in a NaBr solution. (A) Excess number of Na+ and Br– in various bulk [Na+]. Here, the experimental data (blue circles) is extracted from ref (19). (B) Linear density of the bound Na+ and Br– ions predicted by the gMCTBI model for the different [Na+]. (C) Predicted TB Na+-ion average distributions corresponding to the three peaks in (B). (D) Most probable distributions of the TB Na+ ions (purple balls) at [Na+] = 10, 100, and 200 mM. To show the ion distributions, we use VMD and Chimera to show the (C) average and the (D) most probable distributions, respectively.

Similar to the cations, the depletion of Br– ions around the DNA can also be classified into the strong depletion (a deep valley) and weak depletion (a long tail) regions (Figure B). For a higher salt bulk concentration [NaBr], more anions Br– would be excluded from the DNA helix due to the repulsion from the negatively charged backbone, leading to a valley in the linear density profile in the strong depletion region. In the meantime, the stronger charge neutralization (due to more excess Na+ and depleted Br– in the strong accumulation/depletion regions) results in weaker depletion effect for Br– in the weak depletion region (i.e., a higher long tail in the weak depletion region). In Figure B, we find that the difference in the long tails for Br– between high and low salt concentrations (e.g., between [Na+] = 200 and 10 mM) is smaller than that for Na+. Therefore, different from ΓNa, whose trend with [Na+] is mainly governed by the long tail in the linear density profile, the overall depletion (ΓBr) of Br– is dominated by the number of the depleted Br– in the strong depletion region (deep valley), giving rise to the decrease of ΓBr with an increasing [Na+].

As [Na+] changes, the net excess charge Γtot = ∑αZαΓα = ΓNa + |ΓBr| (α = ion species) remains near +46e. Considering that the DNA backbone charge is equal to −46e, the result here suggests that the DNA with the bound ions remains neutral. The MCTBI and the gMCTBI (with and without monovalent ion correlation, respectively) provide similar predictions for the excess and depletion number of ions. The result seems to suggest that the monovalent ion correlation effect might not be strong enough to cause significant differences in the total number of bound/repelled ions. However, previous studies based on the comparisons between MD simulation and NLPB model predictions indicated that the monovalent ion correlation effect could lead to different ion distributions near the surface of the nucleic acids[61] (Section ).

Because the gMCTBI model samples the discrete positions for each TB ion, particularly, the monovalent ion, which is treated with a continuum model in the original MCTBI model, gMCTBI can account for discrete distributions for both the monovalent and multivalent ions. As a result, one of the major improvements in the gMCTBI model is its ability to compute the fluctuations and ensemble (probability) distribution of the monovalent (as well the multivalent) ions [probability p(k) of finding a TB Na+ ion at site k]. Describing ion binding at the level of discrete ion positions for monovalent ions can uncover many ion binding properties unattainable through the original MCTBI model.

The Na+ ions in the strong accumulation region may dominate the nucleic acid–ion interactions because these ions show lower diffusion and closer interaction with the DNA.[43] As shown in Figure B, there are three density peaks in the strong accumulation region. In fact, the three peaks in the linear density profile are mainly from the TB Na+ ions (see Figure S1 in SI). For peak 1 (Figure C), the closest peak to the DNA helix, the (fully hydrated) Na+ ions are mainly distributed inside (and along) the deep, wide major groove of the DNA (see Figure C). For peak 2, in contrast, ions are mainly distributed around the minor groove. These hydrated ions are too bulky to enter the (narrow) minor groove. Physically, ions at peaks 1 and 2 are attracted to the major and minor groove regions due to the high DNA backbone charge densities in these regions. We note that the predicted ion probability distribution also shows scattered high-probability ion binding sites near the major groove. For regions (such as peak 3) further away from the DNA, the distinction between the groove and nongroove regions disappears and TB Na+ ions are more uniformly distributed around the DNA. Here, we should note that since the gMCTBI model does not account for the desolvation effect of the ions, the binding sites of dehydrated Na+ ions and the corresponding density peak are not predicted.

Another important advantage of the gMCTBI model (vs the original MCTBI) is that it can give the most probable (discrete) distribution of the TB monovalent ions (see Figure D). At a low [Na+] (e.g., 10 mM in Figure D), the most probable distribution corresponds to TB Na+ ions distributed along the minor groove, presumably due to the high charge density in the (narrow) minor groove region. At a higher [Na+], such as 100 mM, more Na+ ions are populated and interact with the minor groove. Additionally, at high [Na+] (200 mM in Figure D), in addition to the majority of ions, which are distributed in the minor groove region, there are Na+ ions (circles in the figure) distributed in the major groove region in the most probable distribution. Here, the most probable distribution of TB ions is determined by the many-ion interaction energy of all of the bound ions (including TB and DB ions; see Section 4 in SI). Hence, a minor movement of the TB ions from the sites in the most probable distribution may not cause a dramatic increase in the interaction energy, suggesting that the TB ions may fluctuate around the bound sites of the most probable distribution.

Interplay between Different Monovalent Cations

The gMCTBI model allows us to further investigate the ion binding properties for a nucleic acid in a mixed salt solution, which contains background cations (bc) Na+ and competing cations (cc) such as Li+, Cs+, or TMA+. Figure A–C shows the excess bound ion number around a 24-bp DNA helix in various salt conditions. Experiment[55]–theory comparisons indicate that both the current gMCTBI model and the original MCTBI model can accurately predict the number of excess ions. The model predicts that as the bulk concentration of the competing cations [cc] is increased, more competing cations would bind to the DNA due to the lower entropic cost for ion binding. In the meantime, less background cations (Na+) are bound to the DNA. The slightly reduced depletion of Br– (ΓBr) is accompanied by a slight increase in the total number of excess cations (Γcc + Γbc). The predicted results are similar to those for a single salt solution (see Figure A).

Figure 2

(A–C) Excess number of ions predicted by the gMCTBI (red lines with squares) and the original MCTBI (black lines with triangles) model for a mixed solution with (A) Na+ and Li+, (B) Na+ and Cs+, and (C) Na+ and TMA+. The experiment data come from ref (55). (D–F) Linear density of the ions solved by the gMCTBI model corresponding to the cases in (A)–(C), respectively.

As shown in Figure A–C, for the different salt solutions [Li+]eq ≈ 35 mM (versus 40 mM Na+), [Cs+]eq ≈ 59 mM (versus 50 mM Na+), and [TMA+]eq ≈ 67 mM (versus 50 mM Na+), respectively, the competing and the background cations show a similar number of excess ions for the three solutions. Here, [cc]eq denotes the competing cation concentration under which the competing cation and background cation reach an equal number of excess ions: Γcc = Γbc. Although the different solutions show a similar number of excess ions, smaller cations are more competitive and show higher peaks in the close vicinity of the DNA (and lower peaks and tails in regions away from the DNA; see Figure D–F). The result indicates that the smaller cations are bound to the DNA more tightly. The predictions agree with previous results based on 3D-RISM[62] and MD simulations.[58]

Interplay between Monovalent and Divalent Ions

As shown in Figure A–C for the 24-bp DNA duplex and Figure D for a tRNA in various salt conditions, gMCTBI gives overall improved predictions than the original MCTBI. Furthermore, the results reveal an intriguing interplay between the monovalent cation, the divalent cation, and the anion. For a solution containing background cations and competing cations, an increase in the competing cation concentration would lead to more competing cations and less background cations bound to the nucleic acid.

Figure 3

Ion excess number predicted by the gMCTBI (red lines with squares) and the original (black lines with triangles) MCTBI models for (A) 6 mM background Mg2+ ion mixed with Na+, (B) 5 mM background Mg2+ mixed with Na+, (C) 20 mM background Na+ mixed with Mg2+ for a 24-bp DNA, and (D) 10 mM background Na+ mixed with Mg2+ for a 76-nt tRNA (PDB id: 1TRA(81)). Here, the experimental data are adopted from ref (55) for (A), ref (17) for (B) and (C), and ref (78) for (D).

For a mixed monovalent–divalent salt solution, as shown in Figure A,B (with Mg2+ as the background cation), for a 24-bp DNA helix, the number of excess background cation Mg2+ and competing cation Na+ reaches equilibrium at 54 and 49 mM Na+ in the background of 6 and 5 mM Mg2+, respectively. The total excess charge Γtot keeps the nucleic acid neutralized. For solutions with Na+ as the background cation, as shown in Figure C (for the 24-bp DNA) and Figure D (for a tRNA), the equilibrium occurs at 1.5 and 0.39 mM Mg+ (competing cation) with background cation concentrations of 20 and 10 mM Na+, respectively. The results show quantitatively the monovalent–divalent ion competition.

For a solution with Mg2+ as the background cation (Figure A,B), the depletion number of the anions decreases with the increase in [Na+] of the competing cation, whereas in the cases with Na+ as the background (Figure C,D), the depletion number of the anions increases with [Mg2+]. The above trend of ΓCl is a result of the decrease/increase in the number of bound Mg2+ ions, which dominate over Na+ ions in RNA binding. As more Mg2+ ions are bound to the nucleic acid, less anions will be excluded in the region far away from the nucleic acid, resulting in an increase in ΓCl.

Moreover, even in the solution mixed with Mg2+ ions, the entire nucleic acid-solution system still keeps charge neutrality (see Figure ). In fact, in a small space around a phosphate, we might observe the charge inversion (the change of net charge sign) or even a giant charge inversion[75] (net charge large than +e) if one or two TB Mg2+ ions are bound to a phosphate. However, in other regions beyond the small space above, phosphates in the nucleic acid and the TB cation-attracted anions can neutralize the inverse net charge in the local space around a phosphate. We also note that the charge inversion phenomenon of the entire nucleic acid-solution system can indeed be observed for trivalent or higher-valency cations.[74,76]

Effect of the Monovalent Ion Correlation

For a solution of (single or mixed) monovalent cations, neglecting the monovalent ion correlation effect may still provide good estimations for the ion excess number (see Figures A and 2A–C). However, the explicit sampling of discrete monovalent ion distributions in the gMCTBI model, which accounts for ion correlation effects, can provide more detailed information such as ion spatial distributions (see Figure C,D). In fact, other ion binding properties beyond the excess number of ions suggest that the ion correlation effect is indeed important for the monovalent ions around the nucleic acid even if the solution contains only monovalent ions. For example, the linear density profiles in Figure A shows that monovalent ion–monovalent ion correlation can promote Na+ accumulation near the DNA surface. This is because ion correlation lowers the Coulomb energy through the cooperative organization of the ions and thus induces more ions to bind to the nucleic acid. The finding agrees with the predictions from 3D-RISM and MD simulations,[61] which also suggested that there are more cations near the DNA surface than that predicted by mean-field theory.

Figure 4

(A) Comparison of linear density for Na+ around the 24-bp DNA helix between the predictions from the original MCTBI (black) and the gMCTBI model. (B) Most probable distributions of the background ions and the competing ions: 40 mM Na+ versus 35 mM Li+, 50 mM Na+ versus 59 mM Cs+, and 50 mM Na+ versus 67 mM TMA+. Red, yellow, purple, and green balls represent the Li+, Cs+, Na+, and TMA+ ions, respectively.

The correlation effect for monovalent ions can enhance the ion–ion competition near the DNA surface. Figure B shows the most probable distributions of the TB ions in the presence of competing cations at concentration [cc]eq. In the case of Li+–Na+ competition, Li+ ions are found to be distributed only along the minor groove of the DNA helix. However, for the bulkier competing cations (Cs+ and TMA+), more Na+ ions are found to bind to the DNA (six Na+ ions versus three Cs+ ions and seven Na+ ions versus three TMA+ ions, respectively). Therefore, even though the DNA attracts the same total (excess) number of cations, the smaller cations win the competition (against the larger cation) and are more likely to be found in the close vicinity of the DNA surface. However, as a caveat, we note that ion dehydration is not considered here. A smaller ion may bind more tightly with the DNA, but, on the other hand, it is less likely to dehydrate as water molecules wrap around smaller ions more tightly.[58] In the future, we need to develop a model that can integrate hydration with the correlation effect for bound monovalent cations.

For a Na+–Mg2+ mixed solution, compared with the original MCTBI model, as shown in Figure , the gMCTBI model can give overall improved predictions for the ion excess number in the ion atmosphere. The result suggests the importance to consider the Na+–Mg2+ ion correlation. In fact, previous studies have suggested the importance of such a correlation effect in other macromolecular systems such as protein solutions.[77] For the 24-bp DNA helix system, the original MCTBI model, which ignores the Na+–Mg2+ ion correlation, overestimates the excess bound Na+ ions. The correlation effect is more pronounced for systems involving a larger number of monovalent bound cations, such as solutions with high [Na+] in a fixed Mg2+ background (Figure A,B) or with low [Mg2+] in a fixed Na+ background (Figure C). For the tRNA, Figure D shows that the original MCTBI model underestimates the excess number of Na+ ions due to neglecting the Na+–Na+ correlation effect.

As shown by the linear density in Figure for the 24-bp DNA helix at high (200 mM) Na+ with 5 mM background Mg2+, the original MCTBI overestimates Na+ binding and underestimates Mg+ binding near x = 10 and 15 Å. The result is due to the missing correlation effect for Na+ ions in the original MCTBI model. For Na+ ion binding, although as shown in Figure A, the Na+–Na+ correlation can enhance Na+ binding (through multi-ion cooperative organization), the Na+–Mg2+ correlation (mainly repulsive through discrete charge–charge volume exclusion and Coulombic repulsion) causes a reduction in Na+ binding (see Figure A). Therefore, the original MCTBI model, which ignores the effect of discrete Na+ ions and the Na+–Mg2+ correlation in the TB region, overestimates Na+ ion binding. Furthermore, as shown in Figure B, the reduced Na+ ion binding results in an increased Mg2+ binding. However, in a solution of high [Mg2+] with background Na+ or a dilute [Na+] solution with background Mg2+, ion binding is dominated by Mg2+ and the Na+-involved correlation effect is weaker than the Mg2+–Mg2+ correlation. Therefore, the gMCTBI model and the original MCTBI model, which differ only in the treatment of the monovalent ion-induced correlation effect, give similar predictions for the excess number of bound ions.

Figure 5

(A, B) Linear density of the bound (A) Na+ and (B) Mg2+ ions around a 24-bp DNA helix in the solution with 200 mM Na+ and 5 mM Mg2+. The inset in (B) shows the accumulated Mg2+ ions around the DNA helix corresponding to the two peaks. (C) Most probable distributions of the TB ions in the solution containing 5 mM background Mg2+ ions mixed with [Na+] = 30, 100, and 200 mM Na+. The red and green balls denote the Mg2+ and Na+ ions, respectively

The density plot for Mg2+ (Figure B) shows two peaks. Our predictions for the most probable ion distributions (Figure C) indicate that the two peaks correspond to bound cations in the major and minor grooves, respectively. For low [Na+] (30 mM), the gMCTBI model predicts that only Mg2+ ions are closely distributed around the DNA. As [Na+] increases, less Mg2+ ions are found near the DNA surface. Only at high [Na+], such as 200 mM (>[Na+]eq), Na+ ions are found in the close vicinity of the DNA surface in the most probable ion distribution.

For the tRNA case, as shown in Figure D, for a limited range in [Mg2+] < 1 mM, the gMCTBI model predicts a higher ΓNa and a lower ΓMg than that predicted by the original MCTBI model. The predictions by the gMCTBI model are in better agreement with the experimental data.[78] Unlike the DNA helix, the tRNA has a more complicated compact structure, so we use the distance x to the RNA surface as the x-axis to draw the linear density. According to the linear density for Na+ shown in Figure A, at [Mg2+] = 0.6 mM, the gMCTBI-predicted Na+ linear density is higher than the original MCTBI-predicted result. The reason is that the gMCTBI model can treat the discrete Na+ ions, whereas the previous MCTBI model cannot. For the complicated 3D structure of the tRNA, Na+ ions (green balls in Figure B) can be trapped in the pockets that cannot be occupied by the bulkier Mg2+. The trapped Na+ ions would cause an effective (excluded volume and electrostatic) repulsive force on the TB Mg2+ ions near the pockets. Such scenarios, which can be well predicted by the current gMCTBI model, cannot be treated by the original MCTBI model. In a solution with higher [Mg2+], more TB Mg2+ ions would be found in the TB region. The stronger Mg2+–Na+ (repulsive) correlation effect tends to exclude/displace the trapped Na+ ions (see the case of [Mg2+] = 6 mM in Figure B). As a result, the difference between the gMCTBI and MCTBI predictions becomes smaller.

Figure 6

(A) Linear density of Na+ ions around the tRNA predicted by the gMCTBI (red lines) and the original (blue line) MCTBI models for solutions containing 10 mM background Na+ mixed with 0.6 mM competing Mg2+ (solid lines) and 6 mM competing Mg2+ (dash lines), respectively. (B) Most probable distributions of the TB ions in a solution containing 10 mM background Na+ mixed with Mg2+. The red and green balls denote the Mg2+ and Na+ ions, respectively

Conclusions

To summarize, we have developed a new model, the gMCTBI model, to predict the ion binding properties for nucleic acids. The model can treat discrete ion positions and correlations and ion fluctuations for both multivalent and monovalent ions. Theory–experiment tests for ion binding properties in monovalent and monovalent–divalent solutions support the validity of the gMCTBI predictions. In a solution with (single or mixed) monovalent ions, although the gMCTBI and the mean-field methods predict similar global binding properties (such as the excess number of monovalent ions), the gMCTBI model, which can treat discrete, correlated monovalent ion distributions, may provide more reliable predictions for the linear density for the bound ions and the most probable ion distributions. Furthermore, the gMCTBI model distinguishes from the original MCTBI model in the treatment of monovalent ions. Therefore, the differences between the predictions from the two models reveal significant effects of the monovalent ion-involved correlations.

For a mixed solution of different species of monovalent ions, the gMCTBI model, which accounts for discrete monovalent ions, can account for the ion size effect in ion binding accessibility. Smaller ions can enter narrow grooves that are not accessible to bulkier ions and hence are more competitive in ion binding. For a mixed monovalent–divalent ion solution, the gMCTBI model shows that divalent ions can be more than 9-fold more competitive than monovalent ions. For example, for a 24-bp DNA helix, the model predicts the ion concentrations for reaching equal number of the excess ions: 54 mM Na+ versus 6 mM Mg2+, 49 mM Na+ versus 5 mM Mg2+, and 20 mM Na+ versus 1.5 mM Mg2+.

The interplay between the different species of ions depends on the 3D structure of the nucleic acid. For example, a less bulky (hydrated) monovalent ion can more likely bind to RNA narrow pockets, which may be inaccessible to bulkier divalent ions. Such phenomena become more significant for more complicated nucleic acid structures such as tRNA, which can contain various narrow pocket sites for ion binding. Only at high divalent ion concentrations can the monovalent ions be displaced from the binding pockets (due to the repulsion from the divalent ions).

The framework of the gMCTBI model may be applicable to other systems such as highly charged surfaces and charged proteins. However, as a caveat, we note that the current form of the model ignores anion-involved correlations. Recent experiments[19,55] reported that a monovalent cation and a monovalent anion could form an ionic cluster near the DNA surface, indicating that the ion-involved correlation may be important in ion–nucleic acid interactions. Therefore, future development of the model would include anion correlation effects by explicitly sampling the distribution of discrete anions.

Materials and Methods

Preparations for the Theoretical Predictions

We use the modeling of ion effects for two nucleic acid structures to illustrate the gMCTBI: a 24-base pair (24-bp) B-form DNA and a 76-nt yeast tRNAPhe. The 3D structure of the 24-bp DNA is generated from the latest version (version 2.1) of the program X3DNA,[79] whereas the structure of tRNAPhe is downloaded from the protein data bank (PDB)[80] and its PDB ID is 1TRA.[81] In this study, we use the same sequences as the ones used in the experimental study for the B-DNA helix:[55] S1 = 5′-GGT GAC GAG TGA GCT ACT GGG CGG-3′ and S2 = 5′-CCG CCC AGT AGC TCA CTC GTC ACC-3′. A coarse-grained (CG) charge model is used[72,73] where each phosphorus atom carries an electronic charge (−e) and other atoms are treated as neutral. The CG charge model, however, cannot account for the sequence-dependence of the charge distribution on the nucleic acids. Therefore, the predictions from the gMCTBI model are dependent on the nucleic acid structure but not the sequence. Here, we should note that because the 24-bp DNA carries −46e net charges,[17,19,55] we manually deleted the 5′ terminal phosphates in both chains of the 24-bp DNA helix.

To reduce the boundary effect, the nucleic acid is placed in a large solution box, whose size is larger than 6 times Debye length. To test the model, various salt conditions are used, including single salt and mixed salt (with two cation species) conditions. All solutions contain Na+ either as the single species cation or as the background in a mixed salt solution. For the mixed salt, the other cation species may be monovalent ions, Li+, Cs+, and tetramethylammonium ion (TMA+), or divalent ions Mg2+. The bulk concentrations in the solution satisfy the charge neutrality condition: ∑(Zc)+ = c–, where Z and c denote the charge and the bulk concentration of the cation species i, respectively, and c– is the bulk concentration of an anion (Br– or Cl–). Ions are considered to be fully hydrated in our calculations. For Na+, Mg2+, and Cl–, their (hydrated) ion radii are 3.5, 4.5, and 4.0 Å, respectively. We note that these are the same parameters used in the previous studies.[65,72,73] The Mg2+ and the Na+ ion radii correspond to two and one hydrated shells, respectively.[82,83] Previous experiments[19,55] suggested that as the anion Cl– is replaced with Br–, the ion binding effects remain nearly the same. Thus, we set the effective Br– radius the same (4.0 Å) as that of Cl–. For other alkali cations Li+ (3.1 Å) and Cs+ (4.1 Å), their radii are estimated from the difference in van der Waals radius[84] between the cation (Li+ or Cs+) and Na+ and the hydrated Na+ radius.[85,86] Additionally, the radius of the (hydrated) organic cation TMA+ is assumed to be 4.65 Å.[87]

Modeling the Monovalent Ion Correlation Effect

The gMCTBI model is developed based on the original MCTBI model[72] (see SI for a brief summary of the MCTBI model). The major difference between the gMCTBI model developed here and the original MCTBI model is the different treatments of the monovalent cations.

TB Regions for the Different Types of Cations

The NLPB calculation is employed to roughly estimate the ion concentration around a given nucleic acid. Based on the ion concentration, the ion correlation strength at a given position x can be calculated ashere, superscripts (or subscripts) i and j denote the ion species. In the above equation, i and j can be the same or different ion species. z, kB, T (=25 °C in our calculations), and εw (≈78 at room temperature) are the charge (valency), Boltzmann constant, temperature, and solvent dielectric constant, respectively. aws(x) and c(x) denote the Wigner–Seitz radius[88] and the local ion concentration at position x, respectively. c0 is the bulk concentration in the solution. The region with ion correlation strength C(x) larger than 2.6 (the critical value for the gaseous to liquid phase transition for an ion system[42]) is defined as the TB region.

One of the important features that distinguish the gMCTBI model from the original MCTBI model is the different TB regions. Here, we use a 24-bp B-DNA in a mixed solution of [Mg2+] = 10 mM and [Na+] = 100 mM as an example. In the original MCTBI model, monovalent ions are completely excluded from the consideration of ion correlation, and as shown in Figure A, the TB region is defined only for Mg2+ ions. In contrast, in the gMCTBI model, as shown in Figure B, we distinguish three parts of the TB region according to the different ion species (i.e., different i and j’s in eq ). For a mixed divalent/monovalent (such as Mg2+/Na+) salt solution, the M22 region is the part for Mg2+ only (i = j = Mg2+ in eq ), the M12 region is the part for both Mg2+ and Na+ simultaneously (i = Mg2+ and j = Na+ in eq ), and the M11 region is the part for Na+ only (i = j = Na+ in eq ).

Figure 7

TB regions in the (A) original MCTBI model and (B) gMCTBI model. The 24-bp DNA structure is labeled by blue color. The black, cyan, and pink regions correspond to the TB regions for Na+–Na+, Na+–Mg2+, and Mg2+–Mg2+ correlations, respectively.

The three parts of the TB region are determined separately by eq . Basically, the M22 region is usually the largest because the Mg2+–Mg2+ correlation is the strongest. The M12 region is accessible to both TB Na+ and Mg2+ ions. Because of the weaker Na+–Mg2+ (than Mg2+–Mg2+) correlation, the M12 region is thinner than the M22 region. The M11 region, which may be accessible to Na+ but not to Mg2+ due to the bulky size of a fully hydrated Mg2+, is the closest to the surface of the nucleic acid and the thinnest (<0.5 Å in thickness[42]). In fact, in a dilute monovalent cation solution, the M11 and M12 regions could disappear. To make sure that there is enough space for explicit sampling of the discrete distributions of monovalent cations, we set layers with 0.5 Å minimum thickness for the M11 and M12 regions, respectively. In addition to the radii of atoms and the hydrated ions, the TB Na+ can be found at a distance ∼15 Å from the helix axis. In particular, for the M11 region, if only based on the ion correlation strength calculated from eq , the thickness of the region is much less than 0.5 Å even at high [Na+]. Therefore, in the solution with pure monovalent ions, almost all of the TB ions are located in the M11 region. However, for these TB ions, the correlation exists although it may be weak. The M11 region is set up such that, compared with the implicit method, the explicit treatment for ions can provide a higher-resolution description for the bound ions near the nucleic acid surface.

Partition Function

To take into account the monovalent ion effects, in the gMCTBI model, we compute the partition function (PF) of a nucleic acid-solution in a way such that the different species of ions (monovalent and multivalent ions) are treated equally. Specifically, for a given nucleic acid with Np nucleotides, the PF is calculated fromhere, superscripts (or subscripts) i and j represent the ion species. Nb and Nmax denote the number and the maximum saturated number of TB ions in the TB region, respectively. Similar to the original MCTBI model,[72] we assume that the maximum allowed TB ion numbers correspond to the reversal of the net nucleic acid chargeHere, we note that Nmax and Nmax are determined by the MCID process described in the following text.

Z(Nb, Nb) in eq is the partition function over all of the ion distributions with Nb TB ion i and Nb TB ion j in the TB region. Z(Nb, Nb) is computed ashere, Zid is the PF for the (ideal) solution without the nucleic acid and c0 denotes the bulk concentration of the ion. The product “∏” corresponds to the insertion process of the TB ions. w(n) is the statistical weight arising from the nth step of the insertion and is dependent on the pre-existing ion distribution in the TB region (see eq S1 in SI). ΔGd is the electrostatic free energy for the DB ions (see eq S3 in SI). One of the key ingredients in the model is the sampling of the TB ion distributions using the MCID algorithm, in particular, for the monovalent ions; see below.

Monte Carlo Insertion Deletion (MCID) Algorithm

We consider a divalent/monovalent ion mixed solution. To sample the TB ions, we first insert the ions one by one into the TB region. For each insertion step, we run two-step Monte Carlo sampling: we first randomly select the ion type and then randomly choose a position to insert the selected ion. According to eq , for the nth inserted cation, the free energy change of the whole nucleic acid-solution system can be calculated ashere, n could be n for the species i or n for the species j. The first and second terms describe the interaction energy changes for the TB ions (in the TB region) and for the DB ions (in the DB region), respectively. Here, ΔΔGd(n – 1 → n) = ΔGd(n) – ΔGd(n – 1), where ΔGd(n) and ΔGd(n – 1) are the electrostatic free energies of the DB ions with n and n – 1 TB ions in the TB region, respectively (see eq S3 in SI). To enhance the computational efficiency, for the calculation of ΔΔGd, we assume that a newly inserted TB ion is nonspecifically condensed to the nucleic acid; thus, ΔΔGd is a function of the number of TB ions. This approximation allows us to precalculate ΔΔGd to save the computational time.

We randomly select the species i (or j) for the inserted TB ions according to the following probabilityAfter determining the species of the inserted ion, we randomly place the ion in the TB region according to the following probability (for ion n placed at site k)here, U(k) represents the interaction energy change upon the insertion of an ion n at site k (see eq S2 in SI). m is the number of available (vacant) sites in the M11 and M12 TB regions for monovalent cations and in the M12 and M22 regions for divalent cations. Based on eqs –7, we insert the TB ions one by one until the TB region is saturated. In fact, the maximum number of the TB ions Nmax for species i and Nmax for species j is determined by the total charge limitation (eq ) and the insertion process (eqs –7). Here, eqs and 6 determine the type of the inserted ion (species i or j) and eq gives the site (location) of the inserted ion. Equation determines whether the insertion process should be terminated. If the total charge of the TB ions reaches the maximum saturated value (2Np with Np as the number of nucleotides), the insertion process is terminated and the subsequent deletion process is started.

In the above insertion process, as a new TB ion is inserted, the pre-existing TB ions cannot adjust their distribution. To allow the ions to redistribute, we re-sample the ion distribution by randomly removing ions one by one such that higher-energy ions have a higher probability to be removed. Specifically, the probability for removing ion n′ (at site k) from an n-TB ion system is determined by the energy change ΔU(k) for the processFrom the (n – 1)-TB ion distributions generated in the above deletion procedure, we compute the statistical weight w(n) using eq S1 in SI. We note that unlike the original MCTBI model, the current gMCTBI model considers the sampling of both monovalent and multivalent TB ions through an extra sampling process, namely, the random sampling of the ion species for the ions to be inserted (or deleted) based on eqs and 6.