Secondary structure formation of homopolymeric single-stranded nucleic acids including force and loop entropy: implications for DNA hybridization

Secondary structure formation of homopolymeric single-stranded nucleic acids including force and loop entropy: implications for DNA hybridization

Thomas R. Einert 1Physik Department, Technische Universität München, James-Franck-Straße, 85748 Garching, Germany, Tel.: +49-89-28914337, Fax: +49-89-28914642, 1einert@ph.tum.de    Henri Orland 2Institut de Physique Théorique, CEA Saclay, 91191 Gif-sur-Yvette Cedex, France 2    Roland R. Netz 1Physik Department, Technische Universität München, James-Franck-Straße, 85748 Garching, Germany, Tel.: +49-89-28914337, Fax: +49-89-28914642, 1einert@ph.tum.de 3Fachbereich Physik, Freie Universität Berlin, Arnimallee 14, 14195 Berlin, Germany 3
August 28, 2019
Abstract

Loops are essential secondary structure elements in folded DNA and RNA molecules and proliferate close to the melting transition. Using a theory for nucleic acid secondary structures that accounts for the logarithmic entropy  for a loop of length , we study homopolymeric single-stranded nucleic acid chains under external force and varying temperature. In the thermodynamic limit of a long strand, the chain displays a phase transition between a low temperature / low force compact (folded) structure and a high temperature / high force molten (unfolded) structure. The influence of on phase diagrams, critical exponents, melting, and force extension curves is derived analytically. For vanishing pulling force, only for the limited range of loop exponents a melting transition is possible; for the chain is always in the folded phase and for always in the unfolded phase. A force induced melting transition with singular behavior is possible for all loop exponents and can be observed experimentally by single molecule force spectroscopy. These findings have implications for the hybridization or denaturation of double stranded nucleic acids. The Poland-Scheraga model for nucleic acid duplex melting does not allow base pairing between nucleotides on the same strand in denatured regions of the double strand. If the sequence allows these intra-strand base pairs, we show that for a realistic loop exponent pronounced secondary structures appear inside the single strands. This leads to a lower melting temperature of the duplex than predicted by the Poland-Scheraga model. Further, these secondary structures renormalize the effective loop exponent , which characterizes the weight of a denatured region of the double strand, and thus affect universal aspects of the duplex melting transition.

Keywords:
RNA, DNA, nucleic acids, denaturation, pulling, melting, loop entropy, phase transition
\addunit\molar

M \addunit\calorycal

1 Introduction

Ribonucleic acids continue to stay in the focus of experimentalists and theorists Gesteland2005 . The advance of single molecule techniques Smith1996 ; Liphardt2002 ; Rief1999 ; Maier2000 ; Bockelmann1997 ; Mossa2009 nowadays allows to study single chains of nucleic acids under tension and varying solution conditions and thereby yields unprecedented insights into the behavior and folding properties of these essential molecules. Theory on RNA folding vastly relies on the idea of hierarchical folding proposed by Tinoco et al. Tinoco1971 ; Tinoco1999 , stating that given a sequence (the primary structure), the secondary structure (i. e. the list of all base pairs) forms independently of the tertiary structure (the overall three-dimensional arrangement of all atoms). This is in contrast to the protein folding problem, which does not feature these well separated energy scales between the different structural levels and hence is more involved Finkelstein2004 . The idea of hierarchical folding therefore suggests to solely focus on the secondary structure, i. e. the base pairs, and to neglect the tertiary structure. This constitutes a major simplification and enables to calculate partition functions exactly and to predict the secondary structure formed by a given RNA sequence. De Gennes Gennes1968 was the first to calculate the partition function of an ideal homopolymeric RNA chain by using a propagator formalism and solving the partition function by means of a singularity analysis of the generating functions. Due to his real space approach for an ideal polymer, the loop exponent was fixed at . The loop exponent  characterizes the logarithmic entropy contribution of a loop of length . Ten years later, Waterman and Smith Waterman1978a devised a recursion relation appropriate for the partition function of folded RNA, which now lies at the heart of most RNA secondary structure and free energy prediction algorithms currently used. Since all results obtained for RNA are also valid for DNA on our relative primitive level of modeling, we will mostly explicitly refer to RNA in our paper but note that in principle all our results carry over to single stranded DNA molecules, as well. Subsequently, several theoretical models were developed to study RNA and DNA: those were focused on melting Einert2008 ; Hofacker1994 ; McCaskill1990 , stretching Einert2010 ; Montanari2001 ; Gerland2001 ; Mueller2002 ; Hanke2008 ; Cocco2002a , unzipping Lubensky2002 , translocation Bundschuh2005 , salt influence Einert2010a ; Tan2006 ; Mamasakhlisov2007 , pseudoknots Baiesi2003 ; Orland2002 , and the influence of the loop exponent Einert2008 ; Einert2010 ; Blossey2003 ; Kafri2000 . In this context, an interesting question in connection with the melting of double stranded DNA (dsDNA) arises: Do secondary structure elements form in the single strands inside denatured dsDNA loops or not? Formation of such secondary structures in dsDNA loops would mean that inter-strand base pairing between the two strands – being responsible for the assembly of the double helix – is in competition with intra-strand pairing, where bases of the same strand interact. This question is not only important for the thermal melting of dsDNA but also for DNA transcription, DNA replication and the force-induced overstretching transition of DNA Alberts2002 ; Rief1999 .

In this paper, the influence of the loop exponent on the behavior of RNA subject to varying temperature and external force is studied, which goes beyond our previous work where only the temperature influence was considered Einert2008 . We neglect sequence effects and consider a long homopolymeric, single stranded RNA molecule. A closed form expression for the partition function is derived, which allows to study the thermodynamic behavior in detail. The phase diagram in the force-temperature plane is obtained. We find that the existence of a temperature induced phase transition depends crucially on the value of the loop exponent : at vanishing force a melting transition is possible only for the limited range of loop exponents . is a typical exponent that characterizes the entropy of loops usually encountered in RNA structures – hairpin loops, internal loops, multi-loops with three or more emerging helices. That means that RNA molecules are expected to experience a transition between a folded and an unfolded state. This is relevant for structure formation in DNA or RNA single strands and can in principle be tested in double laser trap force clamp experiments Gebhardt2010 . Our findings also have implications for the denaturation of double helical nucleotides (e. g. dsDNA). Since intra-strand and inter-strand base pairing compete, secondary structure formation of the single strands inside denatured regions of the duplex has to be taken into account in a complete theory of dsDNA melting. In the case where intra- and inter-strand base pairing occurs, the classical Poland-Scheraga mechanism for the melting of a DNA duplex has to be augmented by the single-strand folding scenario considered by us, as the Poland-Scheraga theory is only valid in the case where no intra-strand base pairing is possible. If the intra-strand interactions are strong enough to induce folded secondary structures, we show that the loop exponent governing the entropy of inter-strand loops is renormalized and takes on an effective universal value that only depends on whether the inter-strand loops are symmetric (consisting of two strands of the same length) or asymmetric. The resulting duplex melting transition is universal and turns out to be strongly discontinuous in the first, symmetric case, and on the border between continuous and discontinuous in the second, asymmetric case. In the case when intra-strand base pairing is weak and a single strand is in the unfolded phase, the situation is different and qualitatively similar to the original Poland-Scheraga results, yet with a lower melting temperature. All these effects can be studied experimentally. We make explicit suggestions for dsDNA sequences, with which the formation of intra-strand secondary structures inside inter-strand loops can be selectively inhibited or favored.

2 Derivation of partition function

2.1 Model

Single stranded RNA is modeled as a one-dimensional chain. The basic units are nucleotides with the four bases (cytosine (C), adenine (A), guanine (G), uracil (U)), which are enumerated by an index . A nucleotide can establish a base pair (bp) with another nucleotide via hydrogen bonding, leading to helices and loops as the structural building blocks of an RNA secondary structure, see fig. 1. In agreement with previous treatments, a valid secondary structure is a list of all base pairs with the constraint that a base can be part of at most one pair. In addition, pseudoknots are not allowed, that means that for any two base pairs and with , , and we have either or  Orland2002 . This imposes a hierarchical order on the base pairs, meaning that two base pairs are either nested and part of the same substructure or are independent and part of different substructures, fig. 2. Helix stacking – the interaction of two helices emerging from the same loop – and base triples as well as the overall three-dimensional structure are not considered. Base pairs are stabilized by two different interactions. First, by hydrogen bonds between complementary bases and, second, by the stacking interaction between neighboring base pairs, which are accounted for by the sequence independent parameters and , respectively. A helix with base pairs consequently has the free energy , where is a helix initiation free energy. Therefore, the hydrogen bonding and the stacking interaction can be combined to yield the binding energy per base pair , which we define to be positive. can be measured experimentally by duplex hybridization Xia1998 and contains the binding free energies as well as the extensive part of configurational polymer entropy. Further, the stacking interaction appears as an additional contribution to the helix initiation free energy. We define and describe the binding free energy by a single, sequence independent parameter Einert2008

(1)

is the statistical weight of a bound base pair, is the absolute temperature, and the Boltzmann constant. Therefore this is a model for homopolymeric RNA, which can be realized experimentally with synthetic alternating sequences or . It has also been argued that this homopolymeric model describes random RNA above the glass transition Bundschuh2002a .

The non-extensive contribution of the free energy of a loop is given by

(2)

and describes the entropy difference between an unconstrained polymer and a looped polymer. The loop exponent  is for an ideal polymer and for an isolated self avoiding loop with in dimensions Gennes1979 . However, helices, which emerge from the loop, increase even further. In the asymptotic limit of long helical sections, renormalization group predicts for a loop with emerging helices Duplantier1986 ; Kafri2000 , where in an expansion. One obtains for terminal, for internal loops and for a loop with four emerging helices. For larger the expansion prediction for becomes unreliable. One sees that the variation of over different loop topologies that appear in the native structures of RNA is quite small, which justifies our usage of the same exponent for loops of all topologies that occur in a given RNA secondary structure.

Figure 1: 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 , here . The thick arrows to either side illustrate the force  applied to the 5’- and 3’-end.
Figure 2: The arc diagram is a representation of the secondary structure depicted in fig. 1. A dot represents one base. Solid lines denote the backbone bonds and thick gray lines the non-nested backbone bonds. Dashed arcs denote hydrogen bonds between two bases. A pseudoknot (dotted arc) is recognized here as crossing arcs. If no pseudoknots are present the structure is hierarchical, meaning that substructures are either nested or juxtaposed.

2.2 Canonical partition function

Figure 3: Illustration of the recursion scheme for the canonical partition function in eq. (3). Boxes denote partition functions of substrands (the range is given by the subindices). The numbers inside a box give the number of non-nested backbone bonds. To calculate the partition function of a strand ranging from through , consider the partition function of a strand ranging from through and add base number , which may (right term in first row) or may not (left term in first row) establish a base pair with base number . In the second row is calculated by closing structures with non-nested bonds with a hydrogen bond (dashed line). For homopolymeric RNA, the sequence dependence drops out and only the lengths of the substrands, , , , need to be considered, see eq. (3).

As we neglect pseudoknots, only hierarchical structures are present, which allows to write down a recursion relation for the partition function. Further, as we consider homopolymers and omit sequence effects by using a constant base pairing weight , the system is translationally invariant. Hence, the canonical partition function of a strand ranging from base at the 5’-end through at the 3’-end depends only on the total number of segments and on the number of non-nested backbone bonds . 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, which couples to an external stretching force and which can be observed for example in force spectroscopy experiments Einert2008 ; Bundschuh2005 ; Mueller2003 , see figs. 1 and 2. The recursion relations for can be written as

(3a)
and
(3b)

and is illustrated in fig. 3. The Heaviside step function is if and if . Eq. (3a) describes the elongation of an RNA structure by either adding an unpaired base (first term) or by adding an arbitrary substrand that is terminated by a helix. Eq. (3b) constructs by closing structures with non-nested bonds, summed up in , by a base pair. denotes the number of configurations of a free chain with links and can be completely eliminated from the recursion relation by introducing the rescaled partition function . We set for computational simplicity and combine eqs. (3a) and (3b), which leads to the final recursion relation

(4)

with the boundary conditions , for , , or . The thermodynamic limit of an infinitely long RNA chain is described by the canonical Gibbs ensemble, which is characterized by a fixed number of segments , but a fluctuating number of non-nested backbone bonds . Therefore we introduce the unrestricted partition function

(5)

which contains the influence of an external force  via the statistical weight  of a non-nested backbone bond. For RNA with no force applied to the ends, one has . We model the RNA backbone elasticity by the freely jointed chain (FJC) model, where the weight of a non-nested backbone bond subject to an external force is given by

(6)

here, we introduced the inverse thermal energy and the Kuhn length .

2.3 Grand canonical partition function

For studying the phase transition and the critical behavior, it is useful to introduce the generating function or grand canonical partition function

(7)

where is the fugacity.

Figure 4: Structure of the grand canonical partition function according to eqs. (7) and (11). The grand canonical partition function of the Gibbs ensemble is a sum over all numbers of non-nested backbone bonds (thick gray lines) with statistical weight . Between two adjacent backbone bonds can be either a single nucleotide (dot), with statistical weight , or a structure with weight , whose first and last base are paired. The white squares serve as wildcards for either possibility and have the statistical weight . Thin black lines depict backbone bonds that are part of a helix or loop and have statistical weight .

Performing the weighted double sum on both sides of eq. (4) yields

(8)

which can be solved and one obtains the generating function

(9)

Here we have defined as the grand canonical partition function of RNA structures with zero non-nested backbone bonds, i. e. structures which consist of just one nucleotide or structures where the terminal bases are paired,

(10)

Eq. (9) has an instructive interpretation, which becomes clear by expanding the fraction in a geometric series

(11)

where is the statistical weight for a backbone segment which is not nested. Between two adjacent segments we have the possibility to put either a single nucleotide (with statistical weight ) or a structure whose first and last bases are paired (with statistical weight ). See fig. 4 for an illustration.

In order to determine the function , we compare the coefficients of the power series in in eqs. (7) and (11) and obtain . The lower summation index is due to exchanging the summations in eq. (7), bearing in mind that for . This identity can be inserted into eq. (10) and yields the equation

(12)

which determines . , for , is the polylogarithm Erdelyi1953 . We introduce

(13)

and rewrite eq. (12) as . Eq. (12) has at most two positive and real solutions as can be seen from fig. 5, where we plot the two sides of eq. (12). is a continuous, monotonically increasing function of with as follows from eq. (10). Therefore, only the smallest positive root of eq. (12) yields the correct . For there is always a positive and real solution for . Increasing increases until eventually, at , the real solution for vanishes. Thus, has a branch point at . Depending on the value of the loop exponent , the polylogarithm and (a) are divergent for , (b) are finite, but feature a diverging slope for , or (c) have a finite value and derivative for at , see fig. 5. This will become important later, when the existence of phase transitions is studied.

Figure 5: Graphical solution of the equations, which determine , , . The functions (dashed line) and (solid lines), eq. (13), are plotted for and different values of the fugacity and the loop exponent , (a)  (b)  (c) . Points at which both curves intersect are solutions of eq. (12) and determine . In case of two positive solutions the smaller yields the correct solution as is to be a continuous, monotonically increasing function of with , eq. (10). Points at which both curves are adjacent to each other (open circles) determine the branch point . Points at which (filled circles) determine the position of the pole  in the absence of force, .

2.4 Back-transform to canonical ensemble

Since the thermodynamic limit is defined in the canonical Gibbs ensemble, we now demonstrate how to obtain , eq. (5), from , eq. (9). For large systems, , the canonical partition function is given by the dominant singularity  of , which is defined as the singularity which is nearest to the origin in the complex -plane Flajolet1990 ; Einert2010 . In particular if with independent of , we obtain

(14)

where is the gamma function Abramowitz2002 . Therefore, the Gibbs free energy reads to leading orders in

(15)

In fact, features two relevant singularities. First, the branch point  of , which is independent of , and second a simple pole , where the denominator of vanishes, see eq. (9). Depending on which singularity has the smallest modulus, the molecule can be in different phases. In the following sections it will turn out that the low temperature, compact or folded phase is associated with , whereas the high temperature, extended or unfolded phase is characterized by .

Let us consider the branch point first. It can be seen from fig. 5, that for at least one real solution of eq. (12) exists, where the smaller solution determines . Right at the two solutions merge and the slope of is at the tangent point. This yields the equation for the position of the branch point singularity , which is a function of only,

(16)

The behavior of in the vicinity of the branch point can be obtained by expanding eq. (12) for and

(17)

where we used the short notation and defined

(18)

Due to the exponent in the above equation, the function exhibits a first order branch point at and the grand canonical partition function, eq. (9), scales as

(19)

Together with eq. (14), we obtain the following scaling for the canonical partition function

(20)

which leads to a logarithmic -contribution with universal prefactor to the free energy , in accord with the findings of de Gennes Gennes1968 . It will turn out that eq. (20) describes the low temperature or folded phase of the system.

Now let us consider the pole singularity of the grand canonical partition function. is a function of and and is given by the zero of the denominator of in eq. (9),

(21)

The position of the pole can be evaluated in a closed form expression by plugging eq. (21) into eq. (12) and solving the resulting quadratic equation for . One obtains

(22a)
(22b)

The behavior of in the vicinity of the pole can be obtained by expanding eq. (12) for , , and

(23)

where we used the short notation and introduced

(24)

Therefore, the grand canonical partition function scales as

(25)

and together with eq. (14) we obtain the scaling of the canonical partition function

(26)

Later we will see, that eq. (26) describes the denatured high temperature phase of the system. In contrast to the branch point phase, eq. (20), no logarithmic contribution to the free energy is present.

3 Critical behavior

The two singularities and are smooth functions of external variables such as temperature  or force , which enter via the weight of a base pair , eq. (1), and the weight of a non-nested backbone bond , eq. (6). As the system is described by the singularity, which is closest to the origin, a phase transition associated with a true singularity in the free energy, eq. (15), is possible if these two singularities cross. For that purpose let us shortly review the three constitutive equations eqs. (12),(16), and (21). As observed earlier, the smallest positive root of eq. (12) yields the function . The simultaneous solution of eqs. (12) and (16) yields the branch point , whereas the simultaneous solution of eqs. (12) and (21) yields the pole , which can be expressed in a closed form, see eq. (22).

Figure 6: Illustration of the graphical solution of eqs. (12),(16), and (21) at the phase transition. The solid line sketches , eq. (13), the dashed line the function . The open circles denote points where the conditions for the branch point are met: Curves have common points, eq. (12), and curves are tangent to each other in these points, eq. (16). The black filled circles denote points where , eq. (21). If the black filled circle lies on the dashed curve the conditions for a pole are met, eqs. (12) and (21). For a given temperature and force the fugacity is increased from to the value where either eqs. (12) and (16) hold (open circle on dashed line, folded phase) or eqs. (12) and (21) hold (black filled circle on dashed line, unfolded phase). (a) Illustration of the thermal phase transition at zero force, . For low temperatures the branch point is dominant. Upon increasing the temperature, i. e. decreasing , the branch point and the pole approach each other until they eventually coincide at the melting temperature and cause a phase transition. For there is no branch point anymore and the pole is dominant. (b) Illustration of the force induced phase transition at a temperature . For small forces the branch point is dominant. Upon increasing the force, i. e. increasing , the point where (black filled circle) moves towards the branch point. As the branch point is independent of the force, see eqs. (12) and (16), no observable depends on as long as . At the branch point and the pole coincide and a phase transition occurs. For the pole is dominant.

3.1 Critical point and existence of a phase transition

The critical fugacity and thus the phase transition is defined as the point where the branch point and the pole coincide

(27)

which means that all three eqs. (12),(16), and (21) have to hold simultaneously, see fig. 6. Assuming that a pair exists so that eq. (27) is true, this can be evaluated further by plugging eqs. (22) into eq. (16) and we obtain

(28)

This constitutes a closed form relation between and or, by employing eqs. (1) and (6), the critical temperature  and force . The melting temperature is defined as the critical temperature at zero force.

The order of the branch point exactly at and zero force, , is calculated by expanding eq. (12) in powers of and while keeping fixed. For vanishing force and we obtain

(29)

where we used

(30)

Thus, the asymptotic behavior of the generating function at and zero force is

(31)

and we obtain the modified scaling of the canonical partition function right at the melting point temperature

(32)

We see that the loop statistics are crucial and enter via the loop exponent , which gives rise to non-universal critical behavior.

For finite force, , the branch point is first order and the scaling of the generating function at the critical point reads

(33)

with

(34)

leading to the canonical partition function

(35)

with a scaling independent of the loop exponent . Note that eq. (33) scales as in contrast to eq. (19), which leads to the different scaling of in eq. (35) when compared to eq. (20). In the rest of this section we compare thermal and force induced phase transition and in particular determine the parameter range in which a phase transition is possible.

3.1.1 Thermal phase transition

First, we consider the thermal phase transition without external force, i. e. for . In this case the polylogarithm reduces to the Riemann zeta function, . Since , we find that as long as has a real, positive branch point. This is due to the fact that (see fig. 5), , and that is a monotonically increasing function of , see eq. (10).

No thermal phase transition for

For the function always features a branch point since for . This ensures that for every a is found, where is tangent to , see figs. 5a and 5b. As the branch point is always dominant, we find the universal scaling

(36)

for all temperatures and no phase transition is possible. The RNA chain is always in the folded phase.

No thermal phase transition for

For the function and its derivative are finite for . A sufficient condition for a branch point to exist is that the slope of is greater than  for , see filled circles in fig. 5c, and hence

(37)

This can be achieved always for large enough  as long as the numerator is positive. On the other hand, for , where is the root of

(38)

no branch point exists since the numerator in eq. (37) is negative. That means that for the pole is always the dominant singularity of and the molecule is always in the unfolded state.

Thermal phase transition for at

Only for a thermal phase transition is possible. For , see eq. (28), the molecule is in the folded phase governed by the branch point singularity , which is determined by eqs. (12) and (16). Decreasing , i. e. increasing the temperature, causes the branch point and the pole to approach each other. At the critical point , eq. (28), both singularities coincide and a phase transition occurs. For higher temperatures, , the RNA is unfolded and described by the pole , eq. (22). See fig. 6a for an illustration. It will turn out that the temperature induced phase transition at zero force is very weak and that, in fact, the order of the phase transition is , where is the integer with .

3.1.2 Force induced phase transition

For the force induced phase transition the situation is slightly different as the position of the pole depends on the force, which enters via the weight of a non-nested backbone bond , eq. (22). In contrast, the branch point does not depend on and hence it is constant, eqs. (12) and (16). Therefore, the branch point  and the critical point  coincide and can be determined exactly by the relation .

No force induced phase transition if or

If the molecule is already in the unfolded phase, which can be due to high temperature, , or due to the non-existence of a branch point, , a force induced phase transition is not possible. In these cases the pole always dominates the system, regardless of the value of the applied force.

Force induced phase transition if and

A system below the melting temperature, , is in the folded phase at zero force, . For small forces, i. e. , the system is described by the branch point singularity , which is independent of and hence does not depend on the force, eqs. (12) and (16) and fig. 6b. However, as the pole is a monotonically decreasing function of , the branch point and the pole will eventually coincide at and a phase transition occurs. The critical force fugacity  is defined as the root of eq. (28) for fixed . For the pole is the dominant singularity and governs the system. Note that in contrast to the thermal phase transition, a force induced phase transition is possible even for . It will turn out that the force induced phase transition is second order in accordance with previous results Mueller2003 ; Montanari2001 .

3.2 Global phase diagrams

Eq. (28) determines the phase diagram. In fig. 7a we show the phase transition between the folded and unfolded states of a homopolymeric RNA in the - plane for a few different values of the loop exponent  that correspond to an ideal polymer, , and values between and as they are argued to be relevant for terminal and internal loops of varying topology including the effects of self-avoidance. Below the transition lines in fig. 7a, the chain is in the unfolded (extended) state, above the line in the folded (compact) state. With growing loop exponent, the extent of the folded phase shrinks. In fact, diverges for , where , cf. eq. (38). Thus, for only the unfolded phase is present, as observed earlier in the restricted case of zero force,  Einert2008 . On the other hand, the critical line for and goes down to zero, for , which indicates that if no external force is applied to the molecule, there is only the folded phase for and hence no thermal phase transition is possible. However, applying a sufficient force can drive the molecule into the unfolded state, as derived earlier. Concluding, only for a thermal phase transition, denoted by the filled circles in fig. 7a, is possible. A force induced phase transition is possible whenever and . For small force and , eq. (28) can be expanded around and yields the universal asymptotic locus of the phase transition

(39)
Figure 7: (a) Phase diagram of homopolymeric RNA in the - plane for different values of the loop exponent featuring an unfolded phase (bottom right) and a folded phase (top left). For , the phase boundary approaches and the melting point diverges; therefore only the unfolded phase exists for . For , there is no melting transition at zero force as . Thus if no force is applied, the system is always in the folded phase regardless of the temperature. The molecule can be denatured, though, by applying an external force even for , as can be seen from the phase boundary for . The filled circles denote the thermal denaturation transition point in the absence of an external force for . (b) Phase diagram in the - plane. With zero force, (solid line), the weight  drops to zero for and diverges as . For finite force, (dashed line), a phase transition is possible, even for . (c) Phase diagram in the - plane for . Below the phase boundary the folded state is present, above the unfolded phase. Re-entrance at constant force is observed, as reported by Müller Mueller2003 .

In fig. 7b we show the critical line for two different values of as a function of the loop exponent . In the absence of an external pulling force, i. e. for (solid line), the transition line only occurs in the limited range . It is seen that for loop exponents around the relevant value of , the critical base pairing weight is quite small and of the order of . A base pairing weight smaller than unity corresponds to a repulsive base pairing free energy that is unfavorable. This at first sight paradoxical result, which means that the folded phase forms even when the extensive part of the base pairing free energy is repulsive, reflects the fact that the folded state contains a lot of topological entropy because of the degeneracy of different secondary structures. The consequences for the theoretical description of systems, where single stranded nucleic acids occur, including DNA transcription, denaturation bubbles in dsDNA, untwisting of nucleic acids Leger1999 , and translocation Bundschuh2005 will be briefly discussed in section 4.

The phase diagram, eq. (28), can also be displayed in the - plane by virtue of eqs. (1) and (6) and is shown in fig. 7c. Here, re-entrance at constant force becomes visible, in line with previous predictions by Müller Mueller2003 . Expanding eq. (6) we obtain , for . Eq. (39) yields the scaling of the critical force close to the melting temperature as

(40)

which depends on the loop exponent  and deviates from the predictions by Müller Mueller2003 , who found a universal exponent .

3.3 Thermodynamic quantities and critical exponents

We now consider the thermodynamic and critical behavior of various quantities. An arbitrary extensive quantity with the conjugate field is obtained from the grand potential  via differentiation with respect to and the chemical potential held constant

(41)

To evaluate the behavior of in the thermodynamic limit, , one sets , where is defined as the chemical potential, at which diverges, i. e.  for . Another route to obtain is to conduct the calculation in the canonical ensemble, i. e. , and to use the dominating singularity Einert2010 , where

(42)

see eq. (15). In fact, for eqs. (41) and (42) are equivalent and is associated with the dominating singularity of , namely , which will be shown now. The Gibbs free energy  and the grand potential  are related via a Legendre transform

(43)

Therefore,

(44)

where the two final expressions are evaluated at . While performing derivatives of the dominant singularity and of the function is straightforward for and , see eq. (22), one has to employ implicit differentiation of eq. (16) to obtain the derivative of and , see supplementary material. For the latter case it turns out that eq. (41) is more convenient to work with.

3.3.1 Fraction of paired bases

Figure 8: Fraction of paired bases as a function of temperature for and various at zero force, . A phase transition is observed only for for the range of positive values of considered here (indicated by the filled circle), since for the critical weight of a hydrogen bond is , which can only be obtained for amounting to a repulsive interaction. For no thermal phase transition can be observed at all. The inset shows the third derivative for , which reveals the phase transition.
Figure 9: Fraction of paired bases as a function of force for various and . The phase transition is visible as a kink in the curves and indicated by a filled circle. (a)  and varying . (b)  and varying .

The fraction of paired bases is

(45)

We obtain

(46)

in the folded phase (, ) and

(47)

in the unfolded phase (, ). In fig. 8 the temperature dependence and in fig. 9 the force dependence of is shown. The singularity at the critical point of the thermal phase transition for zero force, , is very weak and becomes visible in the th derivative, with being the integer with , see supplementary material. The th derivative exhibits a cusp, see inset of fig. 8,

(48)

which is characterized by the critical exponent for and for , see supplementary material. The force induced phase transition is continuous, too, yet it exhibits a kink in

(49)

which is characterized by the exponents for and for . Therefore, the force induced phase transition is second order in accordance with previous results Mueller2003 ; Montanari2001 . We note that for , i. e. , a finite fraction of bases are still paired. The situation is different for the force induced phase transition where for .

3.3.2 Specific heat

Figure 10: Specific heat as a function of temperature for different loop exponents . The non-critical behavior of the heat capacity curve depends on the loop exponent only marginally. However, the existence and position of the critical point (indicated by the filled circle) and the critical behavior depend on . The inset depicts the third derivative of the specific heat, , for revealing the phase transition.

The specific heat is defined as

(50)

One observes that the specific heat in fig. 10 exhibits only a very weak dependence on the loop exponent, which stands in marked contrast to the findings for the short explicit sequence of tRNA-phe Einert2008 ; Einert2010a , where a pronounced dependence of the heat capacity on is observed. The non-analyticity of , eq. (48), translates into a divergence of the th derivative of the specific heat at the melting temperature

(51)

with the critical exponent for and for . This singularity is illustrated in the inset of fig. 10 for .

3.3.3 Fraction of non-nested backbone bonds

Figure 11: (a) Fraction of non-nested backbone bonds  as a function of temperature for . For the critical weight of a hydrogen bond is , which can only be obtained for , amounting to a repulsive interaction. Thus, for or all segments are parts of loops or helices and hence . Filled circles indicate the melting temperature . (b) and (c) fraction of non-nested backbone bonds as a function of force for various and . Again, the phase transition is visible as a kink in the curves and is indicated by a filled circle. (b)  and varying . (c)  and varying . Filled circles indicate the position of the phase transition.

The fraction of non-nested backbone bonds is obtained by

(52)

and is

(53)

in the folded phase, as does not depend on , and reads

(54)

in the unfolded phase. For the fraction of non-nested backbone bonds assumes a finite value smaller than one, which again indicates that the denatured phase in our model features pronounced base pairing. However, for large force on invariably obtains . As can be seen nicely in fig. 11, both and feature a kink at the critical point.

3.3.4 Force extension curve

Figure 12: Force extension curve as a function of force for various and . The phase transition is indicated by filled circles and occurs at zero extension and a finite threshold force needed to unfold the compact folded structure. (a)  and varying . (b)  and varying . Additionally the force extension curve of a freely jointed chain is plotted, which is the limiting form for .

The force extension curve is closely related to the fraction of non-nested backbone bonds . The extension per monomer is given by

(55)

Since the Langevin function is a smooth function, the critical behavior of is governed by the behavior of the fraction of non-nested bonds . As can be seen in fig. 12a, the stretching behavior of the Langevin function is approached as the base pairing weight decreases, otherwise pronounced deviations are seen in the force-stretching curves. Also, a finite stretching force to unravel the folded state is needed. From fig. 12 it becomes obvious that the force induced phase transition is second order as the force extension curve exhibits a kink at the critical force denoted by filled circles. The force extension curves are in accordance with previous results Montanari2001 ; Mueller2003 .

4 Implications for DNA melting

How do the previous results impact on the theoretical description of the denaturation of double stranded nucleic acid systems, particularly DNA melting? When double stranded DNA approaches the denaturation transition, more and more inter-strand base pairs break up and loops proliferate. In the traditional theories based on the Poland-Scheraga model Poland1966 ; Poland1966a , the possibility of intra-strand base pairing was not considered. These models are thus accurate for duplexes formed between strands with sequences and , where indeed base pairs (between A and T and between G and C) can only form between the two strands, not within one strand. For the case of duplexes formed by two strands with the sequence or , both intra- and inter-strand base pairs can form and have identical statistical weights. In this case, our model predicts significant modifications for the duplex melting scenario.

The above sequence examples are prototypes for two extreme cases of the general scenario characterized by statistical weights and for intra-strand and inter-strand pairing, respectively. Duplexes formed between and are characterized by , whereas duplexes formed between two or strands are characterized by . Intermediate values of and can be achieved experimentally, for example, by and its complementary sequence , in which case only of the intra-strand base pairs can be of the Watson-Crick type and which would effectively lead to a lower weight of intra-strand base pairs, i. e. . In naturally occurring DNA, a similar situation might be present above the glass transition where, for certain sequences, self-hybridization in a single strand is possible to some extent Bundschuh2002a . The weights of inter-strand and intra-strand base pairs can also be changed by applying an external force or torque on the duplex, for instance in the setup by Léger et al. Leger1999 . We expect the weight of inter-strand base pairs to decrease when the duplex is untwisted. This might lead to denatured regions in the duplex, where secondary structure can form if the sequence allows. We speculate that subsequent pulling first leads to the denaturation of the secondary structures, whose signature would be a threshold force around - Montanari2001 ; Einert2010a , see fig. 12, followed by the over-stretching transition of DNA Einert2010 . According to our previous arguments, a marked dependence on the sequence should be observed.

Figure 13: Illustration of a double stranded DNA molecule near the melting transition. (a) If intra-strand base pairing is possible, in inter-strand loops the vast majority of bases will form intra-strand secondary structure elements. (b) If the sequence does not favor intra-strand base pairs inter-strand loops without secondary structure will form. (c) Illustration of an asymmetric loop.

Let us now review the classical Poland-Scheraga calculation Poland1966 ; Poland1966a and allow – in addition to the standard treatment – for intra-strand base pairing in denatured regions of the duplex, see fig. 13a. The grand canonical partition function of the two-state Poland-Scheraga model reads

(56)

where

(57)

is the grand canonical partition function of a duplex in the bound state with the weight of an inter-strand base pair and now plays the role of the fugacity of a base pair. We distinguish between three different scenarios for a denatured or molten region inside the double strand, which is characterized by the partition function , depending on the sequence and the intra-strand base pairing weight : (i) for , with defined by eq. (28), intra-strand secondary structures in the folded (low temperature) phase are present, (ii) for intra-strand secondary structures in the unfolded (high temperature) phase form, and (iii) for large denaturation bubbles occur without intra-strand base pairing, corresponding to the original Poland-Scheraga model. The grand canonical partition for case (i) is characterized by the branch point singularity, eq. (20),

(58)

with