---
abstract: 'Secondary structure formation of nucleic acids strongly depends on salt concentration and temperature. We develop a theory for RNA folding that correctly accounts for sequence effects, the entropic contributions associated with loop formation, and salt effects. Using an iterative expression for the partition function that neglects pseudoknots, we calculate folding free energies and minimum free energy configurations based on the experimentally derived base pairing free energies. The configurational entropy of loop formation is modeled by the asymptotic expression $ -{c}\ln {m}$, where ${m}$ is the length of the loop and ${c}$ the loop exponent, which is an adjustable constant. Salt effects enter in two ways: first, we derive salt induced modifications of the free energy parameters for describing base pairing and, second, we include the electrostatic free energy for loop formation. Both effects are modeled on the Debye-Hückel level including counterion condensation. We validate our theory for two different RNA sequences: For tRNA-phe, the resultant heat capacity curves for thermal denaturation at various salt concentrations accurately reproduce experimental results. For the P5ab RNA hairpin, we derive the global phase diagram in the three-dimensional space spanned by temperature, stretching force, and salt concentration and obtain good agreement with the experimentally determined critical unfolding force. We show that for a proper description of RNA melting and stretching, both salt and loop entropy effects are needed.'
author:
- 'Thomas R. Einert'
- 'Roland R. Netz'
bibliography:
- 'RNA\_salt.bib'
title: ' Theory for RNA folding, stretching, and melting including loops and salt'
---

Introduction {#sec:introduction}
============

Ribonucleic acid (RNA) is one of the key players in molecular biology and has in the past attracted theoretical and experimental physicists because of its intriguing structural and functional properties. RNA has multiple functions: beyond being an information carrier it has regulatory and catalytic abilities [@Gesteland2005]. Comprehending how RNA folds and what influences the folding process are key questions [@Tinoco1971]. Thus, the reliable prediction of RNA structure and stability under various conditions is crucial for our understanding of the functioning of RNA and nucleic acid constructs in general [@Liedl2007; @Dietz2009].

The influence of temperature and solution conditions on RNA folding stays in the interest of experimental groups. Traditionally the thermal melting of RNA was monitored [[via]{} ]{}differential scanning calometry or UV spectroscopy for the bulk ensemble [@Xia1998; @Mathews1999; @Privalov1978; @Vives2002]. More recently, single molecule pulling and unzipping experiments have been used to unveil the influence of different solution conditions and even determine energy parameters [@Liphardt2001; @Vieregg2007; @Huguet2010].

On the theoretical side, RNA denaturation has been modeled on various levels of coarse graining. Focusing on the secondary structure, namely the base pairs (bp), and omitting tertiary interactions, equilibrium folding and unfolding has been modeled very successfully [@McCaskill1990; @Zuker1981; @Markham2005; @Gerland2001; @Bundschuh2005; @Montanari2001; @Mueller2003; @Dimitrov2004; @Imparato2009a; @Einert2008; @Einert2011a]. In the presence of a logarithmic contribution to the loop entropy, it has been shown that homopolymeric RNA, where sequence effects are neglected, features a genuine phase transition, which can be induced by force or temperature [@Mueller2002; @Mueller2003; @Einert2008; @Einert2011a]. However, the specific sequence influences the stretching response of a molecule, which has been shown by @Gerland2001 [@Gerland2003], yet without considering the logarithmic loop entropy. More detailed insights can be obtained by simulations, which are numerically quite costly, though, when compared to models focusing only on secondary structure. Coarse grained, Go-like simulations of short RNA hairpins allowed to analyze the dynamics of the folding and unfolding process [@Hyeon2005; @Hyeon2006]. Ion specific effects have been studied by performing molecular dynamics [@Auffinger2000] or coarse grained simulations [@Tan2006; @Tan2007; @Tan2008; @Jost2009]. Much less is known about the salt dependence of denaturation transitions of RNA.

While for DNA numerous corrections of the base pairing free energies due to varying salt concentration exist, see [@Owczarzy2004] and references therein, analogous results for the salt dependence of RNA energy parameters are sparse [@Tan2006]. However, molecular biology and biotechnological applications depend on the reliable prediction of RNA stability for different solution conditions.

In this paper we extend these previous works and develop a theory that allows to include all these effects – sequence, salt dependence, logarithmic loop entropy, stretching force – and demonstrate that all are necessary to obtain a complete picture of the thermodynamics of the secondary structure of RNA. Neglecting tertiary interactions, we use a recursion relation, which allows to correctly account for logarithmic and thus non-linear free energy contributions due to the configurational entropy of loops [@Einert2008]. To include the influence of monovalent salt on RNA stability, we model the RNA backbone as a charged polymer interacting [[via]{} ]{}a Debye-Hückel potential and give heuristic formulas for the modification of the loop free energy and the base pairing and stacking free energy parameters. Debye-Hückel is a linear theory, yet we include non-linear effects caused by counterion condensation using Manning’s concept [@Manning1969]. The backbone elasticity of single stranded RNA (ssRNA) is described by the freely jointed chain (FJC) model. Our description allows for a complete description of the behavior of RNA in the three-dimensional phase space spanned by temperature, salt concentration, and external stretching force. We find that for an improved description of RNA melting curves one needs to include both salt effects and loop entropy. Only the combined usage of these two contributions enables to predict the shift of the melting temperature (due to salt) and the cooperativity (due to logarithmic loop entropy), which is illustrated in the case of tRNA-phe. As an independent check we consider the force induced unfolding of the P5ab RNA hairpin and observe good agreement with experimental values with no fitting parameters. The influence of salt is illustrated by melting curves and force extension curves for various salt concentrations. For the P5ab hairpin the phase diagram is determined and slices through the three-dimensional parameter space are shown.

Free energy parameterization
============================

RNA folding can be separated into three steps, which occur subsequently and do not influence each other to a fairly good approximation [@Tinoco1999]. The primary structure of RNA is the mere sequence of its four bases cytosine (C), guanine (G), adenine (A), and uracil (U). Due to base pairing, [[i.e. ]{}]{}either the specific interaction of C with G or the interaction of A with U, the secondary structure is formed. Therefore, on an abstract level, the secondary structure is given by the list of all base pairs present in the molecule. Only after the secondary structure has formed, tertiary contacts arise. Pseudoknots [@Richards1969; @Batenburg2000], helix stacking, and base triples [@Higgs2000] as well as the overall three-dimensional arrangement of the molecule are considered as parts of the tertiary structure. The main assumption of hierarchical folding is, that tertiary structure formation operates only on already existing secondary structure elements [@Tinoco1999]. Although cases are known where this approximation breaks down, it generally constitutes a valid starting point [@Cho2009]. In this paper, where the main point is the influence of the loop entropy and the salt concentration on the secondary structure, we therefore neglect tertiary interactions altogether.

![Schematic representation of the secondary structure of an RNA molecule. Dots represent one base, [[i.e. ]{}]{}cytosine, guanine, adenine, or uracil. Solid lines denote the sugar-phosphate backbone bonds, broken lines base pairs, and thick gray lines the non-nested backbone bonds, which are counted by the variable ${M}$, here ${M}=11$. The thick arrows to either side illustrate the force ${F}$ applied to the 5’- and 3’-end.[]{data-label="fig:1"}](049_rna_secondary_structure_scheme_gs_image)

Given a set of base pairs, the secondary structure consists of helices and loops as the basic structural units, [[cf. ]{}]{}[fig.]{} \[fig:1\]. Since pseudoknots are neglected, every nucleotide can be attributed unambiguously to exactly one subunit. The free energy of a certain secondary structure is then given by the sum of the free energy contributions of the individual structural subunits, as we will detail now.

Free energy of a loop {#sec:free-energy-loop}
---------------------

We model the free energy of a loop consisting of $m$ backbone bonds, see [fig.]{} \[fig:1\], with $$\label{eq:1}
  {{\mathcal G}_{\mathrm{l}}}({m})  = {{\mathcal G}_{\mathrm{l}}^{\mathrm{conf}}}({m}) + {{\mathcal G}_{\mathrm{l}}^{\mathrm{salt}}}({m}) + {{\mathcal G}_{\mathrm{l}}^{\mathrm{init}}}{\;}.$$ The first term is the loop entropy difference between an unconstrained polymer and a ring-like polymer, which is characterized by the loop exponent ${c}$ [@Duplantier1986; @Kafri2000; @Einert2008; @Mueller2002] $$\label{eq:2}
  {{\mathcal G}_{\mathrm{l}}^{\mathrm{conf}}}({m}) = - {\mathrm{k_B}T}\ln {m}^{-{c}}
  {\;},$$ with ${\mathrm{k_B}}$ the Boltzmann constant and ${T}$ the absolute temperature. The loop exponent ${c}$ is ${c}_{\mathrm{ideal}}=3/2$ for an ideal polymer and ${c}_{\mathrm{SAW}} = 1.76$ for an isolated self avoiding loop. Helices emerging from the loop limit the configurational space available to the loop and hence increase ${c}$. One obtains ${c}_1=2.06$ for terminal, ${c}_2=2.14$ for internal loops and ${c}_4=2.16$ for a loop with four emerging helices [@Einert2008]. Since the differences between these exponent values are quite small, we assume a constant loop exponent ${c}= 2.1$ in this paper and only compare with the case of vanishing loop entropy characterized by $c=0$.

The second term in [eq.]{}  describes the free energy difference between a charged ring of length $m{l_{\mathrm{ss}}}$ and a straight rod of the same length due to electrostatic interactions, with ${l_{\mathrm{ss}}}= \unit{6.4}{\angstrom}$ the length of one ssRNA backbone bond [@Tan2008]. The electrostatics are modeled on the Debye-Hückel level [@Kunze2002] $$\begin{aligned}
  \label{eq:3}
  {{\mathcal G}_{\mathrm{l}}^{\mathrm{salt}}}({m}) &= {\mathrm{k_B}T}{l_{\mathrm{B}}}({m}{l_{\mathrm{ss}}}) {{\tau}_{\mathrm{ss}}}^2 \Biggl[
  \ln({\kappa}{m}{l_{\mathrm{ss}}}) - \ln(\pi/2) + \gamma
    - \frac{{\kappa}{m}{l_{\mathrm{ss}}}}{2} \leftidx{_1}{F}{_2}\left(
    1/2, \begin{pmatrix} 1\\3/2 \end{pmatrix}, \left(\frac{{\kappa}{m}{l_{\mathrm{ss}}}}{2\pi}\right)^2
  \right)\notag\\
  & {\quad}{\quad}+ \frac12 \left(\frac{{\kappa}{m}{l_{\mathrm{ss}}}}{\pi}\right)^2 \leftidx{_2}{F}{_3}\left(
    \begin{pmatrix}
      1\\1
    \end{pmatrix},
    \begin{pmatrix}
      3/2\\3/2\\2
    \end{pmatrix}, \left(\frac{{\kappa}{m}{l_{\mathrm{ss}}}}{2\pi}\right)^2
  \right) \notag\\ & {\quad}{\quad}+ \frac{1}{{\kappa}{m}{l_{\mathrm{ss}}}}\left(1-\exp(-{\kappa}{m}{l_{\mathrm{ss}}}) + {\kappa}{m}{l_{\mathrm{ss}}}\Gamma(0, {\kappa}{m}{l_{\mathrm{ss}}})\right) \Biggr] {\;},\end{aligned}$$ with ${l_{\mathrm{B}}}= {e_0}^2/({\mathrm{k_B}T}4\pi{\varepsilon_0}{\varepsilon_{\mathrm{r}}}) $ the Bjerrum length, which in water has a value of roughly $\unit{7}{\angstrom}$, ${\kappa}^{-1} =
\sqrt{{\varepsilon_0}{\varepsilon_{\mathrm{r}}}{\mathrm{k_B}T}/(2{\mathrm{N_A}}{e_0}^2I) }$ the Debye screening length, ${\varepsilon_0}$ the vacuum dielectric constant, ${\varepsilon_{\mathrm{r}}}\approx80$ the relative dielectric constant of water [@Murrell1994], ${I}= 1/2({{\rho}_{\mathrm{a}}}{{z}_{\mathrm{a}}}^2 +
{{\rho}_{\mathrm{c}}}{{z}_{\mathrm{c}}}^2)$ the ionic strength, ${{\rho}_{\mathrm{a}}}/{{\rho}_{\mathrm{c}}}$ and ${{z}_{\mathrm{a}}}/{{z}_{\mathrm{c}}}$ the concentration and the valency of the anions/cations, ${\mathrm{N_A}}$ the Avogadro constant, ${e_0}$ the elementary charge, $\gamma\approx0.58$ Euler’s constant, $\Gamma(a,x)$ the incomplete gamma function, and $\leftidx{_p}{F}{_q}$ the generalized hypergeometric functions [@Abramowitz2002]. To account for modifications of the line charge density ${{\tau}_{\mathrm{ss}}}$ due to non-linear electrostatic effects, we employ Manning’s counterion condensation theory [@Manning1969], predicting $${{\tau}_{\mathrm{ss}}}= \min(1/ {l_{\mathrm{ss}}},  1/({l_{\mathrm{B}}}{{z}_{\mathrm{c}}}))\label{eq:4}
  {\;}.$$ [Eq.]{}  amounts to a ground state approximation of the electrostatic contribution to the free energy of a loop. This is rationalized by the fact that the electrostatic interaction is screened and decays exponentially over the Debye length, which is roughly ${\kappa}^{-1} = \unit{1}{\nano\meter}$ for salt solution. However, typical distances between bases in a loop are of the order of the helix diameter $d=\unit{2}{\nano\meter}$ or larger. Therefore, we expect electrostatic interactions to be basically independent of the global configuration of a loop, which justifies both the ground state approximation and our additivity approximation, where ion effects and conformational contributions decouple, see [eq.]{} . In the supporting material, see [eq.]{} S4, we give an interpolation formula for [eq.]{}  involving no hypergeometrical functions.

The last term in [eq.]{}  is the loop initiation free energy ${{\mathcal G}_{\mathrm{l}}^{\mathrm{init}}}$. As we are employing a logarithmic loop entropy, [eq.]{} , we cannot use the standard value for ${{\mathcal G}_{\mathrm{l}}^{\mathrm{init}}}$, which was extracted from experimental data for a different loop parameterization [@Xia1998; @Mathews1999]. Therefore, a modified value ${{\mathcal G}_{\mathrm{l}}^{\mathrm{init}}}$ is obtained by fitting ${{\mathcal G}_{\mathrm{l}}}({m})$, given by [eq.]{} , to experimental data using ${c}= 2.1$ in ${{\mathcal G}_{\mathrm{l}}^{\mathrm{conf}}}({m})$ and the salt concentration ${\rho}=
\unit{1}{\molar}$ in ${{\mathcal G}_{\mathrm{l}}^{\mathrm{salt}}}({m})$, see [fig.]{} \[fig:2\]a. In this figure we show experimentally determined free energies for terminal, internal and bulge loops as a function of the loop size, which exhibit a dependence on the type of the loop. As an approximation, we do not distinguish between those loop types in the theory and consequently fit a single parameter ${{\mathcal G}_{\mathrm{l}}^{\mathrm{init}}}$ to the data, which turns out to be ${{\mathcal G}_{\mathrm{l}}^{\mathrm{init}}}= \unit{1.9}{\kilo\calory\per\mole}$ for ${T}=
\unit{300}{\kelvin}$, see supporting material section C. In [fig.]{} \[fig:2\]a the fitted ${{\mathcal G}_{\mathrm{l}}}({m})$ for the loop exponent ${c}=2.1 $ is depicted by the solid line; the other lines illustrate the effect of different loop exponents on the loop free energy according to [eq.]{}  using the same value for ${{\mathcal G}_{\mathrm{l}}^{\mathrm{init}}}$. [fig.]{} \[fig:2\]b illustrates the effect of salt on the loop free energy for a given value of ${c}=2.1$.

![a Free energy of a loop as a function of the number of segments $m$ for different loop exponents ${c}= 0,\,1.5,\,2.1$ (lines) and for NaCl concentration ${\rho}=\unit{1}{\molar} $. Symbols denote experimental values for various types of loops (hairpin, bulge, internal) [@Serra1995; @Tan2008] for ${\rho}=
\unit{1}{\molar}$ NaCl. ${{\mathcal G}_{\mathrm{l}}^{\mathrm{init}}}$ is obtained by fitting ${{\mathcal G}_{\mathrm{l}}}({m})$, [eq.]{} , to the experimental data for ${c}=
2.1$ and ${\rho}= \unit{1}{\molar}$. The same salt concentration and the same value for ${{\mathcal G}_{\mathrm{l}}^{\mathrm{init}}}$ is used for plotting the curves with ${c}= 0,\,1.5$. b Salt dependence of the free energy of loops as a function of the number of segments $m$ for different salt concentrations ${\rho}= \unit{1}{\molar},\
\unit{0.1}{\molar},\,\unit{0.01}{\molar}$ and loop exponent ${c}=
2.1$ according to [eq.]{} . []{data-label="fig:2"}](043_bsl_combine_029-030_gs_image.pdf)

Free energy of a helix {#sec:free-energy-helix}
----------------------

The free energy of a helix $$\label{eq:5}
  {{\mathcal G}_{\mathrm{h}}}= {{\mathcal G}_{\mathrm{h}}^{\mathrm{stack}}}+ {{\mathcal G}_{\mathrm{h}}^{\mathrm{init}}}+ {{\mathcal G}_{\mathrm{h}}^{\mathrm{term}}}+ {{\mathcal G}_{\mathrm{h}}^{\mathrm{salt}}}{\;}$$ depends on the sequence $\{b_i\}$, which consists of the four nucleotides $b_i =\mathrm{C,G,A,U}$. The stacking free energy ${{\mathcal G}_{\mathrm{h}}^{\mathrm{stack}}}$ is based on experimentally determined parameters incorporating the base pairing free energy as well as the stacking free energy between neighboring base pairs. In the standard notation, ${{{ g}_{\mathrm{h}}}^{\mathrm{stack}}}[{(b_{i},b_{j})},{(b_{i+1},b_{j-1})}]$ is the contribution of the two neighboring, stacked base pairs ${(b_{i},b_{j})}$ and ${(b_{i+1},b_{j-1})}$ to ${{\mathcal G}_{\mathrm{h}}^{\mathrm{stack}}}$. The explicit values for the enthalpic and entropic parts are given in the supporting material. We use the expanded nearest neighbor model [@Mathews1999; @Xia1998] to calculate the base pairing and stacking contributions of a helical section ranging from base pair ${(i,j)}$ through ${(i+h,j-h)}$ and obtain $$\label{eq:6}
  {{\mathcal G}_{\mathrm{h}}^{\mathrm{stack}}}= \sum_{h'=1}^h{{{ g}_{\mathrm{h}}}^{\mathrm{stack}}}[(b_{i+h'-1}, b_{j-h'+1}),(b_{i+h'}, b_{j-h'})]
  {\;}.$$ The initiation and termination free energies in [eq.]{}  take into account weaker pairing energies of AU or GU base pairs at the ends of the helix. We use the standard literature values for ${{\mathcal G}_{\mathrm{h}}^{\mathrm{init}}}$ and ${{\mathcal G}_{\mathrm{h}}^{\mathrm{term}}}$ [@Mathews1999; @Xia1998] and summarize the explicit values in the supporting material. Increasing the salt concentration increases the stability of a helix: First, counterions condense on the negatively charged backbone and reduce the electrostatic repulsion and, second, the diffuse counterion cloud surrounding the charged molecule screens the interaction. We model the two strands of a helix as two parallel rods at distance $d=\unit{2}{\nano\meter}$ interacting [[via]{} ]{}a Debye-Hückel potential characterized by the screening length ${\kappa}^{-1}$. The electrostatic interaction energy per nucleotide with the other strand is given by $$\label{eq:7}
  {{ g}_{\mathrm{h}}^{\mathrm{DH}}}({\rho})= {\mathrm{k_B}T}{{\tau}_{\mathrm{ds}}}^2{l_{\mathrm{ds}}}^{}{l_{\mathrm{B}}}^{}
  \int_{-\infty}^\infty \frac{\exp(-{\kappa}\sqrt{d^2+z^2})}{\sqrt{d^2+z^2}} dz   =
  2 {\mathrm{k_B}T}{{\tau}_{\mathrm{ds}}}^2{l_{\mathrm{ds}}}^{}{l_{\mathrm{B}}}^{} {{\mathrm{K}}_{0}({\kappa}d)}
  {\;}.$$ ${l_{\mathrm{ds}}}= \unit{3.4}{\angstrom}$ is the helical rise per base pair of double-stranded RNA (dsRNA) and ${{\mathrm{K}}_{0}({\kappa}d)}$ is the zeroth order modified Bessel function of the second kind. Again, we employ Manning’s theory [@Manning1969] to calculate the line charge density ${{\tau}_{\mathrm{ds}}}= \min(1/  {l_{\mathrm{ds}}}, 1/({l_{\mathrm{B}}}{{z}_{\mathrm{c}}}))$. The reference state for the salt correction of the pairing free energy is at temperature ${T}= \unit{300}{\kelvin}$ with monovalent salt concentration ${\rho}= \unit1\molar$, as the experimental pairing free energies ${{{ g}_{\mathrm{h}}}^{\mathrm{stack}}}$ were determined at this concentration. The free energy shift for a helix consisting of $h$ base pairs due to electrostatic interactions is then $$\label{eq:8}
  {{\mathcal G}_{\mathrm{h}}^{\mathrm{salt}}}= h( {{ g}_{\mathrm{h}}^{\mathrm{DH}}}({\rho}) -   {{ g}_{\mathrm{h}}^{\mathrm{DH}}}(\unit{1}{\molar}))
  {\;}.$$

The use of Debye-Hückel theory to incorporate salt effects enables to include the overall dependence on temperature and salt concentration but involves several approximations. First, we are using Manning’s counterion condensation theory to obtain the actual line charge density of ssRNA and dsRNA [@Manning1969]. However, Manning condensation is known to underestimate the line charge at increasing salt concentration and therefore favors the bound state [@Netz2003]. Second, when calculating the electrostatic energy of a loop we effectively use a ground state approximation and neglect conformational fluctuation effects. Third, when two ssRNA strands come together to form a helix, the line charge density increases since the distance between two bases decreases. The salt dependence of the work to decrease the axial distance between two bases from ${l_{\mathrm{ss}}}= \unit{6.4}{\angstrom}$ to ${l_{\mathrm{ds}}}=
\unit{3.4}{\angstrom}$ is neglected. This approximation favors the unbound state. Therefore, it is very important to validate the model we employ, which we do by detailed comparison with experimental data. From the favorable comparison with experiments we tentatively conclude that the various errors partially cancel and the resulting expression for the salt influence is quite accurate. We point out that after determining ${{\mathcal G}_{\mathrm{l}}^{\mathrm{init}}}$ in [eq.]{} , no further fitting is done and only standard literature values are used.

Our theory is able to consider variations of the salt concentration as well as of the temperature, which makes it suitable to study RNA melting at various salt concentrations in a consistent way. However, since our approach is solely based on mean field theory, it will become unreliable in the case of multivalent ions, where correlations become important. Also, ion specific effects, which are important for divalent ions such as Mg^2+^ [@Draper2005], are not considered in our approach.

Response of the molecule to an external stretching force {#sec:response-molecule-an}
--------------------------------------------------------

In atomic force microscope or optical tweezers experiments, it is possible to apply a stretching force ${F}$ to the two terminal bases of the molecule. We model the stretching response of the ${M}$ non-nested backbone bonds, see [fig.]{} \[fig:1\], with the freely jointed chain (FJC) model [@Montanari2001; @Gerland2001; @Gerland2003]. A non-nested bond is defined as a backbone bond, which is neither part of a helix nor part of a loop. It is outside all secondary structure elements and therefore contributes to the end-to-end extension observed in force spectroscopy experiments. The force dependent contribution to the free energy per non-nested monomer is given by $$\label{eq:9}
  {{{ g}_{\mathrm{f}}}^{\mathrm{FJC}}}= {{{\mathcal G}_{\mathrm{f}}}^{\mathrm{FJC}}}/ {M}= -{\mathrm{k_B}T}\frac{{l_{\mathrm{ss}}}}{{b_{\mathrm{ss}}}} \ln
  \left(\frac{\sinh(\beta{b_{\mathrm{ss}}}{F})}{\beta{b_{\mathrm{ss}}}{F}}\right)
  {\;},$$ where $\beta = 1/({\mathrm{k_B}T})$ is the inverse thermal energy and ${b_{\mathrm{ss}}}=
\unit{1.9}{\nano\meter}$ is the Kuhn length of ssRNA [@Montanari2001] (we used the Kuhn length of ssDNA as the corresponding ssRNA data is less certain). The stretching response of one non-nested monomer to an external force is then given by $$\label{eq:10}
  x^{\mathrm{FJC}}({F}) = -\frac{{\mathrm{d}}{{{ g}_{\mathrm{f}}}^{\mathrm{FJC}}}}{{\mathrm{d}}{F}}  = {l_{\mathrm{ss}}}{\mathcal L}(\beta{F}{b_{\mathrm{ss}}}) =
  {l_{\mathrm{ss}}}\left( \coth(\beta{F}{b_{\mathrm{ss}}}) + 1/ (\beta{F}{b_{\mathrm{ss}}})\right)
  {\;},$$ ${\mathcal L}$ is the Langevin function. Electrostatic effects on the stretching response are considered to be small and hence are neglected [@Marko1995; @Netz2001a].

Calculation of the partition function {#sec:calc-part-funct}
=====================================

So far we showed how to calculate the free energy of one given secondary structure. The next step is to enumerate all possible secondary structures and to obtain the partition function, which allows to study the thermodynamics of the system. As we neglect tertiary contacts – and in particular pseudoknots – for any two base pairs ${(i,j)}$ and ${(k,l)}$ with $i<j$, $k<l$, and $i<k$ we have either $i<k<l<j$ or $i<j<k<l$. This allows to derive a recursion relation for the partition function of the secondary structure. In our notation, the canonical partition function ${Q}_{i,j}^{{M}}$ of a sub-strand from base $i$ at the 5’-end through $j$ at the 3’-end depends on the number of non-nested backbone bonds ${M}$ [@Einert2008; @Bundschuh2005; @Mueller2002], see [fig.]{} \[fig:1\]. The recursion relations for ${Q}_{i,j}^{{M}}$ can be written as

\[eq:11\] $$\label{eq:9a}
    {Q}_{i,j+1}^{{M}+1} = \frac{{v_{\mathrm{f}}({M}+1)}}{{v_{\mathrm{f}}({M})}} \left[ {Q}_{i,j}^{{M}} +
      \sum_{k=i+{M}+1}^{j-N_{\mathrm{loop}}}{Q}_{i,k-1}^{{M}} {Q}_{k, j+1}^{0}\right]$$ and $$\label{eq:9b}
    {Q}_{k,j+1}^{0}= \sum_{h=1}^{(j-k-N_{\mathrm{loop}})/2}
    \exp[-\beta{{{\mathcal G}_{\mathrm{h}}}}^{{(k,j+1)}}_{{(k+h,j+1-h)}}]\sum_{m=1}^{j-k-1-2h}
    {Q}_{k+1+h,j-h}^{m}\frac{\exp[-\beta{{\mathcal G}_{\mathrm{l}}}(m+2)]}{{v_{\mathrm{f}}(m)}}
    {\;}.$$

[Eq.]{}  describes elongation of an RNA structure by either adding an unpaired base (first term) or by adding an arbitrary sub-strand $ {Q}_{k, j+1}^{0}$ that is terminated by a helix. [Eq.]{}  constructs $ {Q}_{k, j+1}^{0}$ by closing structures with ${m}$ non-nested bonds, summed up in $
{Q}_{k+1+h,j-h}^{{m}}$, by a helix of length $h$. $N_{\mathrm{loop}}=3$ is the minimum number of bases in a terminal loop. ${v_{\mathrm{f}}({M})}$ denotes the number of configurations of a free chain with ${M}$ links and drops out by introducing the rescaled partition function ${\tilde Q}_{i,j}^{{M}}={Q}_{i,j}^{{M}}/{v_{\mathrm{f}}({M})}$ and will not be considered further since its effects on the partition function are negligible. ${{{\mathcal G}_{\mathrm{h}}}}^{{(k,j+1)}}_{{(k+h,j+1-h)}}$ is the free energy of a helix beginning with base pair ${(k,j+1)}$ and ending with base pair ${(k+h,j+1-h)}$ according to [eq.]{} . ${{\mathcal G}_{\mathrm{l}}}(m+2)$ is the free energy of a loop consisting of $m+2$ segments as given by [eq.]{} . ${{\mathcal G}_{\mathrm{l}}}$ and ${{\mathcal G}_{\mathrm{h}}}$ contain all interactions discussed in the previous section. [Eq.]{}  allows to compute the partition function in polynomial time (${\mathcal O}({N}^4)$). Further, our formulation allows to treat non-linear functions for ${{\mathcal G}_{\mathrm{l}}}({m})$ and ${{\mathcal G}_{\mathrm{h}}}(h)$; for instance, ${{\mathcal G}_{\mathrm{l}}}({m})$ is strongly non-linear by virtue of [eqs.]{}  and .

The unrestricted partition function of the entire RNA, where the number of non-nested backbone bonds ${M}$ is allowed to fluctuate, is given by $${Z}_{{N}} = \sum_{{M}= 0}^N \exp[-\beta {{{ g}_{\mathrm{f}}}^{\mathrm{FJC}}}{M}]{\tilde Q}_{0,N}^{{M}}
  \label{eq:12}$$ and contains the influence of force [[via]{} ]{}${{{ g}_{\mathrm{f}}}^{\mathrm{FJC}}}$ defined in [eq.]{} . The partition function ${Z}_{{N}}$ contains all secondary structure interactions, but neglects pseudoknots and other tertiary interactions. As has been argued before, this approximation is known to work very well [@Tinoco1999] and yields reliable predictions for the stability of nucleic acids [@Gruber2008].

Using the same ideas, we determine the minimum free energy (mfe) and the mfe structure. The mfe structure, is defined as the secondary structure, which gives the largest contribution to the partition function. Since it cannot be derived from the partition function itself, it has to be determined from a slightly modified set of recursion relations, see supporting material.

Salt dependence of melting curves {#sec:salt-depend-melt}
=================================

![ Melting curve of the 76 bases long tRNA-phe of yeast; the minimum free energy structure at ${\rho}= \unit{1}{\molar}$, ${T}= \unit{300}{\kelvin}$, ${c}= 2.1$ is shown as an inset. Symbols denote experimental melting curves for NaCl concentrations ${\rho}=
    \unit{20}{\milli\molar}$ (squares) and $\unit{150}{\milli\molar}$ (circles) [@Privalov1978]. Our predictions for different salt concentrations are depicted by the dashed ($\unit{20}{\milli\molar}$), dash-dotted ($\unit{150}{\milli\molar}$), and solid ($\unit{1}{\molar}$) lines. The respective arrows indicate melting temperatures obtained by experiments of another group [@Vives2002] at the same salt concentration ${\rho}=\unit{150}{\milli\molar}$ (left arrow) and ${\rho}=\unit{1}{\molar}$ (right arrow). The dotted line shows our prediction for ${\rho}= \unit{150}{\milli\molar}$ and ${c}= 0$ and exemplifies that a non-zero loop exponent is responsible for rendering the transition more cooperative, in closer agreement with experiment; for $\unit{150}{\milli\molar}$ and ${c}= 0$ the melting temperature is at higher temperatures since the energy parameters are optimized for ${c}= 2.1$. The gray dash-dotted curve is the prediction of the Vienna package, which uses a linearized multi-loop entropy corresponding to ${c}= 0$ and ${\rho}=
\unit{1}{\molar}$. This is to be compared to our prediction for ${c}= 2.1$ and ${\rho}= \unit{1}{\molar}$: while the melting temperatures are similar, the cooperativity, [[i.e. ]{}]{}the widths of the peaks are different due to different loop exponents. []{data-label="fig:3"}](031_rna_salt_melting_curve_tRNA-phe_Privalov_gs_image)

In this section we calculate melting curves for different salt concentrations by applying [eqs.]{}  and , which include our salt dependent free energy parameterization. In [fig.]{} \[fig:3\] we compare experimental results [@Privalov1978; @Vives2002] with our predictions for the heat capacity of yeast tRNA-phe; the sequence is given in the supporting material section D. The heat capacity is readily obtained by $${C}= {T}\frac{\partial^2{\mathrm{k_B}T}\ln {Z}_N}{\partial{T}^2}
  {\;},\label{eq:13}$$ where ${Z}_{N}$ is the unrestricted partition function of the RNA at zero force, [eq.]{} . In all our calculations, we use the same literature parameter set for the stacking and pairing free energy ${{{ g}_{\mathrm{h}}}^{\mathrm{stack}}}$. No additional fit parameter enters except the loop initialization free energy ${{\mathcal G}_{\mathrm{l}}^{\mathrm{init}}}$, which is determined in [fig.]{} \[fig:2\]a from a separate experimental data set. The salt dependence of the experimentally observed melting temperatures is reproduced well, compare [fig.]{} \[fig:3\]. The arrows indicate additional experimental results [@Vives2002] for the melting temperature for ${\rho}= \unit{150}{\milli\molar}$ and $\unit{1}{\molar}$, which again coincide with our prediction. We also plot a calculated melting curve for loop exponent ${c}= 0$ and NaCl concentration ${\rho}= \unit{150}{\milli\molar}$, which exhibits a far less cooperative transition than observed in the corresponding curve with ${c}= 2.1$. Finally, we compare our prediction for ${\rho}=
\unit{1}{\molar}$ and ${c}= 2.1$ with the prediction of RNAheat in the Vienna Package [@Hofacker1994] for ${\rho}= {1}{\molar}$, which uses a linearized multi-loop entropy amounting to ${c}= 0$ in our framework. The predicted melting temperatures are almost identical. However, the widths of the peaks in both melting curves differ and our melting profile for ${c}= 2.1$ is more peaked. Taking all these observations together leads to the conclusion that only a combined use of logarithmic loop entropy (characterized by a non-zero loop exponent) and salt dependent free energy corrections leads to a correct prediction of melting curves. The additional features in the experimental data, [[e.g. ]{}]{}the shoulder at lower temperatures and the increased width of the experimental curves might be attributed to tertiary structure rearrangements, which are not captured by our approach, or to melting occurring in multiple stages.

Salt dependence of stretching curves {#sec:salt-depend-stretch}
====================================

Apart from temperature, force is an important variable to study denaturation of RNA molecules [@Li2008; @Liphardt2001; @Manosas2006; @Seol2007; @Tinoco2004; @Woodside2006; @Woodside2006a; @Hyeon2005; @Gerland2001; @Kumar2008a; @Mueller2002; @Montanari2001]. In [fig.]{} \[fig:4\] we show the salt dependence of stretching curves for yeast tRNA-phe. The stretching curves have been obtained by describing the force response of the ${M}$ non-nested backbone bonds, see [fig.]{} \[fig:1\], with the freely jointed chain (FJC) model, see [eq.]{} , $$\label{eq:14}
  x({F}) =  {\mathrm{k_B}T}\frac{\partial \ln {Z}_N}{\partial{F}} = {\mathrm{k_B}T}\frac{\partial\ln{Z}_N}{\partial{{{ g}_{\mathrm{f}}}^{\mathrm{FJC}}}}\frac{\partial{{{ g}_{\mathrm{f}}}^{\mathrm{FJC}}}}{\partial{F}} = M x^{\mathrm{FJC}}({F})
  {\;},$$ where we used the expectation value of the number of non-nested backbone segments $$\label{eq:15}
  {M}= -{\mathrm{k_B}T}\frac{\partial\ln{Z}_N}{\partial{{{ g}_{\mathrm{f}}}^{\mathrm{FJC}}}}
  {\;}.$$

![Salt dependence of stretching curves of tRNA-phe for different salt concentrations ${\rho}=
  \unit{20}{\milli\molar},\,\unit{150}{\milli\molar},\,\unit{1}{\molar}$. Increasing salt concentration stabilizes the secondary structure due to screening of the electrostatic interaction. The dotted line is the theoretical prediction for the force extension curve of a freely jointed chain, [eq.]{} . The deviation of the predicted curves for RNA from the FJC curve is due to the presence of secondary structure. The observed plateau force is due to the rupture of the secondary structure. We show the force extension curves in a a linear plot and b a double logarithmic plot, indicating that the force extension curve is linear in the low-force regime, before the secondary structure is ruptured apart.[]{data-label="fig:4"}](044_rna_combine_34-35_gs_image)

As for the melting curves, one observes that increasing salt concentration stabilizes the structure, leading to higher unfolding forces. All curves converge in the large force limit to a freely jointed chain of the length of the whole RNA molecule $(N-1){l_{\mathrm{ss}}}$, where $N = 76$ is the number of bases in the chain. The deviation for small forces from this theoretical prediction is due to the secondary structure of RNA, which is present at small forces and which becomes disrupted at forces ${F}\gtrsim \unit{\text{3-7}}{\pico\newton}$. In [fig.]{} \[fig:5\]a we show the force extension curve of the P5ab hairpin [@Liphardt2001]; the sequence is given in the supporting material section D. Apart from the salt dependence of the force extension curve, one observes that the unzipping of the helix occurs in two stages. This is seen best by considering the fraction of non-nested segments and its derivative, [fig.]{} \[fig:5\]b. The first stage is a smooth unzipping of the first three base pairs up to the bulge loop visible as a shoulder at ${F}\approx{8}{\pico\newton}$ in the derivative. The second stage is a sharp transition, where the rest of the hairpin unzips. In [fig.]{} \[fig:5\]c we show mfe predictions for the secondary structure at different forces for ${\rho}= \unit{1}{\molar}$ NaCl. For ${F}<
\unit{5}{\pico\newton}$, we predict correctly the experimentally observed native state with all base pairs intact [@Liphardt2001]. For forces ${F}\approx\unit{8}{\pico\newton}$, an intermediate state appears, where the first three base pairs are unzipped up to the bulge loop. Denaturation is observed for ${F}\gtrsim
\unit{14}{\pico\newton}$. The native structure of the P5ab hairpin contains the stacked pairs ${{\mathrm{GG}}\atop{\mathrm{AA}}}$ – bp[(17,42)]{} and bp[(18,41)]{} [@Liphardt2001]. For this stack, no free energy parameters are available and we use the parameters for the stack ${{\mathrm{GG}}\atop{\mathrm{UU}}}$, instead. However, other parameterizations for this stack work equally well and reproduce the experimental transition force within errors, see [fig.]{} \[fig:9\].

![a Salt dependence of stretching curves of the 56 bases long RNA hairpin P5ab [@Liphardt2001] for different salt concentrations ${\rho}=
  \unit{20}{\milli\molar},\,\unit{150}{\milli\molar},\,\unit{1}{\molar}$. Increasing salt concentration stabilizes the secondary structure due to screening of the electrostatic interaction. The dotted line is the force extension curve of an FJC, [eq.]{} . b The fraction of non-nested segments ${M}/{N}$ as a function of force. One observes that the unzipping of the hairpin occurs in two stages, which is visible as a shoulder for $\unit{5}{\pico\newton}\lesssim{F}\lesssim
  \unit{6-10}{\pico\newton}$ (exact values depend on the salt concentration) and a successive cooperative transition. The inset shows the derivative ${\mathrm{d}}{M}/({N}{\mathrm{d}}{F})$ for ${\rho}=
  \unit{1}{\molar}$, where the first transition is visible as a shoulder at ${F}\approx\unit{8}{\pico\newton}$. The sharp peak at ${F}=\unit{11}{\pico\newton}$ is the rupture of the complete helix. c Predicted minimum free energy structures of the hairpin P5ab at different forces, see also the supporting material section A. For ${F}< \unit{5}{\pico\newton}$ the hairpin is in the native state with all base pairs intact. At ${F}\approx\unit{8}{\pico\newton}$ the first helix, consisting of three base pairs and bounded by the bulge loop, is ruptured. This causes the first smooth transition. Forces ${F}\gtrsim \unit{14}{\pico\newton}$ lead to the unzipping of the whole hairpin in a very cooperative fashion.[]{data-label="fig:5"}](045_rna_combine_32-33-36_gs_image)

![The effect of different parameterizations for the free energy parameters for ${{\mathrm{GG}}\atop{\mathrm{AA}}}$ on the denaturation curve is marginal. Here, the fraction of non-nested backbone bonds is plotted against the force for the P5ab hairpin and ${T}= \unit{298}{\kelvin}$, ${\rho}= {250}{\milli\molar}$, ${c}= 2.1$. The solid line is obtained with the parameters used for the P5ab hairpin in the rest of this paper ${{{ g}_{\mathrm{h}}}^{\mathrm{stack}}}({{\mathrm{GG}}\atop{\mathrm{AA}}}) =
  {{{ g}_{\mathrm{h}}}^{\mathrm{stack}}}({{\mathrm{GG}}\atop{\mathrm{UU}}})$. The dotted line is obtained by using ${{{ g}_{\mathrm{h}}}^{\mathrm{stack}}}({{\mathrm{GG}}\atop{\mathrm{AA}}}) =
  {{{ g}_{\mathrm{h}}}^{\mathrm{stack}}}({{\mathrm{UU}}\atop{\mathrm{AA}}})$, whereas the dashed line is for ${{{ g}_{\mathrm{h}}}^{\mathrm{stack}}}({{\mathrm{GG}}\atop{\mathrm{AA}}}) =
  0$. All three curves coincide and differ only slightly at the transition, exhibiting only marginally different transition forces, which all agree with the experimentally observed unfolding force within errors@Liphardt2001. The values of the free energy parameters are given in the supporting material.[]{data-label="fig:9"}](123_rna_salt_compare_different_parameterizations_for_GAGA_in_P5ab_gs_image)

Phase diagrams of RNA hairpin P5ab
==================================

With the tools established in the previous sections, we are now able to study phase diagrams of RNA. We consider the P5ab hairpin, which is a well studied system [@Liphardt2001; @Gerland2003; @Hyeon2005; @Wen2007; @Cocco2003]. In [fig.]{} \[fig:6\]b the phase diagram in the ${F}$-${\rho}$ plane is shown for ${T}=
\unit{298}{\kelvin},\,\unit{300}{\kelvin},\,\unit{320}{\kelvin}$ and ${c}= 2.1$. The phase boundary is defined as the force where half of the helical section is unzipped. For the definition of the phase boundary, we exclude the three unpaired bases at the 5’- and the four bases at the 3’-end, see [fig.]{} \[fig:5\]c, and use the condition ${M}- 7 = ({N}-7)/2$. This threshold value of ${M}/{N}$ is depicted by an arrow in [fig.]{} \[fig:6\]a. Below the phase boundary, the hairpin is stable, above the molecule is denatured. In [fig.]{} \[fig:8\]b we additionally include the experimental results by @Liphardt2001 agreeing nicely with our results. It is important to note, that this transition is not a phase transition in the strict statistical mechanics sense, but just a crossover. A true phase transition is defined as a non-analyticity of the free energy, which can only occur for an infinite system with long-range interactions [@Einert2008]. The three-dimensional phase space we are considering is spanned by temperature, force, and salt concentration. In [figs.]{} \[fig:7\] and \[fig:8\] we show slices in the ${F}$-${T}$ and in the ${T}$-${\rho}$ plane. The phase boundary for the ${F}$-${T}$ plane is determined the same way as in the ${F}$-${\rho}$ plane, yet with varying temperature and fixed salt concentration. The phase boundary in the ${T}$-${\rho}$ plane is determined differently: heat capacity curves as a function of temperature are calculated for different salt concentrations. The position of the peaks in the heat capacity curves (one is depicted by an arrow in [fig.]{} \[fig:8\]a) determine the phase diagram in [fig.]{} \[fig:8\]b. Therefore, slight differences between the phase diagrams in [figs.]{} \[fig:6\], \[fig:7\] on the one hand and [fig.]{} \[fig:8\] on the other hand may arise.

We observe that for large salt concentrations, the denaturation forces and temperatures are rather independent of the salt concentration, see [figs.]{} \[fig:6\] and \[fig:8\]. Only when the Debye screening length ${\kappa}^{-1}$ is of the order of the typical length scale of RNA, which is the case for ${\rho}\lesssim\unit{100}{\milli\molar}$, a marked dependence on the salt concentration is observed.

![a Fraction of non-nested segments of the P5ab hairpin as a function of force for different salt concentrations and constant temperature ${T}=
  \unit{298}{\kelvin}$. The position of the crossover, which is defined as the point where ${M}- 7 = ({N}-7)/2$, [[i.e. ]{}]{} ${M}/{N}=0.56$ (indicated by the arrow), determines the phase diagram. b Phase diagram of the P5ab hairpin in the ${F}$-${\rho}$ plane for different temperatures ${T}=
  \unit{298}{\kelvin},\,\unit{300}{\kelvin},\,\unit{320}{\kelvin}$. Below the curve the RNA is in the hairpin phase, above the RNA is denatured. The symbol at ${\rho}= \unit{250}{\milli\molar}$, ${F}=\unit{13.3}{\pico\newton}$, and ${T}= \unit{298}{\kelvin}$ denotes the experimental data by @Liphardt2001 and coincides with our prediction.[]{data-label="fig:6"}](046_rna_combine_37-38_gs_image)

![a Fraction of non-nested segments of the P5ab hairpin as a function of force for different temperatures and constant salt concentration ${\rho}= \unit{250}{\milli\molar}$. The position of the crossover (arrow, ${M}/{N}= 0.56$) determines the phase boundary. With increasing temperature a decrease of the denaturation force is observed. Above the melting temperature ${{T}_{\mathrm{m}}}\approx\unit{358}{\kelvin}$ the molecule is always in the denatured state. b Phase diagram of the P5ab hairpin in the ${F}$-${T}$ plane. Below the curve the RNA is in the native hairpin phase, above the RNA is denatured. The symbol denotes experimental values [@Liphardt2001].[]{data-label="fig:7"}](047_rna_combine_39-40_gs_image)

![a Heat capacity curves for different salt concentrations as a function of temperature. The peak position moves to higher temperatures with increasing salt concentration. The positions of the peaks, denoted exemplarily for the curve with ${\rho}=\unit{20}{\milli\molar}$ by the arrow, determine the phase diagram. b Phase diagram of the P5ab hairpin in the ${T}$-${\rho}$ plane. Below the curve the RNA is in the native hairpin phase, above the RNA is denatured.[]{data-label="fig:8"}](048_rna_combine_41-42_gs_image)

Conclusions {#sec:conclusions}
===========

We construct a theory for RNA folding and melting that includes the effects of monovalent salt, loop entropy, and stretching forces. Our theory is based on salt and temperature dependent modifications of the free energies of RNA helices and loops that include electrostatic interactions on the linear Debye-Hückel level – augmented by Manning condensation – and conformational fluctuation effects [[via]{} ]{}the asymptotic, non-linear expression for the entropy of loop formation. Decreasing salt concentration is shown to generally destabilize RNA folds and to lower denaturation temperatures and forces. The predictions are in good agreement with experimental data as shown for two different scenarios, namely the heat capacity curves for the thermal denaturation of tRNA-phe and the response of the P5ab RNA hairpin to an external pulling force.

Due to the usage of the linear Debye-Hückel approximation in conjunction with the Manning condensation concept, our approach is limited to monovalent salt and neglects ion-specific effects. Electrostatic nonlinear and correlation effects could in principle be taken into account by more advanced modeling using variational approaches [@Netz2003], while ion-specific effects could be straightforwardly included using effective interactions between different ions and RNA bases [@Schwierz2010]. More complex phenomena involving multivalent ions such as Mg^2+^ could in principle be modeled by allowing for a few tertiary contacts, which is left for future studies.

We find that for a proper description of RNA melting curves, correct modeling of the loop entropy is crucial. A non-zero loop exponent leads to an increased cooperativity of the melting transition and thus makes the heat capacity curve narrower in good agreement with experimental results. We conclude that for a correct description of RNA denaturation thermodynamics, both loop entropy and salt effects are important and should be included in standard structure and melting curve prediction software.

Acknowledgements {#sec:acknowledgements}
================

Financial support comes from the DFG [[via]{} ]{}grant NE 810/7. T.R.E. acknowledges support from the Elitenetzwerk Bayern within the framework of CompInt.
