Electrostatic free energy landscapes for nucleic acid helix assembly.

Metal ions are crucial for nucleic acid folding. From the free energy landscapes, we investigate the detailed mechanism for ion-induced collapse for a paradigm system: loop-tethered short DNA helices. We find that Na+ and Mg2+ play distinctive roles in helix-helix assembly. High [Na+] (>0.3 M) causes a reduced helix-helix electrostatic repulsion and a subsequent disordered packing of helices. In contrast, Mg2+ of concentration >1 mM is predicted to induce helix-helix attraction and results in a more compact and ordered helix-helix packing. Mg2+ is much more efficient in causing nucleic acid compaction. In addition, the free energy landscape shows that the tethering loops between the helices also play a significant role. A flexible loop, such as a neutral loop or a polynucleotide loop in high salt concentration, enhances the close approach of the helices in order to gain the loop entropy. On the other hand, a rigid loop, such as a polynucleotide loop in low salt concentration, tends to de-compact the helices. Therefore, a polynucleotide loop significantly enhances the sharpness of the ion-induced compaction transition. Moreover, we find that a larger number of helices in the system or a smaller radius of the divalent ions can cause a more abrupt compaction transition and a more compact state at high ion concentration, and the ion size effect becomes more pronounced as the number of helices is increased.

INTRODUCTION

Nucleic acid (DNA and RNA) structures and structural changes are critical for biological functions (1,2). Due to the polyanionic nature of the nucleic acid molecules, the thermal stability and folding kinetics of nucleic acids are strongly dependent on the ionic condition of the solution (1–14). A rudimentary question in nucleic acid folding is how the negatively charged helices pack against each other to form a compact 3D structure, and how different ions and ion concentrations affect the structure and stability of the helix assembly. Experiments have suggested that multivalent ion-mediated force plays an important role in DNA condensation (15–28) and compaction of other polyelectrolyte molecules (29–32). Moreover, extensive experiments have demonstrated the particularly efficient role of multivalent ions in RNA compaction (33–39). For example, it was found that [Mg2+] in millimolar (mM) or [Na+] in molar (M) can stabilize a compact intermediate state (33–44), and the compact state induced by Mg2+ is more compact and ordered (39,41).

Recent experimental studies have significantly advanced our understanding of RNA-ion interactions. For example, several experiments demonstrated that cations with higher valency and smaller size are more efficient in RNA-folding and such ions could also slow down the folding by stabilizing possible misfolded states (8,33–36). Moreover, a series of experiments indicated that charge neutralization of the nucleotide backbone may be responsible for the initial compaction of RNA (39–43). A recent small angle X-ray scattering (SAXS) measurement for a simple two-helix system suggested the existence of an attractive inter-helix force at high [Mg2+] and a repulsive force at very low [Mg2+] (45). Furthermore, a different SAXS experiment for loop-tethered helix system showed that the inter-helix force at high [Mg2+] is not sufficient to cause a collapsed state for the system (46). Based on the (mean-field) non-linear Poisson-Boltzmann (PB) theory, a thermodynamic framework for Mg2+–RNA interaction suggested that the predominant electrostatic contribution to RNA-folding stability comes from the diffusively bound ions (9,46–48).

Inspired by the experiments on the ion effect in nucleic acid folding (15–43,45,49), we investigate the ion-dependence of the folding free energy landscape for a paradigm system: two or three loop-tethered helices immersed in an ionic solution. The motivation of the present study is to quantitatively understand (i) the physical mechanism of the ion-induced compaction of the helix assembly and (ii) how the length and the electrostatic properties of loops, the multiplicity of the helices, and the ion size determine the stability and structural changes of the system. The free energy landscape, which gives the complete free energy distribution in the conformational space, enables detailed study of the physical mechanism of the ion-induced helix–helix assembly, such as the stability and the compactness of the collapsed state and the cooperativity of the compaction transition.

One of the key issues in the present study is the ion-mediated helix–helix interactions (19,20,28,33). Previous computer simulations (50–54) and analytical calculations (55–60) for similar systems reveal the existence of helix–helix attraction in the presence of multivalent ions (31,51). However, these previous studies are mostly based on simplified nucleic acid models and do not take into account the effects of the tether loops and the added salt in the solution. The mean-field PB theory (61–65) predicts no helix–helix attraction even for multivalent ion solutions. The counterion condensation theory (66) can predict the inter-helix attraction arising from the increased ion entropy (67–69) for both monovalent and multivalent ions.

Recently, we developed a statistical mechanical model, the tightly bound ion (TBI) model (70–72). The advantage of the model is that it can account for correlations and fluctuations, i.e. different ion-binding modes, for the bound ions. The model is based on realistic helical structure of nucleic acids and takes into account the effect of the added salt in the supporting solutions and the realistic ion sizes (70–72).

MATERIALS AND METHODS

Structural model

We focus on two model systems: two helices tethered by a semiflexible loop [Figure 1A and (49)] and three helices tethered by two loops (Figure 1B). To simplify the calculation, we treat only symmetric helix–helix configurations and assume that the axes of the helices are co-planar (2D). As a result, the system can be characterized by two structural parameters: the angle θ and the distance x between the helices. The helices are assumed to adopt the B-form structure and are described by the grooved primitive model (70–73), where the charged phosphate groups, the neutral groups and the ion-inaccessible helix core regions are represented by charged spheres, neutral spheres and ion-inaccessible cylinders, respectively. The coordinates and sizes of the spheres and the cylinders are determined from the experimental atomic coordinates of a B-form helix; see Supplementary Data. We choose the above described two- (or three-) helix system because (i) the system is simple thus ideal for detailed analysis and (ii) there are available experimental results for a similar system (49).

Figure 1

Symmetric configurations for two (A) and three (B) DNA helices tethered by one (or two) semiflexible loop(s) can be characterized by two structural parameters: the angle θ between the helical axises and the end-end distance x. The DNA helices are produced from the grooved primitive model (see Materials and Methods) (70,72,73).

Free energy landscape

For a configuration described by (x, θ) for the nucleic acid model, the free energy G (x, θ) is calculated as where GE (x, θ) is the electrostatic free energy for the helices immersed in ion solutions and Gloop (x, L) (L = loop length) is the entropic free energy of the loop(s).

We compute GE (x, θ) from TBI model (70–72), which explicitly accounts for the correlation and fluctuation of the bound ions, and we obtain Gloop (x, L) from analytical expressions from polymer theories (74): where Ploop (x, L), the probability of finding a loop with end-end distance x is given by an analytical expression derived from the polymer theory (74): where A is thenormalization factor. We note that calculations based on a different analytical formula for Ploop (x, L) (75) would give very close results for our calculations.

A loop of polyethylene or PEG chain used in the experiment (49) is neutral and its persistence length is ion-independent. In contrast, a realistic polynucleotide loop is polyanionic and its persistence length lP (and the loop entropy) is sensitive to the ionic condition (76–79). We account for the ion-dependence of the entropy of a polynucleotide loop through an empirical ion-dependent persistence length derived from the experimental measurements (76–79).

The distribution of free energy G (x, θ) (Equation 1) for different (x, θ) states gives the free energy landscape of the system. We explore the landscape for conformations described by x ∈ [22 Å, 60 Å] and θ ∈ [0°, 180°]. Such an ensemble of conformations would include a broad range of different types of conformations, including the compact, the random and the extend conformations. We choose 22 Å as the lower limit for x because this is the closest inter-axis distance of approach between two helices. In the actual depiction of the landscape, we scan the (x, θ) space with increments of 2 Å and 20° for x and θ, respectively. For each given (x, θ) state, we enumerate all the accessible 3D conformations through the rotation of the planar rigid body with the fixed x and θ. We use such an ensemble of 3D conformations to compute the free energy G (x, θ).

TBI model

For a given helix–helix configuration (x, θ), we calculate the electrostatic free energy using the TBI model. The development of the TBI model is motivated by the need to accurately account for Coulombic and excluded volume correlations between bound ions and between bound ions and helices as well as the fluctuations of the ion-binding modes (70–72); see Supplementary Data and Refs (70–72) for details.

In the TBI model, according to the ion–ion correlation strength, we classify the tightly bound (= strongly correlated) and the diffusive (= weakly correlated) ions/regions (70–72). We treat diffusive ions using non-linear PB and treat the tightly bound ions separately by exhaustively enumerating all the possible distributions of the tightly bound ions (= ion-binding modes). Each mode M is described by a set of numbers with each number m equal to the number of the tightly bound ions surrounding the ith phosphate. For each mode, we compute the partition function Z by evaluating the energy and entropy of the diffusive ions, the tightly bound ions, and the interactions between the diffusive and the tightly bound ions (70–72). In the calculation, the many-particle Coulombic and excluded volume correlations between the tightly bound ions and between the tightly bound ions and the helices are accounted for. Sum over all the possible ion-binding modes gives the total partition function for a given nucleic acid structure, from which the electrostatic free energy is obtained. We assume Na+ and Mg2+ ions are hydrated and have radii 3.5 Å and 4.5 Å (70–72), respectively. When solving the non-linear PB, we add a thin charge-free layer of one cation radius to the molecular surface to account for the volume exclusion of the cations (61,70–72). The dielectric constant ɛ is assumed to be 12 [experimental value (80)] for the helix interior and 78 (of bulk water at 25°C) elsewhere in the solution.

SAXS profile

SAXS can effectively reveal the shape and size of macromolecules (42–49,81,82). To compute the SAXS profile I (Q) with the grooved primitive helix model, we use the following equation for I (Q) (83): where M0 is the total number of the charged and neutral spheres in the model, and r is the distance between two spheres. The scattering vector Q is equal to 4πsin(θ)/λ, where θ is the Bragg angle and λ is the wavelength of the radiation. f is the scattering factor (84) , where Q is in the range between 0 and Qmax, the largest scattering vector in reciprocal space (84). For simplicity, we use the same scattering factor f (Q) for all spheres in the model (84). Qmax is taken as 3 Å−1 in the calculations. Our tests with different Qmax (≳2 Å−1) give quite robust results.

RESULTS AND DISCUSSIONS

We consider a system of loop-tethered DNA helices immersed in an ionic solution. Each helix is assumed to have a finite length of 10 bp. To simplify the analysis, we assume that the helices are symmetrically oriented, so the system can be characterized by two structural parameters: angle θ and distance x between the helices (see Figure 1). The stability of a state described by (x, θ) is determined by its free energy G (x, θ) (see Equation 1). G (x, θ) for all the possible (x, θ) states with x ∈ [20 Å, 60 Å] and θ ∈ [0°, 180°] gives the free energy landscape.

From the free energy landscape G (x, θ), we can identify all the possible stable states, which are represented by the minima on the landscape, and the transitions between the different stable states. We characterize the compactness of the helix assembly by a Boltzmann-weighted mean value for the distance r between the centers of the helix pair: As described in Methods, in our calculation, we choose x ∈ [22 Å, 60 Å] and θ ∈ [0°, 180°] with increments of 2 Å and 20° for x and θ, respectively.

According to the compactness of the structure, we distinguish three types of states: (i) the extended state E for configurations with large x and large θ, (ii) the random fluctuation state R for configurations with random x and θ and (iii) the collapsed state C for configurations with x ∈ [22, 26]Å and θ ∈ [0°, 40°] (69). Our calculation shows that the collapsed state C has a mean compactness of , which agrees with the upper limit for the helix–helix separation for a collapsed (aggregated) DNA array (18–28). It is important to note that such an electrostatically collapsed state is less compact than the final RNA native structure (44), where further compaction can be induced by specific tertiary interactions that are ignored in the present study.

The key to understanding ion-induced structural transition is the free energy landscape G (x, θ). As shown in Equation 1 in ‘Methods’, G (x, θ) of the system consists of two parts: (i) the electrostatic free energy GE (x, θ) defined as the loop-free free energy for the helices and the ionic solution and (ii) the loop-free energy Gloop (x) for the loop conformational entropy. Before including the loop effect, we first focus on the effect of the electrostatic free energy landscape GE (x, θ).

Without the tether loop

Na+ ions

In Figure 2A–D we show density plots for the electrostatic free energy landscapes GE (x, θ) for two helices immersed in a Na+ solution. At very low [Na+] (∼1 mM), the favorable low-free energy state is the extended state E with the largest x and θ so that the repulsive helices can maximally avoid each other. As [Na+] is increased, the favorable state would ‘smear out’ to reach the region of smaller x and θ as shown in Figure 2B and C for 0.01 M and 0.1 M NaCl, respectively. When [Na+] becomes very high, such as 1 M, the favorable conformations can spread over nearly the whole region of x and θ. In Figure 3A, the solid lines show the compactness as a function of [Na+]. We find that the increase of [Na+] can induce the compaction process. The rapid decrease in occurs for [Na+] 0.3 M. For [Na+] >0.3 M, the helix–helix system approaches the random fluctuation state R. Physically, as [Na+] increases, with the decrease in the entropic penalty for ion-binding, more Na+ ions can bind to the helices, resulting in a weakened repulsion between the helices. The predictions show that Na+ can cause a random fluctuation state at high [Na+] (42–49), and the random fluctuation state is more compact than the electrostatically repulsive unfolded state at low [Na+]. Such high concentration Na+-induced compaction has been observed in experiments for RNAs (33–43). For larger RNA molecules, experiments also show that Na+ can induce the partial formation of some native tertiary contacts (33–43).

Figure 2

The normalized probability exp[−G (x, θ)/k] for loop-free helix–helix system. Here G (x, θ) is the electrostatic free energy landscape for different (x, θ) configurations in a solution of NaCl (A–D), MgCl2 (E–H), or divalent ion (M2+) with (small) radius = 3.5 Å (I–L). The red and blue colors represent the low and high free energies, respectively.

Figure 3

(A and B) Ion concentration-dependence of the compactness (= the mean distance between helix centers) for two helices without loop (solid lines), with neutral polyethylene loop (dotted lines) and polynucleotide loop (dashed lines) in NaCl (A) and MgCl2 (B) solutions. (C and D) The compactness for three helices. Results for different end-to-end distance maximum (= loop length L): xmax = 60 Å, 44 Å and 28 Å (from the upper to lower) are plotted, respectively. In A–D, the gray solid lines denote the compactness of the random fluctuation states R for different xmax's and gray dashed lines denote the compactness of the collapsed state C. The compactness is computed from the random averaging over x ∈ [22,26]Å and θ ∈ [0°, 40°] (69) for state C and x ∈ [22 Å, xmax] and θ ∈ [0, θmax] (θmax = 180° for a two-helix system and 120° for a three-helix system) for state R. We note that for state C, and the result is in agreement with the helix–helix separation for an (least compact) aggregated DNA array (18–28).

In each panel of Figure 3, we show three curves, each corresponding to a different constraint for the maximum helix–helix distance xmax. Physically, xmax is determined by the length of the linker (loop) that tethers the helices. The results for different xmax's indicate how the spatial constraint from the finite loop length affects the two-helix electrostatic free energies and the resultant compactness of the system. In general, for low ion concentration, the two-helix configuration tends to become extended. So the compactness is more sensitive to xmax (= the upper limit of the extension). Larger xmax causes a more extended state. For high ion concentration, the system becomes more compact and thus the compaction is less sensitive to xmax. As a result, larger xmax causes a more abrupt (ion-induced) change in the compactness, i.e. a larger xmax causes a sharper compaction transition.

Mg2+ ions

For [Mg2+] below 1 mM, as [Mg2+] is increased, due to the stronger charge screening and consequently the weakened repulsion, helices become more compact as shown by the shift of the free energy minima from Figure 2E to F.

For [Mg2+] above 1 mM, more ions are bound and the correlation effect becomes strong. The bound ions can self-organize to reach correlated low-free energy configurations. Such correlated energies decrease rapidly as helices approach each other, leading to an inter-helix attractive force (45,72). Higher [Mg2+] causes more bound ions and stronger inter-helix attraction (72), and hence a more compact state (as shown by the shift of the free energy minima from Figure 2F to H and the change of the compactness in Figure 3B).

Because the inter-helix attraction is strongest for two parallel helices (53), the free energy minima occur at θ = 0 in Figure 2F–H. As shown in Figure 3B, the Mg2+-induced compact state is more compact than the random fluctuation state R but less compact than the collapsed state C. Such effect is more notable for larger xmax; see Figure 3 for the curves for three different xmax. Moreover, for very close approach between the helices, the parallel helices would repel each other due to the Coulombic repulsions between helices and crowded bound ions (72). In such close-packing case, the helices could be stabilized through rotation in order to minimize the repulsion; see the (rotated) local minimum at (θ, x) ≃ (40°, 22 Å) in Figure 2H. As [Mg2+] exceeds 0.1 M, the helices may be over-neutralized by the Mg2+ ions and the compaction becomes saturated and close to the random fluctuation state (R) (especially for smaller xmax), due to the Coulombic and excluded volume repulsions between the bound Mg2+ ions (27,28,72).

Both Na+ and Mg2+ can facilitate helix–helix compaction, but through different mechanisms. High concentration of Na+ causes compaction through strong ionic screening for helix charges, while Mg2+ ions can cause the compaction through helix–helix attraction (72). Moreover, Na+-induced compact state is heterogeneous as shown in Figure 2D while Mg2+-induced compact states are ordered as shown in Figure 2G and H.

The overwhelming efficiency of Mg2+ over Na+ is much more pronounced for tertiary than secondary structure folding. For example, we find that for a 9 bp single helix, 10 mM Mg2+ and 1 M Na+ can provide the same stability (71). For helix–helix assembly as a tertiary fold, however, our results show that a 10-fold dilute [Mg2+], i.e. 1 mM [Mg2+], would be sufficient to provide the same stability as 1 M Na+. This suggests that Mg2+ tends to promote the formation of the tertiary structure so the free energy can be more efficiently lowered by the Mg2+ ions. The result is in accordance with the experiments on Na+ and Mg2+-mediated RNA compaction (33–43). It is important to note that the ‘anomalous’ efficiency of the Mg2+ cannot be merely attributed to the more efficient charge screening, as shown by the PB-predicted electrostatic free energy landscape for Mg2+ in Supplementary Figure 9, which gives a disordered and much less compact state even at high [Mg2+].

It is important to note that, however, the Mg2+-mediated helix–helix interaction is not strong enough and x0, the inter-helix distance at the free energy minimum, is not small enough to cause a collapsed state; see Supplementary Figure 10. In Supplementary Figure 10, we show the magnitude of ion-mediated helix–helix interaction in NaCl and MgCl2 solutions. For example, ΔGE = GE − GE (extended state) for 0.1 M Mg2+ is the lowest at (x0, θ) ≅ (32 Å, 0°) with ΔGE ∼ −3.3k. As the helices further approach each other (with x < x0), ΔGE for the parallel helices increases, and the minimum ΔGE is shifted from parallel helices (θ = 0) to a rotated helix–helix configuration with θ ≠ 0 in order to avoid electrostatic repulsion between the helices. For example, for x = 24 Å, the minimum ΔGE (∼−2.4k) occurs for a V-shaped conformation with θ = 40°. The rotation angle is larger for smaller x. For example, the minimum ΔGE (most stable configuration) occurs at θ ≃ 60° for x = 20 Å. When Mg2+ exceeds 0.1 M, the Mg2+-mediated attraction becomes weaker due to the increased likelihood for over-neutralization (70,72). For example, for [Mg2+] = 0.6 M, the free energy minimum occurs at (x0, θ0) ≃ (32 Å, 0°) with ΔGE ∼ −2.7k, which is weaker than −3.3k, the minimum free energy for 0.1 M Mg2+. Similar multivalent ion-induced destabilization has been observed experimentally (22,27,28,30,53). For example, packed DNA assembly has been found to be dissolved by high multivalent ion concentration (22,27,28). The predicted ΔGE is slightly stronger than the experimentally derived value (49). The slightly overestimated stability in the model may be caused by the restriction of (2D) co-planar and symmetric conformations for the two-helix assembly and the possible change in the hydration state of the helices and ions.

With a neutral tether loop

To investigate the loop effect, we consider a system of two helices tethered by a polyethylene glycol (PEG) loop with persistence length lP ∼ 3.8 Å (49,85). We obtain the total free energy G(x, θ) from Equation 1 with the loop-free energy Gloop(x, L) determined from Equation 2 for each given configuration (x, θ) for two helices tethered by the loop. We first compare our prediction with the experimental measurements for a similar nucleic acid model (49). To model the short loop in the experiment (49), we assume loop length L = 26 Å. The loop size defines the maximum helix–helix end-end distance: x ≤ xmax = 26 Å. We calculate the SAXS profiles from Equation 4 for each given (x, θ). The Boltzmann-weighted average by factor e− for all the possible states (x, θ) gives the mean SAXS profile: where the conformational sampling covers the full conformational space described by x ∈ [22 Å, 60 Å] and θ ∈ [0°, 180°] with increments of 2 Å and 20°, respectively.

Short neutral loop and comparisons with experiment

As shown in Figure 4, the predicted SAXS profiles for low (20 mM) and high (1.2 M) [Na+] and for high (0.6 M) [Mg2+] are close to the experimental profiles for a similar nucleic acid model (49). For comparison, in Figure 4A we also show the predicted profile for 0.1 M Mg2+. The state in 0.1 M Mg2+ is (slightly) more compact than that the random state in 0.6 M Mg2+ but much less compact than the collapse state.

Figure 4

The SAXS profiles for the two-helix model with neutral PEG loop for different ionic conditions: diamond 20 mM NaCl; plus 0.6 M MgCl2 and square 1.2 M NaCl. (A) The predicted SAXS profiles weighted by the free energy G(x, θ) (see Materials and Methods); (B) The experimental SAXS profiles (49). The dashed lines, solid lines and dotted lines represent the standard SAXS profiles for the collapsed C, random R, and extended E states of the tethered helices, respectively.

Figure 4 shows close SAXS profiles for 1.2 M [Na+] and 0.6 M [Mg2+], while Figure 2 shows very different free energy landscape for Na+ and Mg2+ ionic conditions. How can such seemingly contrasting results be consistent with each other? The SAXS profile I*(Q) is a weighted average over all the available conformations (see Equation 5). As shown in Figure 2D and H as examples, for a given loop length of 26 Å, the accessible conformational space on the free energy landscape is limited to the region with x ≤ 26 Å. Though the full (unrestricted) landscape would certainly yield drastically different I*(Q), the restricted accessible region of the landscape may give similar I*(Q). Physically, within the restricted conformational space, the conformational distribution for high [Na+] and high [Mg2+] are both quite inhomogeneous over a wide-range of the angle θ between the helices.

It is important to note that the above prediction-experiment comparisons involve several limitations of the present theory: (i) The two-helix configuration is assumed to be symmetric and planar. Including the non-planar/non-symmetric configurations would likely weaken the predicted helix packing stability and compactness. (ii) The loop (linker) in the present model connects two helices at the helical axes instead of the helix edges in the experimental model (49). Such difference causes no significant effects for long loops (linkers) but can play a role for shorter loops.

Furthermore, the length of L = 26 Å used in the theoretical model and the experimental system may cause additional effects. As shown in Figure 2A–H, the most stable (lowest free energy) configuration occurs at inter-helix distance >26 Å. As a result, all the configurations considered in the model and in experiment (with x ≳ xmax = L = 26 Å) would result in ΔG higher than the lowest free energy. Therefore, the short loop would effectively reduce the electrostatic force in the compaction process. Moreover, short loops would significantly restrict the accessible conformations of the extended state, resulting in a reduced conformational difference between the compact and the extended states, and hence a reduced sharpness of the compaction transition. For the same reason, longer loops can lead to sharper compaction transition; see Figure 3.

In summary, our results agree with the experimental finding that the Mg2+-mediated electrostatic interaction cannot induce the electrostatically collapsed state C for the two-helix system. Although our predicted helix–helix attraction is slightly stronger than the value derived from experiments, the predicted compact state is only slightly more compact than the random state but much less compact than the collapsed state C. A longer loop causes a sharper compaction transition. Furthermore, as we discuss below, a polynucleotide loop (instead of an electrically neutral loop), which has an ion-dependent flexibility and entropy, would also lead to a sharper transition. Nevertheless, despite the different types of loops considered, we find the compact state is much less compact than the collapsed state.

Long neutral loop

To further investigate the loop effect, we consider longer polyethylene (-CH2-) tether loops (49) with lP ∼ 6 Å (86,87). As shown in Figure 3 for different loop length (maximum helix–helix end-to-end distance), the loop promotes compaction in both Na+ and Mg2+ solutions in order to gain larger loop entropy, which occurs for smaller x. For example, at 10 mM [Mg2+] in Figure 2F, our free energy landscape shows that the inter-helix distance x for the most stable state decreases from 34 Å (without loop) to 32 Å with a loop of length L = 60 Å. As the loop length is decreased, the loop entropy effect is weakened because of the reduced conformational entropy of the loop (74,75).

As a caveat, we note that if the loop is too short, the helices would be restricted to the small x region to form a rotated state, which is less compact (larger ) than the parallel state, as discussed above.

With a polynucleotide loop

As described in Materials and Methods, a realistic polynucleotide loop is polyanionic. Unlike neutral loops, such as polyethylene or PEG chains, polynucleotide persistence length lP and the loop entropy are sensitive to the ionic condition (76–79). Based on the experimental data for the ion-dependence of lP (76–79), we obtain an empirical expression for lP as a function of the ionic strength I (in M) in the solutions (see Figure 5). The above expression is mainly derived from the experimental data for Na+ due to the limited experimental data for lP in Mg2+ solution. As an approximation, we use Equation 7 for both Na+ and Mg2+ solutions.

Figure 5

The persistence length lP of a single-stranded polynucleotide chain as a function of ionic strength I. The solid line is the fitted function Equation 7. The symbols are the experimental data: filled square (76), filled diamond (77) and filled triangle (78) are for Na+ solutions and open square (79) is for mixed Na+/Mg2+ solution.

As shown in Figure 3A and B for the ion-dependence of the compactness, for low ion concentration the loop is stiff (with large lP, so the system is much more extended than the neutral loop case, especially for large loops. For high ion concentration, the loop is flexible (with small lP) (76), so the helices collapse to a state as compact as the neutral polyethylene loop case. Therefore, the polynucleotide loop greatly enhances the sharpness of the compaction transition. Such effect is much more pronounced for Mg2+ due to the stronger ion concentration-dependence of lP. The Mg2+-induced compaction transition occurs in the [Mg2+] range of 1–10 mM. In addition, since Equation 7, which is fitted from the experimental data mainly for Na+, may underestimate the [Mg2+]-dependence of lP (79), the realistic [Mg2+]-induced transition could be even sharper than that predicted here.

With more helices and loops

Since realistic RNA-folding often involves multiple (>2) helices, it would be useful to investigate how the multiplicity of helices affects the compaction process. We consider a simplified system consisting of three helices tethered by two loops as shown in Figure 1B. A rigorous treatment for the system requires a full account of all the possible bound ion configurations for the three-helix system for all the possible (x, θ) configurations, which would be extremely demanding in computational time. So we parse the problem into two simplified sub-problems. We first treat the full (x, θ) configurational ensemble by assuming pairwise additive electrostatic free energy for the three-helix system: , where is the electrostatic free energy for a system of two helices i and j. We then examine the effect of the three-helix correlation by considering the system of three parallel helices (θ = 0), where the enumeration for all the possible triple-helix (rather than pairwise) correlated bound ion modes can be performed.

Pairwise approximation

With the pairwise (additivity) approximation, for each given symmetric three-helix configuration (x, θ), the ion-dependence of the electrostatic free energy as well the loop-free energy Gloop(x, θ) ≃ 2Gsingle loop(x) are ‘amplified’. Therefore, the three-helix system shows a qualitatively similar but quantitatively much stronger ion-dependence than the two-helix system; see Figure 3A–D. For example, at low [Mg2+], the extended state becomes more extended due to the enhanced helix–helix repulsion and the effect from the extended rigid loops, and at high [Mg2+], the compact state becomes more compact due to the stronger electrostatic attraction and the enhanced loop entropy effect.

Triple-helix correlation

By considering the interactions between the three helices (and their associated bound ions), we compute the triple-helix correlated electrostatic free energy landscape ΔGE(x, 0°) = GE(x, 0°) − GE(60 Å, 0°). As shown in Figure 6, the pairwise approximation slightly underestimates and overestimates the stability of the free energy minima for [Mg2+] below 1 mM and above 10 mM, respectively.

Figure 6

The electrostatic free energy ΔGE(x, 0°) = GE(x, 0°)−GE(60 Å, 0°) (in k) as functions of the helix–helix separation x for three parallel DNA helices (θ = 0) at different Mg2+ concentrations: 0.1, 1, 10 and 100 mM (from the upper to lower). Solid lines: calculated from the pairwise additivity approximation. Dotted lines: calculated with three-helix correlation in the TBI model.

For high [Mg2+], the helices are highly neutralized. For example, at 0.1 M Mg2+, tightly bound Mg2+ can neutralize about 92% of the phosphate charges for separated DNA helices. Compared to the two-helix case, the freedom of the bound ions for the helix triad is much more restricted due to the ion–helix and ion–ion Coulombic and excluded volume interactions. Such triple-helix restriction tends to prevent the (frustrated) bound ions from self-organizing to reach the lowest energy states that would be realized for a two-helix system. As a result, the resultant most compact state is less stable and slightly less compact than that predicted from the pairwise additive model. On the other hand, for lower [Mg2+], the number of bound ions is smaller and the charge neutralization is weaker, so the bound Mg2+ ions can self-organize relatively freely to reach the correlated low-energy states (88). In such case, due to the higher charge density for the three-helix case, the system can reach an energy lower than that predicted from the additive two-helix model (88).

Our above predictions agree with the previous analysis based on simplified models (68,88), which showed that the non-additivity effect is not strong for three-body system, and the non-additivity can result in slightly stronger attraction than the pairwise approximation for low ion concentration (68,88).

Ion size effect

Ion size has been shown experimentally to play important roles in RNA compaction (35). To investigate the size effect, we compare the free energy landscapes and the compactness for two divalent ions of different sizes: Mg2+ with radius ∼4.5 Å (70,72) and divalent ion (denoted by M2+) with a smaller radius 3.5 Å. Due to the lack of the information about the ion size-dependence of the loop persistence length, we would focus on the loop-free system.

Two-helix system

By comparing the free energy landscapes in Figure 2E–L and the compactness in Figure 7A for Mg2+ and the small divalent ions, we find that the smaller ions can induce stronger helix–helix attraction and the resultant stable states are more compact (72). Therefore, small ions are more efficient to induce the compaction. Physically, smaller ions can bind to the helices at a smaller ion-helix distance. Stronger ion–helix interaction leads to a higher ion-binding affinity and hence a larger number of bound ions, which, as a result, causes a stronger helix–helix attraction (72) and thus a more compact state.

Figure 7

Compactness of two loop-free helices in Mg2+ (radius ≃ 4.5 Å) and smaller divalent ion (M2+, radius = 3.5 Å) solutions, respectively. Results for different end-to-end distance maximum: xmax = 60 Å, 44 Å and 28 Å (from the upper to lower) are plotted separately. The gray lines denote the compactness of the random fluctuation states R for different xmax's.

Three-helix system

As shown in Figure 7B for the ion-dependence of the compactness, the higher efficiency of the smaller ions is much more pronounced for the three-helix system. Notably, the three-helix system with high concentration of smaller ions can reach a very compact state—the collapsed state C. Such a collapsed state can never be reached by a two-helix system or a three-helix system in the solution of the (more bulky) Mg2+ ions. Physically, this is caused by the stronger helix–helix attraction for smaller ions and the effect is further ‘amplified’ by the larger number of helices.

CONCLUSIONS AND DISCUSSIONS

Using the TBI model (70–72) to account for the correlation and fluctuation effects of bound ions, we explore the electrostatic free energy landscape for nucleic acid helix–helix assembly for several paradigm systems: two and three helices with and without tether loop(s). The computed free energy landscapes lead to several predictions for the ion-dependent compaction of the helices.

Adding Na+ to the solution would induce the nucleic acid helices to undergo a transition from an extended state to a more compact random fluctuation state. However, Na+-mediated electrostatic interactions cannot cause the helices to form the collapsed state even at very high concentrations. This prediction is in accordance with the experiments on Na+-mediated RNA compaction (42–49).

Mg2+ can mediate helix–helix attraction, which causes the helix–helix assembly to form a state more compact than the random fluctuation state but less compact than the collapsed state. Such effect would be reduced for shorter helix–helix linker. Mg2+ ions are much more efficient than Na+ to stabilize nucleic acid helix–helix assembly. For example, 1 mM Mg2+ can induce a state equally compact as that induced by 1 M Na+. The maximum compact state occurs at [Mg2+ around 0.1 M.

The tether loop(s) can facilitate the compaction transition. While an electrically neutral loop enhances compaction for both low and high ion concentrations, a polynucleotide loop can significantly enhance the sharpness of the ion-induced compaction transition due to the ion-dependence of loop flexibility. Such effect is more pronounced for Mg2+ (than Na+) and for larger loops.

Multiplicity of the helices and the loops can significantly enhance the sharpness of the compaction transition and the compactness of the compact state.

Smaller ions are more efficient to induce the electrostatic compaction and such effect is much more pronounced for a multiple helix assembly.

Our treatment here is limited by several simplifications and approximations. First, we have neglected the possible partially or completely (dehydrated) ion site-binding (9,89) and possible preferential interactions between ions and specific groups of the nucleic acids (1). Ion site-binding to specific pockets with very strong negative electric potential can help stabilize specific tertiary structures (46,47) and can be essential for RNA catalytic function (90). The preferential ion-binding to specific groups, e.g. Mg2+ binds to phosphate groups or Mn2+ binds to both phosphates and bases (1), would affect the spatial distribution of ions on nucleic acid surface and contribute to the ion-mediated helix–helix interactions. However, the overall contribution of ion site-binding to the global stability may be small, because the large dehydration penalty may reduce the likelihood of specific binding of the dehydrated ions (46,47).

Second, in order to simplify the computation, we assume each helix pair to be symmetric and planar. Such a system covers a full spectrum of typical configurations, including compact, random and extended states. Thus, our predicted free energy landscape reveals the main feature of the ion-dependent helix assembly. We have also ignored the effect of different rotational phase of the helices about their axes. Previous analysis showed that the free energy landscape is not sensitive to such rotations if the helix–helix separation is large (e.g. x ≥ 28 Å) (91). However, for two closely approached helices (x ≤ 24 Å), such rotation would play an important role because of the detailed geometrical compatibility between helices and the related charge distributions (69,92). For very closely approached helices, parallel configuration is unfavorable because of the strong electrostatic repulsion. In the present model, the helices would be stabilized through a rotation to form a V-shaped structure. In realistic full 3D conformational space, the two helices can be stabilized through off-plane rotations to form a typical X-shaped structure (69,92,93) with two possible packing modes: groove-into-groove mode where the major groove of one helix fits in the major groove of the other helix (69), and backbone-into-groove mode where the nucleotide strand of a helix fits in the major groove of the other helix (92). Helix rotation about its axis is important for the formation of such X-shaped structures (69,92,93).

Third, we have ignored the sequence effect and specific tertiary interactions, which are essential for RNAs to form specific functional structures. In this study, we only investigated the ion electrostatic effects on the electrostatic collapse process. The quantitative modeling of further folding steps requires the inclusion of the sequence-dependent tertiary interactions.

Finally, we have neglected the detailed loop structures and possible deformation and unfolding of helices at different ionic conditions. Nevertheless, the predicted free energy landscapes, the detailed compaction mechanisms for Na+ and Mg2+ ions, and the effects of different loops, multiplicity of helices and ion size may be useful for understanding the ion-dependent nucleic acid helix–helix assembly.

SUPPLEMENTARY DATA

Supplementary Data are available at NAR online.