Thermal denaturation of fluctuating finite DNA chains: the role of bending rigidity in bubble nucleation
aa r X i v : . [ c ond - m a t . s o f t ] S e p Thermal denaturation of fluctuating finite DNA chains:the role of bending rigidity in bubble nucleation
John Palmeri, Manoel Manghi and Nicolas Destainville
Laboratoire de Physique Th´eorique, Universit´e de Toulouse, CNRS, 31062 Toulouse, France (Dated: 18 September 2007)Statistical DNA models available in the literature are often effective models where the base-pairstate only (unbroken or broken) is considered. Because of a decrease by a factor of 30 of the effectivebending rigidity of a sequence of broken bonds, or bubble, compared to the double stranded state, theinclusion of the molecular conformational degrees of freedom in a more general mesoscopic model isneeded. In this paper we do so by presenting a 1D Ising model, which describes the internal base pairstates, coupled to a discrete worm like chain model describing the chain configurations [J. Palmeri,M. Manghi, and N. Destainville, Phys. Rev. Lett. , 088103 (2007)]. This coupled model isexactly solved using a transfer matrix technique that presents an analogy with the path integraltreatment of a quantum two-state diatomic molecule. When the chain fluctuations are integratedout, the denaturation transition temperature and width emerge naturally as an explicit functionof the model parameters of a well defined Hamiltonian, revealing that the transition is driven bythe difference in bending (entropy dominated) free energy between bubble and double-strandedsegments. The calculated melting curve (fraction of open base pairs) is in good agreement with theexperimental melting profile of polydA-polydT and, by inserting the experimentally known bendingrigidities, leads to physically reasonable values for the bare Ising model parameters. Among thethermodynamical quantities explicitly calculated within this model are the internal, structural, andmechanical features of the DNA molecule, such as bubble correlation length and two distinct chainpersistence lengths. The predicted variation of the mean-square-radius as a function of temperatureleads to a coherent novel explanation for the experimentally observed thermal viscosity transition.Finally, the influence of the DNA strand length is studied in detail, underlining the importanceof finite size effects, even for DNA made of several thousand base pairs. Simple limiting formulæ,useful for analyzing experiments, are given for the fraction of broken base pairs, Ising and chaincorrelation functions, effective persistence lengths, and chain mean-square-radius, all as a functionof temperature and DNA length. PACS numbers: 87.10.+e General theory and mathematical aspects, 87.15.Ya Fluctuations, 82.39.Pj Nucleicacids, DNA and RNA bases
I. INTRODUCTION
The stability of double-stranded DNA (dsDNA) at physiological temperature is due to the self-assembly of its basepairs: self-assembly within a same strand via base-stacking interactions between neighboring bases; and self-assemblyof both strands via hydrogen bonds between pairs of complementary bases. These interactions, however. are on theorder of magnitude of a few k B T (thermal energy) [1, 2, 3] and thermal fluctuations can lead, even at physiologicaltemperature, to local and transitory unzipping of the double strand (see e.g. [4] or the reviews [5, 6]). The cooperativeopening of a sequence of consecutive base pairs leads to denaturation bubbles which are likely to play a role from abiological perspective, since they may participate in mechanisms such as replication, transcription or protein binding.For example, it has been proposed [7] that transcription start and regulatory sites could be related to DNA regionswhich have a higher probability of promoting bubbles. Indeed, the energy needed to break an adenine-thymine (A-T)base pair ∼ k B T , connected by two hydrogen bonds, is smaller than the energy needed to break a guanine-cytosine(G-C) one ∼ k B T (3 hydrogen bonds) [1, 2]. At the same temperature, A-T rich sequences present a priori morebubbles than G-C rich ones, even though sequence effects on the occurrence of bubbles are more complex than a simpleexamination of local A-T base abundance [1, 7]. In addition to the sequence, the fraction of denaturation bubbles in vitro naturally depends on temperature, as well as on the ionic strength of the solution [5, 8]. In particular, themelting temperature T m , above which bubbles proliferate and the two strands completely separate, depends on bothsequence and ionic strength. Another parameter that affects the melting temperature is the length of the doublestrand [9, 10]. This is not a purely academic debate because short DNA strands (a few tens of base pairs) are involvedin DNA chip experiments where the hybridization process is precisely affected by temperature in a way depending onstrand length and sequence (see [11] and references therein).Although the intracellular unwinding of DNA is due to active and enzymatic processes by imposing unwindingtorsional stresses [12], the thermally induced denaturation of purified DNA in solution has led to an intensive studyof DNA thermal denaturation [4, 5, 6, 13, 14]. Mesoscopic models have been proposed to account for the thermody-namical properties of denaturation bubbles in DNA. The first models were Ising-like two-state models, where the basepairs can be open or closed (see [4, 5] and references therein). In the simple base-pair model, the Ising parametersare the base-pair chemical potential and the so-called cooperativity parameter, which accounts for the energetic costof a domain wall. This type of one-dimensional Ising model is exactly soluble for homopolymers and for randomsequences [4]. More sophisticated effective Ising models, including the 10 kinds of base-pair doublets, lead to betteragreement with experimental data [5]. Poland and Scheraga [14], following the work of Zimm [15], included “loop en-tropy” in these models, i.e. the entropic cost of closing large loops formed by destacked single-stranded DNA (ssDNA)in bubbles [4, 14]. This Poland-Scheraga model is at the core of the DNA melting simulator MELTSIM algorithmfirst developed by Blake et al. [16, 17]. The intra-loop and inter-loop self-avoidance corrects the loop statistics andrefines the sharpness of the denaturation transition [18, 19]. To get a denaturation transition in these models, an ef-fective temperature dependent Ising chemical potential must be inserted by hand. Recently, Peyrard et al. developednon-linear phonon models where the shape of the interaction potential between base pairs is more precisely taken intoaccount (see [20, 21, 22] and the review [23]). By inserting estimated microscopic parameters, they have shown thattheir model can both lead to a denaturation transition, analogous to interface unbinding, and be useful for studyingbubble dynamics. The predicted transition temperature is, however, very sensitive to the model parameter values,and if physically reasonable values are used [21, 24, 25], T m appears to be much too high and the transition widthmuch too large.More generally, the theoretical study of DNA denaturation can in principle be tackled on a least four levels ofinvestigation determined by the amount of detail included, going from (1) quantum ab initio approaches and (2)classical all atom molecular dynamical simulations [26], through (3) effective mesoscopic approaches coupling chainconformational and base-pair degrees of freedom [27], to (4) effective statistical models for base-pairs alone ([4, 14, 23]and references therein). In theory, it is possible to move up one step in this hierarchy by integrating out the subsetof the degrees of freedom that do not appear at the higher level, giving rise to en effective free energy at each levelof description. Our purpose here is to show, via a minimal model for DNA homopolymers, that if one starts at thethird level and integrates out the chain degrees of freedom, one arrives at a physically coherent level (4) explanationfor the DNA melting transition: a bending free energy driven denaturation transition emerges naturally due to theentropic lowering of the energetic barrier for bubble nucleation. Working at level (3) also has the added advantageof allowing access to the statistics of the chain degrees of freedom (effective persistence length, mean-square-radius,etc.), something that obviously is not possible at level (4).To this end, we have recently proposed a model [28], which considers not only the internal coordinates in termsof Ising spin variables describing the open or closed states, but also external coordinates, the chain tangent vectors,which determine the chain configuration and depend sensitively on chain stiffness. Indeed, ssDNA is two orders ofmagnitude more flexible than dsDNA at normal salt concentration. We have shown that this difference in bendingrigidity provides a novel explanation for the bubble mechanism formation. Further evidence for the importanceof this bending heterogeneity include phenomena such as cyclization, loop formation, and packaging of DNA intonucleosomes [29] (where denaturation bubbles facilitate bending of the otherwise rigid polymer DNA in structureswhere it coils up with curvature radii down to 10 nm, despite a persistence length of dsDNA equal to 50 nm). Ourmodel incorporates precisely this dependence of the polymer bending rigidity on the state of neighboring base pairs.It is the discretized version of a continuous model [28], and couples explicitly an Ising model, describing the internaldegrees of freedom (open or closed), and a Heisenberg or discrete worm like chain (DWLC) one, accounting forthe rotational degrees of freedom between successive monomers of the DNA chain. Its originality lies in the factthat the internal 2-state and external bending degrees of freedom are treated on an equal footing and therefore therenormalized Ising parameters obtained by integrating out the chain can be exactly calculated within the scope of themodel. The melting temperature T m naturally emerges and, together with the transition width ∆ T m , are explicitlywritten as a function of the bending rigidities and strand length, which are experimentally known, and bare Isingparameters. In addition, the effective Ising properties (fraction of broken bases and correlation length) as well asthe chain ones (persistence length and mean square radius) can be computed, allowing in principle direct comparisonwith experiments. An important feature of this effective Ising model is that the end monomers see an effectivechemical potential that is lower than the interior one; this end-interior asymmetry leads in a natural way to a chainlength dependence for T m , the transition width, and other statistical quantities. A similar model had been previouslyproposed by one of us in 2D [30], but its application to the study of denaturation bubbles in dsDNA was not madeexplicit. Recently, a similar approach has been considered in the context of the dsDNA stretching transition [31].The addition in the energy functional of the term corresponding to the external force prevents an exact solution ofthe model in Ref. [31] and an approximate variational scheme had to be implemented.Beyond dsDNA or dsRNA, our coupled model can be used to describe the properties of any two-state biopolymer,as soon as the local bending rigidity depends on the local states. As already mentioned, the transition from B-to S-form of dsDNA in force experiments has been investigated in this framework [31]. The helix-coil transition inpoly-peptides can also be described by such a theory because the α -helix configuration is much more rigid than therandom one [9, 14].The present paper is a detailed account of the results summarized in a Letter [28]. In section II, we present thecoupled classical Ising-Heisenberg model, which, as we show here, can be used to describe DNA thermal denaturationand write the partition function in terms of a transfer matrix, usual in one-dimensional statistical systems. We notein passing that the coupled Ising-Heisenberg model presented here displays a rich array of behavior and therefore maybe of interest in other contexts, such 1D classical spin chains or 0D quantum rotators describing a diatomic moleculewith internal states.Because the full transfer matrix method for the coupled model leads to relatively complex calculations, we firstshow that the model can be reduced to two effective Ising models, a path that provides a great deal of physicalinsight: indeed, these effective models allow the calculation of the free energy, as well as the Ising and chain end-to-end tangent-tangent correlation functions in terms of an effective Hamiltonian with temperature-dependent Isingparameters. A detailed solution of the two effective models is provided in section III. Section IV is devoted to thecalculation of the Ising and chain correlation quantities, using the effective Ising models, as well as a novel correlatormixing both Ising and chain variables. The next sections, V and VI, present in detail the full transfer matrix approach,which leads to the complete calculation of chain correlations and end-to-end distance; a visual interpretation of theexpression for the 2-point correlation functions leads naturally to an analogy with a quantum diatomic molecule. Ourtheory is compared to experimental denaturation profiles of synthetic DNA in section VII and finite-size effects, whichare experimentally relevant, are thoroughly examined in section VIII. Finally, our concluding remarks are given insection IX, where we also summarize our principal theoretical results of greatest interest for interpreting experiments.The principal symbols used in this work are defined and catalogued in Table I at the end of the paper. II. DISCRETE CHAIN MODEL
We model dsDNA as a discrete chain of N monomers (links), each monomer can be in one of two different states,U and B, which denote, respectively, unbroken and broken bonds. The local chain rigidity depends on the nearbylink types. A denaturation bubble is thus formed by a consecutive sequence of B type monomers. The chain’sconformational properties are determined by the set of N unit link tangent vectors { t i : i = 1 , . . . , N } with k t i k = 1.For simplicity, the monomer length, a , is taken to be the same for both U and B (for modeling DNA stretchingtransitions it is necessary to introduce different monomer lengths [9, 31]). The position of the end of the i th link inthe 3D embedding space is X i = X + a P ij =1 t j , where X is an arbitrary starting point. The end-to-end vectoris R = a P Nj =1 t j . The link states are denoted by the value of an Ising variable σ i = ± N variables { σ i , t i } . For a chain in 3D eachlink vector can be expressed in spherical coordinates as t i = (sin( θ i ) cos( φ i ) , sin( θ i ) sin( φ i ) , cos( θ i )) and can thereforebe defined by the azimuthal and polar angles φ i and θ i , denoted together by the solid angle Ω i = ( θ i , φ i ). The energy H [ σ i , t i ] of a state is taken to be H [ σ i , t i ] = N − X i =1 ˜ κ i +1 ,i (1 − t i +1 · t i ) − N − X i =1 " ˜ Jσ i +1 σ i + ˜ K σ i +1 + σ i ) − ˜ µ N X i =1 σ i . (1)The angle γ i,j between two tangent vectors is given bycos γ i,j = t i · t j = sin( θ i ) sin( θ j ) cos( φ i − φ j ) + cos( θ i ) cos( θ j ) . (2)The first term in H is the bending energy of a DWLC with a local rigidity ˜ κ i +1 ,i , having the dimension of energy, thatdepends on the neighboring values of the Ising variables. We have taken, without any loss of generality, the minimum FIG. 1: Illustration of the different Ising parameters appearing in the Hamiltonian βH . Open (closed) base pairs are coded bya spin σ = +1 ( − σ are setto +1. The first line shows the cost, 2 J , of a domain wall. The second line indicates the energy, 2 µ , required to open a basepair. The third line gives the difference in stacking energy between a segment of dsDNA and a denaturated one: dark (light)blue cigars indicate stacked state in dsDNA (in ssDNA) and the absence of dots indicates the destacking of adjacent base pairs,which is already taken into account by the J contribution. of the bending energy to be zero independent of the values of the Ising variables. The second and third terms makeup the energy of the Ising model [4, 5], H I ≡ H I ( ˜ J , ˜ K, ˜ µ ), illustrated in Fig. 1.The term in ˜ J accounts for the local destacking energy (2 ˜ J ) of a domain wall (where σ i passes from one valueto another). The term in ˜ K accounts for the difference in stacking energy between a segment of dsDNA and of adenaturation bubble. The last term gives the energy (2˜ µ ) required to create a link in the state σ i = − βH [ σ i , Ω i ], where β = 1 / ( k B T ), which thus contains thedimensionless parameters κ i +1 ,i ≡ β ˜ κ i +1 ,i , J ≡ β ˜ J , K ≡ β ˜ K , and µ ≡ β ˜ µ . The local rigidity is κ i +1 ,i = κ U for U − U n . n .κ B for B − B n . n .κ UB for U − B or B − U n . n . (3)where “U-U n.n.”, etc. denotes nearest neighbor link types. In terms of the Ising field σ i κ i +1 ,i = 14 ( κ U + κ B − κ UB ) σ i +1 σ i + 14 ( κ U − κ B )( σ i +1 + σ i ) + 14 ( κ U + κ B + 2 κ UB ) . (4)We identify the B state with two identical non-interacting single DNA strands, ssDNA, each with a local rigidityequal to κ B /
2. From Eqs. (1) and (4), the result of the coupling between Ising and tangent variables can alreadybe predicted: both the destacking and stacking parameters ˜ J and ˜ K will be modified by the two first terms inEq. (4) which have exactly the same functional form in σ i , whereas ˜ µ will remain unchanged. Moreover, when allthe bending rigidities are equal, κ U = κ B = κ UB , the two first terms in Eq. (4) disappear and the Hamiltonian (1)decouples into a pure Ising Hamiltonian and a pure DWLC one (isomorphic to a 1D classical Heisenberg model formagnetism [32]): H [ σ i , t i ] = H DWLC [ t i ] + H I [ σ i ]. In the language of magnetism the model studied here is a classicalcoupled Heisenberg-Ising spin chain. Although pure effective Ising models have been used extensively to model helix-coil and denaturation (melting) transitions in macromolecules, it was necessary to introduce phenomenologically aneffective temperature-dependent chemical potential to obtain a melting transition. A key feature of the coupled modelused here is that a melting transition will emerge naturally in the effective Ising model obtained by integrating outthe chain conformational degrees of freedom.The quantities that we will use to study the behavior of the coupled system are the mean of the internal statevariable (“magnetization” in spin language) c ≡ N N X i =1 σ i , (5)the local state average, h σ i i , correlation functions for the Ising variables, h σ i + r σ i i , and the chain tangent vectors, h t i + r · t i i , and the mean square radius R ≡ h R i / . These correlation functions measure the extent of cooperativityexhibited by the coupled system: e.g., the size of the B (U) domains below (above) the melting temperature, and thelength scale on which the chain remains rigid. The concentration of U and B links is given by ϕ B ( N, T ) = 1 − h c i ( N, T )2 = 1 − ϕ U ( N, T ) . (6)Once the type of homopolymeric DNA is chosen, the bare Ising parameters and the chain bending rigidities can beconsidered fixed, and therefore h c i , ϕ U , and ϕ B become functions of the experimental control parameters, namelytemperature, T , and chain length, N . For a pure U chain ϕ U = 1 and ϕ B = 0; for a pure B chain ϕ U = 0 and ϕ B = 1.A chain with a finite concentration of bubbles (B links) will have ϕ B > T m < ∞ , if itexists, will be defined by ϕ U ( T m ) = ϕ B ( T m ) = 1 / O = O [ σ i , Ω i ]that depends on the fluctuating degrees of freedom, [ σ i , Ω i ], is given by hOi ≡ (4 π ) N Z X { σ i = ± } N Y i =1 Z d Ω i π O [ σ i , Ω i ] e − βH [ σ i , Ω i ] , (7)where Z = (4 π ) N X { σ i } N Y i =1 Z d Ω i π e − βH [ σ i , Ω i ] (8)is the partition function. The partition and correlation functions for the coupled system can be calculated usingtransfer matrix techniques. For example, the partition function can be written as Z = (4 π ) N X { σ i } N Y i =1 Z d Ω i π h V | σ ih σ | ˆ P (Ω , Ω ) | σ i · · · h σ N − | ˆ P (Ω N − , Ω N ) | σ N ih σ N | V i , (9)where the transfer kernel that appears N − P (Ω i , Ω i +1 ) = (cid:18) e κ U [cos( γ i +1 ,i ) − J + K + µ e κ UB [cos( γ i +1 ,i ) − − J e κ UB [cos( γ i +1 ,i ) − − J e κ B [cos( γ i +1 ,i ) − J − K − µ (cid:19) . (10)It is written in the canonical base | U i = | + 1 i and | B i = | − i of the U and B states. The end vector | V i = e µ/ | U i + e − µ/ | B i (11)enters in order to take care of the free chain boundary conditions. Different boundary conditions could be easilyhandled in a similar way, for instance for closed (open) ends, | V i = | U i ( | B i ) and all the following results notexplicitly using Eq. (11) remain valid.Before presenting the full transfer matrix method, we first show that the partition function and averages of anyquantities depending only on the Ising variables can be obtained by examining the effective Ising model obtainedby integrating over the chain conformational degrees of freedom (the link tangent vectors) in Eq. (9). The problemreduces to that of an effective Ising model with an “effective free energy” H (0)I , eff containing renormalized parameters.This method works because, for the coupled Ising-chain system, the rotational symmetry is not broken (absence ofa force term ∝ t i · ˆz in the Hamiltonian [9, 31]). Hence the matrix obtained by integrating the kernel ˆ P (Ω i , Ω i +1 )in Eq. (9) is the same for any site i . We thus are able to carry out the solid angle integrations in sequential fashionby using the ( i + 1) th tangent vector as the polar axis for the i th solid angle integration. The solid angle integratedtransfer matrix is ˆ P (0)I , eff = Z d Ω i π ˆ P (Ω i , Ω i +1 ) = (cid:18) e − G ( κ U )+ J + K + µ e − G ( κ UB ) − J e − G ( κ UB ) − J e − G ( κ B )+ J − K − µ (cid:19) (12)where G ( κ ) is the (dimensionless) Helmholtz free energy of a single joint (two-link) subsystem with rigidity κ (eitherU-U, B-B, U-B or B-U): G ( κ ) = − ln (cid:26)Z d Ω4 π exp[ κ (cos( θ ) − (cid:27) = κ − ln (cid:20) sinh( κ ) κ (cid:21) , (13)an increasing function of κ , varying linearly with κ for κ ≪ T ) and as ln(2 κ ) for κ ≫ T entropydominated regime where the spin wave approximation for the chain degrees of freedom is valid). The effective transfermatrix ˆ P (0)I , eff can be written in Ising form using the renormalized Ising parameters, J , K and the prefactor exp( − Γ ),all depending on chain rigidities: ˆ P (0)I , eff = e − Γ (cid:18) e µ + K + J e − J e − J e − µ − K + J (cid:19) (14) J ≡ J −
14 [ G ( κ U ) + G ( κ B ) − G ( κ UB )] (15) K ≡ K −
12 [ G ( κ U ) − G ( κ B )] (16)Γ ≡
14 [ G ( κ U ) + G ( κ B ) + 2 G ( κ UB )] . (17)In the limit of high temperature, the chain tangents will be completely decorrelated and the uninteresting renormal-izations of J and K arise solely from κ i +1 ,i (Eq. 4), in agreement with Eqs. (15) and (16), as revealed by the lineardependence of G ( κ ) on κ in this limit. It is rather in the opposite limit of low temperatures and therefore stronglycorrelated chain tangents that the bending entropy driven denaturation transition arises.The full partition function, Z = Z (0)I , eff , in Ising transfer matrix notation, Z (0)I , eff = (4 π ) N X { σ i } h V | σ ih σ | ˆ P (0)I , eff | σ i · · · h σ N − | ˆ P (0)I , eff | σ N ih σ N | V i , (18)can be rewritten explicitly in terms of an effective Ising free energy, H (0)I , eff : Z (0)I , eff = (4 π ) N e − ( N − X { σ i } e − βH (0)I , eff [ σ i ] (19)where H (0)I , eff = H I ( ˜ J , ˜ K , ˜ µ ): βH (0)I , eff = − J N − X i =1 σ i +1 σ i − K N − X i =1 ( σ i +1 + σ i ) − µ N X i =1 σ i = − J N − X i =1 σ i +1 σ i − L N − X i =1 ( σ i +1 + σ i ) − µ σ + σ N ) (20)with L ≡ µ + K = µ + K − ∆ G UB / , (21)and ∆ G UB ≡ G ( κ U ) − G ( κ B ). Because H (0)I , eff depends on the temperature, it cannot be considered as an Isingstate energy, but rather as an effective free energy obtained by integrating out the chain subsystem (cf. discussionconcerning levels of theoretical study in the Introduction).When two links open the renormalized stacking energy of the links is K , which is smaller than K by the differencein bending free energy between U-B and B-B joints. The second form of βH (0)I , eff in Eq. (20) shows that the effectivechemical potential of one interior base pair is renormalized to L . If the gain in the one link bending free energy ingoing from U to B, ∆ G UB , is greater than the intrinsic energy, 2(˜ µ + ˜ K ), needed to break a closed interior bond,then the effective interior joint chemical potential L can become negative, signaling a change in stability of U andB states. The end links ( i = 1 , N ), however, feel a different chemical potential, µ + K /
2, which is larger than L inthe case of interest ( κ U > κ B ). This end-interior asymmetry, along with the extra bubble initiation energy due to thesecond domain wall, are reflected in the difference between the “effective free energy” needed to create an n -bubbleat a chain end, β ∆ G ( n )end = 2 J − K + 2 nL (22)and in the chain interior β ∆ G ( n )int = 4 J + 2 nL , (23)which is higher than for an end link by 2 J + K . As will be confirmed in section VIII, this difference in effective freeenergy suggests that at sufficiently low T bond melting will begin at the chain ends. Written out in greater detail,Eq. (23) leads to∆ G ( n )int = 4 ˜ J + 2 nk B T L = 4 ˜ J + 2 n (˜ µ + ˜ K ) − k B T [ G ( κ U ) + G ( κ B ) − G ( κ UB )] − nk B T ∆ G UB . (24)In the uncoupled limit ( κ U = κ B = κ UB ) only the first two (temperature independent) terms survive. Althoughthis end-interior asymmetry plays no role in the limit of an infinite chain ( N → ∞ ), it will play an essential rolein determining how bond melting varies with temperature and bond location (melting maps) and how the meltingtemperature varies with chain size. These finite size effects are discussed in detail in section VIII.Because the renormalized Ising parameters depend on the chain parameters, the coupled model does not in generalhave the same behavior as the uncoupled one. Indeed, since G ( κ ) is an increasing function of κ and κ U ≫ κ B fordsDNA, if the difference between κ U and κ B is sufficiently large, then at a certain temperature, T ∞ m , the effectiveinterior bubble chemical potential, L , can vanish. Provided that this temperature be sufficiently low for thermaldisorder to be weak and end effects due to the finite size of the chain to be unimportant, the state of the systemwill flip from nearly pure U for T < T ∞ m , where L >
0, to nearly pure B for
T > T ∞ m , where L < T = T ∞ m , there will be, on average, as many closed as open bonds and ϕ B, ∞ = ϕ U, ∞ = 1 /
2, where ϕ B, ∞ ( T ) ≡ lim N →∞ ϕ B ( N, T ), etc. This transition (strictly speaking a crossover), which can be extremely sharp under somecircumstances (see below), can be interpreted as a melting or denaturation transition. From the above analysis, wesee clearly that it is the difference in free energy between U-U and B-B joints, ∆ G UB , that drives the melting transition.We will see below, moreover, that if κ U and κ B are much greater than one, then the spin wave approximation is valid,and the entropy term in ∆ G UB dominates. We will see in section VII that this is actually the case for real DNA.By calculating the average of the product of the cosines of the N − Q i cos( γ i +1 ,i ) and using the sametechnique used above for Z (0)I , eff , we can define another effective model, but now with partition function Z (1)I , eff = Z * N − Y i =1 t i +1 · t i + . (25)This expression can be written in Ising transfer matrix form as Z (1)I , eff = (4 π ) N X { σ i } h V | σ ih σ | ˆ P (1)I , eff | σ i · · · h σ N − | ˆ P (1)I , eff | σ N ih σ N | V i , (26)where ˆ P (1)I , eff = Z d Ω i π cos( γ i +1 ,i ) ˆ P (Ω i , Ω i +1 ) = e − Γ (cid:18) e µ + K + J e − J e − J e − µ − K + J (cid:19) (27) J ≡ J −
14 [ G ( κ U ) + G ( κ B ) − G ( κ UB )] (28) K ≡ µ −
12 [ G ( κ U ) − G ( κ B )] (29)Γ ≡ . (30)The function G ( κ ) = − ln { R d Ω4 π cos( θ ) exp[ κ (cos( θ ) − } is related to the tangent-tangent correlation functionbetween two adjacent monomers (isolated 2-link sub-system) with rigidity κ : h t · t i − link = h cos( θ ) i − link = exp [ − G ( κ ) + G ( κ )] = u ( κ ) = coth( κ ) − /κ, (31)which is the Langevin function [33]. It increases with κ , varying as κ/ κ ≪ − /κ for κ ≫
1. Thisasymptotic behavior corresponds to G ( κ ) − G ( κ ) varying as ln(3 /κ ) for κ ≪ κ for large κ . For a purechain (U or B) h cos( θ ) i − link is equal to the nearest-neighbor tangent correlation function h t i +1 · t i i = exp( − /ξ p )with persistence length ξ p ( κ ) = − / ln[ u ( k )] = [ G ( κ ) − G ( κ )] − . (32)In the high κ (spin wave) approximation, we obtain ξ p ≃ κ , as expected.The partition function can be rewritten explicitly in terms of an effective Ising free energy, H (1)I , eff : Z (1)I , eff = (4 π ) N e − ( N − X { σ i } e − βH (1)I , eff [ σ i ] , (33)where βH (1)I , eff = − N − X i =1 (cid:20) J σ i +1 σ i + L σ i +1 + σ i ) (cid:21) − µ σ + σ N ) . (34)By repeatedly using the vector identity ( a · c )( b · d ) = ( a · d )( b · c ) − ( a × b ) · ( c × d ) for a · b = t i +1 · t i , etc., alongwith the property that averages of cross products, h t i +1 × t i i , are zero (and that | t i | = 1), the partition function Z (1)I , eff can be written as the product of the end-end tangent-tangent correlation function, h t · t N i , and the effectiveIsing partition function, Z (0)I , eff : Z (1)I , eff = h t · t N i Z (0)I , eff . (35)When κ U = κ B = κ UB = κ , we recover the pure chain tangent-tangent correlation function: h t · t N i = Z (1)I , eff Z (0)I , eff → (cid:26) exp[ − G ( κ )]exp[ − G ( κ )] (cid:27) N − = exp [ − ( N − /ξ p ( κ )] . (36)Coming back to the difference in (dimensionless) free energy ∆ G UB , we can show that at room temperature it isdominated by its entropic part. Indeed, G ( κ ) can be split into an average (dimensionless) energy and average entropycontribution, G ( κ ) = E ( κ ) − S ( κ ) /k B , where E ( κ ) = β ∂G /∂β = κ ∂G ( κ ) /∂κ . Hence, the average energy of atwo-link system can be written in terms of G and G : E ( κ ) = κ [1 − exp ( − G ( κ ) + G ( κ ))] = κ [1 − exp ( − /ξ p ( κ ))] = κ [1 − u ( κ )] , (37)and therefore ∆ E UB ≡ E ( κ U ) − E ( κ B ) = ( κ U − κ B ) − [ κ U u ( κ U ) − κ B u ( κ B )]. Because the function u ( κ ) tendsto 1, ∆ E UB →
0, and therefore for temperatures low enough for the spin wave approximation to be valid for bothU-U and B-B links, we see that ∆ G UB ≃ − ∆ S UB /k B . Indeed in this approximation, the Hamiltonian is Gaussianand equipartition of energy occurs: h ˜ E i ∼ β − ⇒ h E i = β h ˜ E i ∼
1. In this case the melting transition, if it exists,will be driven overwhelmingly by the difference in entropy between U-U and B-B joints. As done for ∆ G UB , thedifference in (dimensionless) free energy ∆ G UB can be split into an average (dimensionless) energy and average entropycontribution using the same formulæ as above. If L = 0 at a certain temperature, T ∞ , then we can expect anothertype of “melting” transition, now driven by the free energy difference ∆ G UB ≡ G ( κ U ) − G ( κ B ). We return to thispoint below.It should be noticed that the above processus could in principle be carried on (with increasing difficulty) to calculatehigher order correlation functions quantities. Hence these effective Ising models give information on multi-pointtangent-tangent correlation functions. In the next section we obtain the solutions to the two effective Ising models. III. SOLUTION OF THE TWO ISING MODELS
The effective Ising partition and correlation functions can be obtained using well-known Ising transfer matrixtechniques [9]. In order to treat in parallel Z (0)I , eff and Z (1)I , eff , we introduce the index l = 0 , Z ( l )I , eff . This index will be useful in Section V where we introduce the transfer matrix approach. Weneed the eigenvectors and eigenvalues of the transfer matrices: ˆ P ( l )I , eff | ψ ( l ) i = λ l | ψ ( l ) i , where the expressions for ˆ P ( l )I , eff are given in Eqs. (12,27). The eigenvalues are λ l, ± = e J l − Γ l n cosh( L l ) ± (cid:2) sinh ( L l ) + e − J l (cid:3) / o (38)and they obey the inequality λ l, + > λ l, − . The two orthonormal eigenvectors, | ψ ( l ) i , are | l, + i = 1 √ γ l e J l (cid:0) a l | U i + a − l | B i (cid:1) and | l, −i = 1 √ γ l e J l (cid:0) a − l | U i − a l | B i (cid:1) , (39)where γ l = (cid:2) sinh ( L l ) + e − J l (cid:3) / and a l = e J l [sinh( L l ) + γ l ] / . (40)The transfer matrices can be expanded in terms of the eigenvectorsˆ P ( l )I , eff = X τ = ± λ l,τ | l, τ i h l, τ | . (41)Using the decomposition of the Ising (2 ×
2) identity matrix, ˆ I I = P τ | l, τ ih l, τ | the orthonormality of the eigenvectors | l, τ i for each value of l [i.e., h l, τ | l, τ ′ i = δ τ,τ ′ ], and the forms Eqs. (18,26) for the effective Ising partition functions,we then find Z ( l )I , eff = (4 π ) N h V | h ˆ P ( l )I , eff i N − | V i = (4 π ) N X τ = ± λ N − l,τ h V | l, τ i . (42)The full partition is given by Z = Z (0)I , eff . The matrix elements, h V | l, τ i , entering Eq. (42) can be obtained explicitlyfrom Eqs. (11), (39), and (40).When κ U = κ B = κ UB , the system decouples, but the pure Ising model with temperature-independent parameters,˜ J , ˜ K and ˜ µ , exhibits neither a second order phase transition at a finite temperature in zero “field” (˜ µ + ˜ K = 0), nora melting transition at a finite temperature for ˜ µ + ˜ K >
0. Indeed, there can be no melting transition because theinequality ϕ B < / µ + K ) ≫ e − J , and thesystem is ordered with ϕ U ≃ ϕ B ≃
0. As the temperature is raised, denaturation bubbles are thermally excited,with a cross-over when sinh( µ + K ) ≃ e − J , or roughly β ′ = ( k B T ′ ) − ≃ (˜ µ + ˜ K + 2 ˜ J ) − . At higher temperatures(sinh( µ + K ) ≪ e − J ), the average “magnetization” monotonously approaches a completely thermally disordered statewith h c i = 0 and ϕ U = ϕ B = 1 /
2. Note that this regime is not reached for DNA since as we will see below, T ′ ≃ T m and the model is certainly no longer valid for such high temperatures. By contrast, the coupled Ising-chain model willexhibit a very different behavior, with a finite temperature melting transition. In the following we implicitly assumethat all temperatures of interest obey T ≪ T ′ .Although the eigenvectors | l, τ i are orthogonal in τ for the same value of l , this is not necessarily the case for differentvalues of l (as we will see below, this is a consequence of the difference in rotational symmetry). In general, dependingon the values of the temperature-dependent effective Ising parameters, L l and J l , the eigenvectors are complicatedmixtures of the canonical basis states, | U i and | B i . These eigenvectors and their corresponding eigenvalues can,however, be simplified in two important limits: • For sufficiently low or high temperatures, below or above the transition temperature, T ∞ l (at which L l vanishes),the inequality sinh ( L l ) ≫ e − J l is obeyed (the experimental melting temperature for infinite chains is thus T ∞ m ≡ T ∞ ). As a consequence, the off-diagonal (domain-wall or “tunneling”) terms in ˆ P ( l )I , eff can be neglectedand the eigenvectors reduce asymptotically to the canonical ones, with the mapping depending on the sign of L l : | l, + i ≃ | U i and | l, −i ≃ −| B i for L l > | l, −i ≃ | U i and | l, + i ≃ | B i for L l <
0. In this limit of strongcooperativity, the eigenvalues reduce to λ l, ± ≃ exp ( J l ± | L l | − Γ l ) for sinh ( L l ) ≫ e − J l (43)and therefore the pure U state is strongly favored if L l > L l <
0, because λ l, + λ l, − ≃ exp (2 | L l | ) for sinh ( L l ) ≫ e − J l . (44)By introducing the following eigenvalues for pure U and pure Bλ l,U ≡ exp[ J + µ − G l ( κ U )] (45) λ l,B ≡ exp[ J − µ − G l ( κ B )] (46)and using the definitions of J l , L l , and Γ l , these limiting forms for the eigenvalues can be further simplified(sinh ( L l ) ≫ e − J l ): λ l, + ≃ λ l,U λ l, − ≃ λ l,B for L l > λ l, + ≃ λ l,B λ l, − ≃ λ l,U for L l > . (47) • For temperatures at or very near the transition temperature, however, the opposite inequality sinh ( L l ) ≪ e − J l holds and L l can be set to zero in Eqs. (38)–(40): the eigenvectors then reduce to symmetric and anti-symmetric,superpositions of the canonical basis vectors, | U i and | B i : | l, ±i ≃ √ | U i ± | B i ) for sinh ( L l ) ≪ e − J l (48)with eigenvalues λ l, ± ≃ e − Γ l ( e J l ± e − J l ) for sinh ( L l ) ≪ e − J l . (49)0 FIG. 2: “Energies” ε l, ± = − ln( λ l, ± ) for l = 0 , µ = 4 .
46 kJ/mol, ˜ J = 9 .
13 kJ/mol and˜ K = 0). We observe that far from the two transition temperatures, the eigenvalues reduce to the limiting forms Eq. (47). Theinset is a zoom close to T ∞ m showing the level repulsion between the branches (0 , ± ). A similar level repulsion occurs near T ∞ between the branches (1 , ± ). In this limit of weak cooperativity the eigenvalue ratio is approximately λ l, + /λ l, − ≃ coth( J l ) and the behaviorof the system is dominated by the domain walls. The average behavior for large N , which is governed by theground symmetrical state, shows vanishing average for the mean of the Ising state variable (or magnetizationin spin language) (for l = 0 and N → ∞ , ϕ U ≃ ϕ B ≃ /
2, since there is no spontaneous second order phasetransition for the 1D zero field Ising model).The exact results for the eigenvalues and eigenvectors interpolate smoothly between the above simplified resultsobtained far from and close to T ∞ l (Fig. 2).Using Eq. (11), these limiting forms can be used to obtain simple approximations for the following important matrixelements : h V | , + i = e µ/ , T < T ∞ m √ µ/ , T = T ∞ m e − µ/ , T > T ∞ m , (50) h V | , −i = − e − µ/ , T < T ∞ m √ µ/ , T = T ∞ m e µ/ , T > T ∞ m . (51)The limiting forms for h V | , −i reveal that this matrix element is negative for low temperature and positive at T ∞ m and therefore must pass through zero at a temperature T ∗ lower than T ∞ m . This special temperature will be studiedin detail in the section concerning finite size effects (Section VIII).If a melting transition exists at a finite temperature, T ∞ m , then L will go from a positive value below T ∞ m , throughzero at the transition, then to a negative value above T ∞ m . This temperature dependence for L is similar to thetime dependence of the uncoupled energy levels in the Landau-Zener problem, a quantum 2-state dipole system in anelectrical field varying linearly with time [34]. It is not surprising, therefore, that the 2 branches for the “adiabatic”states of the Landau-Zener problem are equivalent, as the time t varies from −∞ to + ∞ , to the “eigenenergies”, ε l, ± ≡ − ln( λ l, ± ), of the states | l, ±i presented above, as the temperature varies from below the transition to above(with the same type of limiting behavior near and far from the transition temperature, T = T ∞ l , equivalent to t = 0 inthe quantum Landau-Zener problem, see Fig. 2). As in the quantum mechanics of diatomic molecules, eigenenergiespossessing the same rotational symmetry (same l ) cannot cross (level repulsion), although they do reach a point ofclosest approach at T ∞ m . As an illustration we show in Fig. 2 the Landau-Zener diagrams, ε l, ± ( T ) for l = 0 , T = T ∞ m and T = T ∞ for states with same l and the possibility of level crossing for stateswith different values of l .From the full partition function, we define the (dimensionless) Helmholtz free energy per Ising variable of thecoupled system, F = − N − ln Z . The average value of the Ising state variable (or “magnetization”) can then be1obtained for finite chains using h c i = − ∂F/∂µ , from which ϕ U and ϕ B can be deduced. The expression for h c i simplifies in the limit N → ∞ , because only the largest of the eigenvalues, λ , + , entering in the l = 0 effective Isingpartition function survives: h c i → N →∞ h c i ∞ ≡ − ∂f∂µ = sinh( L )[sinh ( L ) + e − J ] / (52)where f = lim N →∞ F = − ln λ , + . Equation (52) can then be used to find ϕ U, ∞ and ϕ B, ∞ .If L vanishes at a temperature, T ∞ m , low enough for the e − J term in the denominator to be sufficiently small,then the system will undergo a sharp melting transition: h c i ∞ will jump sharply from +1 for T < T ∞ m (pure U state)to − T > T ∞ m (pure B state). The size of e − J term in Eq.(52) will determine the width of the transition region,∆ T ∞ m ≡ (cid:12)(cid:12)(cid:12)(cid:12) ∂ h c i ∞ ∂T (cid:12)(cid:12)(cid:12)(cid:12) − T ∞ m ≃ k B [ T ∞ m ] ˜ µ exp[ − J ( T ∞ m )] (53)which is exponentially small in J ( T ∞ m ) when J ( T ∞ m ) ≫ Z (1)I , eff , we can define the free energy per Ising spin, F (1) = − N − ln Z (1)I , eff . A quantity analogous to the average Ising “magnetization”, h c i , for this partition function canthen be obtained for finite chains using h c (1) i = − ∂F (1) ∂µ = h c ( t · t N ) ih t · t N i , (54)from which ϕ (1) U = (1 + h c (1) i ) / ϕ (1) B can be obtained. From the definition of h c (1) i , we see that it is a mixedcorrelation function that describes how the average system internal state, c , is correlated with its external state (chainconfiguration) via the end-end tangent-tangent correlation function. In the limit of an infinite chain h c (1) i reduces toan expression analogous to Eq. (52): h c (1) i → N →∞ h c (1) i ∞ = − ∂f (1) ∂µ = sinh( L ) (cid:2) sinh ( L ) + e − J (cid:3) / , (55)where f (1) = lim N →∞ F (1) . This limiting formula is similar to the one obtained for h c i in Eq. (52), which implies thatif L vanishes at a temperature, T ∞ , low enough for the e − J term in the denominator to be sufficiently small, thenthe system can undergo a second type of “melting” transition: h c (1) i will jump sharply from +1 for T < T ∞ (pure c (1) “U” state) to − T > T ∞ (pure c (1) “B” state). For T ∞ m < T < T ∞ , h c (1) i ≃ +1 while h c i ≃ −
1. This impliesthat in this temperature range h c ( t · t N ) i ≃ −h c ih t · t N i . This counter-intuitive result is another manifestation ofthe non-trivial coupling between internal and external degrees of freedom. The value of the e − J term determinesagain the width of the transition region:∆ T ∞ ≡ (cid:12)(cid:12)(cid:12)(cid:12) ∂ h c (1) i ∞ ∂T (cid:12)(cid:12)(cid:12)(cid:12) − T ∞ ≃ k B [ T ∞ ] ˜ µ exp[ − J ( T ∞ )] . (56) IV. ISING STATE VARIABLE–ISING AND CHAIN CORRELATION FUNCTIONS
The average value of the local Ising variable h σ i i and the 2-point h σ i + r σ i i correlation function can be calculated bystarting from the expression Eq. (9) for the partition function and using the property that an insertion of a term σ j in the sum of products is equivalent to the insertion of the Pauli matrix in the canonical basis,ˆ σ z = (cid:18) − (cid:19) (57)at the j th position in the product of transfer matrices defining the partition function, Eq. (42). This comes from h σ | ˆ σ z | σ ′ i = σδ σ,σ ′ and the equality σ j h σ j | ˆ P (0)I , eff | σ j +1 i = X σ = ± h σ j | ˆ σ z | σ ih σ | ˆ P (0)I , eff | σ j +1 i = h σ j | ˆ σ z · ˆ P (0)I , eff | σ j +1 i . (58)2We then find h σ i i = (4 π ) N Z h V | h ˆ P (0)I , eff i i − ˆ σ z h ˆ P (0)I , eff i N − i | V i (59) h σ i + r σ i i = (4 π ) N Z h V | h ˆ P (0)I , eff i i − ˆ σ z h ˆ P (0)I , eff i r ˆ σ z h ˆ P (0)I , eff i N − r − i | V i . (60)By the same method used for reducing the partition function, we finally obtain h σ i i = (4 π ) N Z X τ ,τ h V | , τ i λ i − ,τ h , τ | ˆ σ z | , τ i λ N − i ,τ h , τ | V i (61) h σ i + r σ i i = (4 π ) N Z X τ ,τ ,τ h V | , τ i λ i − ,τ h , τ | ˆ σ z | , τ i λ r ,τ h , τ | ˆ σ z | , τ i λ N − r − i ,τ h , τ | V i (62)with τ i = ± . The matrix elements appearing in the above expressions can be found explicitly using Eqs. (11), (39),and (40).The Pauli matrix ˆ σ z , which can be interpreted as a quantum mechanical dipole moment operator, is diagonal inthe canonical basis, | U i = | + 1 i and | B i = | − i [see Eq. (57)], but not necessarily in the basis that diagonalizesˆ P ( l )I , eff . Indeed, in the l = 0 basis we haveˆ σ (0) z = (cid:18) h c i ∞ (1 − h c i ∞ ) / (1 − h c i ∞ ) / −h c i ∞ (cid:19) . (63)On the one hand, the transfer matrix ˆ P ( l )I , eff mixes the canonical basis states, which explains the complicated repre-sentation of the effective Ising partition function, Eqs. (18) and (26), state variable average, Eq. (61), and correlationfunction, Eq. (62) in this basis. On the other hand, in the basis that diagonalizes the transfer matrix, the propagationbetween measurements is simple (no mixing), but now, in general, the “dipole” measurement process, correspondingto ˆ σ z , mixes such states. By directly summing Eq. (61) over i and using the matrix elements of Eq. (63), we cancalculate h c i ( N, T ): h c i ( N, T ) = h c i ∞ (cid:20) − R V R V + e ( N − /ξ I (cid:21) + 2 R V p − h c i ∞ (cid:2) − e − ( N − /ξ I (cid:3) N (cid:2) R V e − ( N − /ξ I (cid:3) (cid:2) − e − /ξ I (cid:3) , (64)where R V ≡ h V | , −ih V | , + i (65)and ξ I is the Ising correlation length ξ I = 1 / ln( λ , + /λ , − ) , , (66)the typical size of minority B (U) domains below (above) T m . Although h V | , −i can be positive or negative (andeven zero), h V | , + i is for physical reasons strictly positive, because both | V i and | , + i are linear combinations ofthe canonical basis states with strictly positive coefficients of proportionality [see Eqs. (39), (40), (50), and (51)]. Theratio R V (which can therefore be negative, zero, or positive) is thus always well defined.The above expression, Eq. (62), for h σ i + r σ i i can be interpreted, using “path integral” imagery, as a quantummechanical measurement process over an imaginary time period of N steps of duration δ . This interpretation is basedon the 1D classical Ising representation of the partition function of a 0D quantum 2 state system [35]: the transfermatrix becomes the quantum propagator, ˆ P (0)I , eff ↔ exp( − δ ˆ H/ ~ ), where ˆ H is the quantum Hamiltonian, and theeigenvalues, λ , of the transfer matrix are related to ε , the energy eigenvalues of the Hamiltonian via ε ↔ − ln λ . Ingeneral the two states for the quantum system are coupled by a non-nonzero tunneling amplitude, which correspondsto the off-diagonal (domain-wall) terms of the transfer matrix. For instance, following Eq. (62), the system is preparedin the initial state | V i and evolves N − r − i time steps under the dynamics determined by the propagator ˆ P (0)I , eff ,until a measurement of the dipole moment is performed, determined by the action of ˆ σ z . The state that comes outof the measurement then evolves r time steps until a second measurement of the dipole moment is performed. Thestate that comes out then evolves i − h σ i + r σ i i is thus the normalizedamplitude that the system returns to the initial state | V i at the end of this double measurement process.3In the limit N → ∞ , the results Eqs. (61)-(62) simplify because we keep only the leading order terms (largesteigenvalues) which sets τ = +: h σ i i → N →∞ h c i ∞ + R V (1 − h c i ∞ ) / exp[ − ( i − /ξ I ] (67) h σ i + r σ i i → N →∞ X τ ,τ h V | , τ ih V | , + i (cid:18) λ ,τ λ , + (cid:19) i − h , τ | ˆ σ z | , τ i (cid:18) λ ,τ λ , + (cid:19) r h , τ | ˆ σ z | , + i , (68)where we have used h σ i i ∞ ≡ h c i ∞ . Using the above results we obtain the limiting form for h c i ( N, T ) when N → ∞ : h c i ( N, T ) → N →∞ h c i ∞ + 2 N R V p − h c i ∞ − e − /ξ I , (69)In the double limit N, i → ∞ , meaning that we ignore the influence of end-monomers, expressions Eqs. (67) and(68) reduce to the simpler cyclic boundary condition forms: h σ i i → N,i →∞ h c i ∞ (70) h σ i + r σ i i → N,i →∞ h c i ∞ + (cid:0) − h c i ∞ (cid:1) exp( − r/ξ I ) . (71)Using the limiting values obtained earlier, Eq. (49), we find that at the melting temperature, T ∞ m , ξ I = − / ln[coth( J )]. When e − J ( T ∞ m ) ≪ ξ I ( T ∞ m ) ≃ e J ( T ∞ m ) / ≫
1. We shall see in section VII that ξ I can be extremely large, but finite, at T ∞ m , where it reaches its maximum value. Moving away from T ∞ m in bothdirections, we find that ξ I ≃ / (2 | L | ) decreases as | T − T ∞ m | increases. When ξ I ( T ∞ m ) ≫
1, the width of the thetransition, Eq. (53), can be rewritten as ∆ T ∞ m ≃ k B [ T ∞ m ] / [˜ µξ I ( T ∞ m )]. Because the system is translationally invariantin the limit N, i → ∞ and h ( σ i + r − h σ i + r i )( σ i − h σ i i ) i = h σ i + r σ i i − h σ i i , Eqs. (70)-(71) show that for |h c i ∞ | 6 = 1 theIsing correlation length measures the range of correlations between spatially separated deviations of the local Isingspin from the average value.Using the same Ising model transfer matrix techniques employed above for the partition function, we can calculatethe chain end-end tangent-tangent correlation function, h t · t N i , which is related to the effective partitions functions Z ( l )I , eff , l = 0 , h t · t N i = Z (1)I , eff Z (0)I , eff = P τ λ N − ,τ h V | , τ i P τ λ N − ,τ h V | , τ i = h V | , + i exp (cid:2) − ( N − /ξ p , + (cid:3) + h V | , −i exp (cid:2) − ( N − /ξ p , − (cid:3) h V | , + i + h V | , −i exp [ − ( N − /ξ I ] (72)where the chain persistence lengths are defined by ξ p , ± = 1 / ln( λ , + /λ , ± ) . (73)This result indicates clearly that in general h t · t N i depends on three distinct characteristic lengths: ξ I and ξ p , ± .In order to better understand the physical content of Eq. (72) (and later results), it is useful to derive simplifiedlimiting forms for the matrix elements and persistence lengths appearing therein. Using the same technique employedto obtain the expressions for h V | , ±i shown in Eqs. (50) and (51), simple limiting forms can be derived for h V | , ±i in the temperature range of experimental interest, T < T ∞ , leading to h V | , + i ≃ e µ/ and h V | , −i ≃ − e − µ/ .Using the limiting values obtained earlier for the eigenvalues, Eq. (49), we find the following limiting forms for thetwo chain persistence lengths:1 /ξ p , + ≃ (cid:26) /ξ pU , T < T ∞ m /ξ I + 1 /ξ pU , T ∞ m < T < T ∞ and 1 /ξ p , − ≃ (cid:26) /ξ I + 1 /ξ pB , T < T ∞ m /ξ pB , T ∞ m < T < T ∞ . (74)The limiting forms for ξ p , + in the range T ∞ m < T < T ∞ and ξ p , − in the range T < T ∞ m have a simple physicalexplanation: the effective persistence lengths, ξ p , ± , of minority domains (B, or − , below T m and U, or +, above T m ) tend to the typical minority domain size, ξ I , when these domains behave as rigid rods ( ξ pU ≫ ξ I for minority Udomains and ξ pB ≫ ξ I for minority B domains).Is is interesting to note that the various correlation lengths can be identified with differences between theeigenenergies appearing in the Landau-Zener diagram (Fig. 2): ξ I = 1 / ( ε , − − ε , + ), ξ p , + = 1 / ( ε , + − ε , + ), and4 ξ p , − = 1 / ( ε , − − ε , + ). Hence we already observe in this diagram that ξ I reaches its maximum at T ∞ m which is thepoint of closest approach of the branches (0 , ± ). In the limit of large N the expression for h t · t N i substantiallysimplifies and depends on only one persistence length, ξ p , + : h t · t N i → N →∞ h V | , + i h V | , + i exp (cid:2) − ( N − /ξ p , + (cid:3) . (75)The value of N for which the limiting form Eq. (75) starts to be a good approximation to Eq. (72) depends on thetemperature via the weights, h V | l, + i and the characteristic lengths, ξ I and ξ p , ± . V. FULL TRANSFER MATRIX APPROACH
To calculate the general chain tangent-tangent correlation function, h t i · t i + r i , for the coupled model, we need tointroduce the more powerful (and more abstract) transfer kernel method. This method will also shed additional lighton the origin of the effective Ising models obtained above by first integrating out the chain degrees of freedom. Tocalculate the partition and correlation functions using this method, we need to solve a spinor eigenvalue problem inorder to find the eigenfunctions and eigenvalues of the transfer kernel: ˆ P | ˆΨ i = λ | ˆΨ i , or more explicitly X σ ′ = ± Z d Ω ′ π ˆ P σ,σ ′ (Ω , Ω ′ )Ψ σ ′ (Ω ′ ) = λ Ψ σ (Ω) , (76)where | ˆΨ i = Ψ +1 (Ω) | U i + Ψ − (Ω) | B i . (77)For the pure Ising model the eigenvalues and eigenvectors can be labeled by the index τ = ± . For the purechain model with rigidity κ and no applied stretching force (like the 1D classical Heisenberg model in zero field) theeigenfunctions, ψ lm (Ω) = √ πY lm (Ω), are proportional to the spherical harmonics, Y lm (Ω), which are indexed by theinteger pair ( l, m ), with l = 0 , , . . . , + ∞ and m = − l, . . . , + l . Furthermore, the eigenvalues for the pure chain model, λ l , are indexed only by l , because they are degenerate in m : λ l = e − κ κ l (cid:18) κ ddκ (cid:19) l (cid:20) sinh( κ ) κ (cid:21) = (cid:16) π κ (cid:17) / e − κ I l +1 / ( κ ) (78)where I l +1 / ( κ ) is the modified Bessel function of the first kind. These eigenvalues take on the values λ = e − κ sinh( κ ) /κ = exp[ − G ( κ )] and λ = λ u ( κ ) = exp[ − G ( κ )] for l = 0 , l = 0 , l, m ; τ ) used for the pure Ising and pure chain models: h Ω | ˆΨ l,m ; τ i = ψ lm (Ω) | l, τ i (79)and the eigenvalues, λ l,τ by ( l, τ ) (degenerate in m ). The eigenvalues and kets, | l, τ i , which are independent of thesolid angle, Ω, must be determined by solving the eigenvalue Eq. (76). In general, the eigenvalues and eigenvectorsfor the coupled system are not, however, simply direct products of the corresponding eigenvalues and eigenfunctionsof the uncoupled Ising and chain systems. By solving the eigenvalue equation Eq. (76), we find that the kets, | l, τ i ,appearing in the full eigenspinor, and the eigenvalues λ l,τ , have already been introduced and obtained for l = 0 and 1[Eqs. (38)-(40)] in the calculation of the effective Ising partition functions, Eq. (42). If we define G l ( κ ) = − ln λ l , thenthe same formulæ, Eqs. (38)–(40), apply for the kets | l, τ i and the eigenvalues λ l,τ in the general case l = 0 , , . . . , + ∞ , τ = ± . The eigenspinors are orthonormal: h ˆΨ l,m ; τ | ˆΨ l ′ ,m ′ ; τ ′ i = h l, τ | l ′ , τ ′ i Z d Ω4 π ψ ∗ l ′ m ′ (Ω) ψ lm (Ω) = δ ll ′ δ mm ′ δ ττ ′ . (80)Once we have the eigenvalues and orthonormal eigenfunctions, we can express the transfer kernel in an abstractoperator notation as ˆ P = + ∞ X l =0 + l X m = − l X τ = ± λ l,τ | ˆΨ l,m ; τ ih ˆΨ l,m ; τ | (81)5and then use the orthonormality of the eigenspinors, as well as the decomposition of unity,ˆ I = X σ Z d Ω4 π | σ Ω ih σ Ω | = + ∞ X l =0 + l X m = − l X τ = ± | ˆΨ l,m ; τ ih ˆΨ l,m ; τ | (82)to calculate the quantities of interest in a straightforward way. As a check on the method, we can, for example,recalculate the partition function using the following expression: Z = (4 π ) N X { σ i = ± } N Y i =1 Z d Ω i π h V | σ Ω ih σ Ω | ˆ P | σ Ω i · · · h σ N − Ω N − | ˆ P | σ N Ω N ih σ N Ω N | V i , (83)or in kernel product form Z = (4 π ) N h V | ˆ P N − | V i = (4 π ) N X l,m ; τ h V | ˆΨ l,m ; τ i λ N − l,τ . (84)Because the end vector | V i contains only the rotational ground state, | ˆΨ , ± i (i.e., l = 0 , m = 0), its matrixelements with the eigenspinors of the transfer kernel simplify to h ˆΨ l,m ; τ | V i = δ l δ m h , τ | V i . (85)Inserting this expression for the matrix element into Eq. (84) leads immediately to the result, Eq. (42), obtainedpreviously for l = 0: Z = Z (0)I , eff = (4 π ) N X τ = ± h V | , τ i λ N − ,τ . (86)In order to calculate the correlation function h t i · t i + r i , we could use t i · t i + r = 13 +1 X m = − ψ ∗ m (Ω i + r ) ψ m (Ω i ) (87)which can be obtained from Eq. (2) and the definition of the spherical harmonics. Thanks, however, to rotationalsymmetry, the average value of t i · t i + r simplifies to h t i · t i + r i = 3 h t i,z · t i + r,z i = h ψ (Ω i ) ψ (Ω i + r ) i (88)where ψ (Ω) = √ θ ). The tangent-tangent correlation function can be written in a kernel product form similarto the expression for the partition function, Eq. (83), with the difference being that we must now insert the projection(or dipole) operator along the z -axis, ˆ Z = cos( θ ), related to ψ (Ω) in the j = i and j = i + r positions. This operator,which is diagonal in the canonical basis, | σ Ω i , has the following matrix elements: h σ i Ω i | ˆ Z | σ i +1 Ω i +1 i = 1 √ ψ (Ω i ) δ (Ω i +1 − Ω i ) δ σ i +1 σ i (89)from which the following equality can be established ψ (Ω i ) h σ i Ω i | ˆ P | σ i +1 Ω i +1 i = √ X σ Z d Ω4 π h σ i Ω i | ˆ Z | σ Ω ih σ Ω | ˆ P | σ i +1 Ω i +1 i = √ h σ i Ω i | ˆ Z · ˆ P | σ i +1 Ω i +1 i . (90)In operator product form, using Eqs. (83,87) the correlation function then becomes h t i · t i + r i = 3(4 π ) N Z − h V | ˆ P i − ˆ Z ˆ P r ˆ Z ˆ P N − r − i | V i . (91)In the basis that diagonalizes the transfer kernel, ˆ P , the operator ˆ Z , which is not diagonal, has the following matrixelements: h ˆΨ l ′ ,m ′ ; τ ′ | ˆ Z | ˆΨ l,m ; τ i = h l ′ , τ ′ | l, τ i δ mm ′ " δ l ′ ,l − (cid:18) l − m l − (cid:19) / + δ l ′ ,l +1 (cid:18) ( l + 1) − m l + 1) − (cid:19) / (92)6which, aside from the factor h l ′ , τ ′ | l, τ i , is the well known selection rule for quantum dipole transitions, i.e., ∆ l = ± m = 0 (in, for example, the Stark effect [36]). Although h l, τ ′ | l, τ i = δ ττ ′ the matrix element h l ′ , τ ′ | l, τ i is notnecessarily zero for l = l ′ and τ = τ ′ , because in this case the matrix element is between states of different rotationalsymmetry. This result indicates that the measurement of the z -axis dipole moment can also induce a change in theinternal state, τ , of the system.Equation (85) shows that in the expression for h t i · t i + r i , the ˆ Z projection operator can only connect a ( l = 0 , m = 0)rotational state with an ( l = 1 , m = 0) one, or vice-versa (as in the Stark effect for the 1s state of the hydrogen atom).To evaluate h t i · t i + r i using Eq. (91) we therefore only need one matrix element, h ˆΨ , τ ′ | ˆ Z | ˆΨ , τ i = 1 √ h , τ ′ | , τ i . (93)By inserting the decomposition of unity between each matrix factor in Eq. (91) and using the orthonormality of theeigenspinors, we obtain h t i · t i + r i = (4 π ) N Z X τ ,τ ,τ h V | , τ i λ i − ,τ h , τ | , τ i λ r ,τ h , τ | , τ i λ N − r − i ,τ h , τ | V i . (94)When i = 1 and r = N −
1, using again the decomposition of unity in the | l, τ i space, we recover our previous resultfor h t · t N i , Eq. (72).The above expressions for h t i · t i + r i , Eqs. (91)–(94) can also be interpreted, using the “path integral” representationof a quantum statistical partition function, as a quantum mechanical measurement process over an imaginary timeperiod of N steps. This interpretation is based on the (
1D classical Ising representation ) ⊗ (
1D classical Heisenberg )representation of the partition function of a 0D quantum diatomic molecule, modeled as a 2 state rigid rotator. Thesystem is prepared in an initial state | V i that is in a mixture of the internal states, τ = ± , but in the sphericallysymmetric rotational ground state. This initial state evolves N − r − i time steps under the dynamics determinedby the propagator ˆ P , until a measurement of the dipole moment along the z -axis, determined by the action of ˆ Z , isperformed. The projections of | V i onto the rotational ground states ( l = 0 , m = 0), h , ±| V i , evolve in a simple waybecause these states are associated with eigenspinors of the transfer kernel (each time step gives rise to an eigenvaluefactor, λ ,τ ). This dipole measurement excites the system from the rotational ground state to the ( l = 1 , m = 0),non-spherically symmetric, first rotational excited state, along with a possible transition in the internal state ofthe diatomic molecule. The state that emerges from the measurement then evolves r time steps, until a secondmeasurement of the dipole moment is performed. This measurement de-excites the system from the ( l = 1 , m = 0)excited state back down to the ( l = 0 , m = 0) ground state or up to the ( l = 2 , m = 0) excited state [see Eq. (92)]again with a possible change in internal state. The state that comes out of this second measurement then evolves i − h t i · t i + r i is then thenormalized amplitude that the system returns to the initial state | V i at the end of this double dipole measurementprocess [because the final state and the ( l = 2 , m = 0) excited state are orthogonal, no l = 2 matrix elements appearin Eq. (94)]. By resolving the initial and final states, both equal to | V i , into | , ±i components, the sum in Eq. (94)is over all possible “time” sequences involving | , ±i .In order to study the limiting forms of h t i · t i + r i , we write h t i · t i + r i = P τ ,τ ,τ C ( τ , τ , τ ) (cid:16) λ ,τ λ ,τ (cid:17) i − (cid:16) λ ,τ λ ,τ (cid:17) r (cid:16) λ ,τ λ , + (cid:17) N − P τ (cid:16) λ ,τ λ , + (cid:17) N − h V | , τ i , (95)where we have introduced the joint amplitude C ( τ , τ , τ ) = h V | , τ ih , τ | , τ ih , τ | , τ ih , τ | V i . (96)When all the bending rigidities are equal to κ , we recover, as expected, the pure discrete wormlike chain result, h t i · t i + r i = exp[ − r/ξ p ( κ )] with ξ p given in Eq. (32).In the limit N → ∞ , we keep only the leading order term ∝ λ N − , + in Z and the surviving terms in the numerator2-point correlation function ( τ = +) to find: h t i · t i + r i → N →∞ X τ ,τ C ′ ( τ , τ ) (cid:18) λ ,τ λ , + (cid:19) i − (cid:18) λ ,τ λ , + (cid:19) r (97)7where C ′ ( τ , τ ) = h V | , τ ih , τ | , τ ih , τ | , + ih , + | V i . (98)Using the effective chain persistence lengths introduced previously in Eq. (73), we can express Eq. (97) in a physicallymore transparent form: h t i · t i + r i → N →∞ X τ = ± exp (cid:2) − r/ξ p ,τ (cid:3) { C ′ (+ , τ ) + C ′ ( − , τ ) exp [ − ( i − /ξ I ] } (99)with ξ p , − < ξ p , + and ξ I the Ising correlation length already introduced in Eq. (66).In the double limit N, i → ∞ , the dependence on the chain ends disappears again and the expression Eq. (99)simplifies to h t i · t i + r i → N,i →∞ X τ = ± h , τ | , + i exp (cid:0) − r/ξ p ,τ (cid:1) , (100)which reveals the importance of the two persistence lengths, ξ p , ± , and the two “transition probabilities”, h , ±| , + i ,for going from the ground state | , + i to the first rotational excited state, | , ±i , with or without a change in internalstate τ . In the temperature range of experimental interest, T < T ∞ , h , + | , + i and h , −| , + i can, to an excellentapproximation, be set equal to ϕ U, ∞ and ϕ B, ∞ , respectively. When this last result is used in conjunction with thelimiting forms for ξ p , ± , Eq. (74), a useful approximation is obtained for Eq. (100), valid for T < T ∞ : h t i · t i + r i ≃ N,i →∞ ϕ U, ∞ exp (cid:0) − r/ξ p , + (cid:1) + ϕ B, ∞ exp (cid:0) − r/ξ p , − (cid:1) . (101)For N, i → ∞ and short distances, r ≪ ξ p , − , we find the limiting linear behavior in r : h t i · t i + r i ≃ N,i →∞ − r/ξ p eff , CF , (102)where 1 /ξ p eff , CF ≡ h , + | , + i /ξ p , + + h , −| , + i /ξ p , − ≃ ϕ U, ∞ /ξ p , + + ϕ B, ∞ /ξ p , − (103)is an effective persistence length for the correlation function (CF) at short distances that clearly reveals the importanceof the shortest persistence length, ξ p , − , under these conditions. In the triple limit N, i, r → ∞ , only one term survives: h t i · t i + r i → N,i,r →∞ h , + | , + i exp (cid:0) − r/ξ p , + (cid:1) . (104)Exactly at T ∞ m , the above expression Eq. (100) simplifies to h t i · t i + r i T ∞ m ≃ N,i →∞
12 [exp ( − r/ξ pU ) + exp ( − r/ξ pB )] , (105)which for r ≪ ξ pB reduces to h t i · t i + r i T ∞ m ≃ − r /ξ pU + 1 /ξ pB ) . (106)Because the inverse persistence lengths enter into Eq. (102), the short distance limiting behavior will tend to bedominated by the shortest one, ξ p , − above T ∞ m , where h , −| , + i > h , + | , + i . Below T ∞ m , however, there will bea competition between the weights h , −| , + i < h , + | , + i and the persistence lengths, 1 /ξ p , − > /ξ p , + .The conditions under which these limiting expressions are valid approximations depends critically on the weightsappearing in the above expressions. If, for example, h , + | , + i ≃ ϕ U, ∞ is sufficiently small compared with h , −| , + i ≃ ϕ B, ∞ at a certain temperature, then the “subdominant” term in Eq. (100), ∝ h , −| , + i (possessingthe smaller persistence length) may actually be dominant over a wide range of r values, as shown in Fig. 3.Indeed, due to the coupling between bending and internal states of DNA, for realistic parameter values (cf. Sec-tion VII), the respective weights h , + | , ±i associated with each correlation length change abruptly at T ∞ m : below T ∞ m , we have h , + | , + i ≃ h U | U i = 1 and h , + | , −i ≃ h U | B i = 0, thus ξ p , + ≃ ξ pU . For T ∞ m < T < T ∞ , we find8 FIG. 3: (a) Variation of the transition probabilities h , ±| , + i with temperature (red solid line for + and dashed blue linefor − ). For T < T ∞ , h , + | , + i ≃ ϕ U, ∞ and h , −| , + i ≃ ϕ B, ∞ . One observes that for T ∞ m < T < T ∞ , h , −| , + i = 1which underlines the relevance of ξ p , − in this temperature range (with parameter values ˜ µ = 4 .
46 kJ/mol, ˜ J = 9 .
13 kJ/mol and˜ K = 0, see section VII). (b) Tangent-tangent correlation function given by Eq. (100) ( N, i → ∞ ) for 3 different temperatures:just before the transition (dotted blue line), controlled by ξ p , + ≃ [1 /ξ pU + 1 /ξ I ] − ≃ ξ pU ; slightly above T ∞ m (solid red line),where the correlation length ξ p , − ≃ ξ pB is dominant for r < r ∗ ≃
20; and after the transition ( T ∞ m < T < T ∞ ) where thecorrelation length ξ p , + ≃ ξ I disappears in favor of ξ p , − ≃ ξ pB (dashed green line).FIG. 4: Variation with temperature of the various correlation lengths appearing in the model results: the Ising correlationlength, ξ I (solid red line); persistence lengths of the coupled system, ξ p eff (appearing for long chains) and ξ p , ± (dashed dottedline respectively blue and green); and of the pure chains, ξ pU,B (dashed purple lines) (in units of a ). At T ∞ m , the Ising correlationlength is peaked but finite, which is related to the point of closest approach of the two (0 , ± ) branches (zoom of Figure 2), andthe effective persistence length, ξ p eff , rapidly crosses over from ξ pU to ξ pB (parameter values ˜ µ = 4 .
46 kJ/mol, ˜ J = 9 .
13 kJ/moland ˜ K = 0). h , + | , + i ≃ h B | U i = 0 and h , + | , −i ≃ h B | B i = 1 which implies ξ p , − ≃ ξ pB . For higher temperatures, T > T ∞ ,the respective weights get swapped again, but now ξ p , + ≃ ξ p , − . These considerations lead us to introduce a criticaldistance, r ∗ , at which the two terms in Eq. (100) are equal: r ∗ ≡ ξ p , − − ξ p , + ! − ln (cid:18) h , + | , + i h , −| , + i (cid:19) ≃ ξ pB ln (cid:18) ϕ B, ∞ ϕ U, ∞ (cid:19) , (107)where in arriving at the last approximation we have used limiting forms that are valid when T < T ∞ and assumed that ξ pB ≪ ξ pU , ξ I (cf. Fig. 4). When r < r ∗ then the correlation function Eq. (100) is dominated by the shortest persistencelength, ξ p , − , and when r > r ∗ the correlation function is dominated by the longest one, ξ p , + . For sufficiently longchains and temperatures close enough to T m , the inequality N ≫ r ∗ holds and this cross over should be clearly visible(see Fig. 3).9 VI. MEAN-SQUARE END-TO-END DISTANCE
To calculate the mean-square end-to-end distance of the chain, we use the two-point correlation function obtainedabove: (cid:18) Ra (cid:19) = N X i,j =1 h t i · t j i = N + 2 N − X i =1 N − i X r =1 h t i · t i + r i . (108)When all the bending rigidities are equal to κ , we recover, as expected, the pure discrete wormlike chain result: R ( κ ) = a N W N ( u ( κ )) where W N ( z ) = 1 + z − z − zN − z N (1 − z ) (109)with u ( κ ) defined in Eq.(31). In the limit N → ∞ , R ( κ ) → N →∞ a N e − /ξ p ( κ ) − e − /ξ p ( κ ) . (110)More specifically, three distinct regimes can be identified: R ( κ ) → a N, a N ξ p ( κ ) ,a N , ξ p ( κ ) ≪ ≪ N ≪ ξ p ( κ ) ≪ N < N ≪ ξ p ( κ ) (freely jointed Gaussian)(effective spin wave Gaussian)(rigid) . (111)Because ξ p ( κ ) = ξ p ( β ˜ κ ) is decreasing function of temperature, the pure DWLC will go from the rigid to the effectiveGaussian to the freely jointed Gaussian regime as the temperature is raised.For the coupled model the double summation in Eq.(108) can also be carried out and we find (cid:18) Ra (cid:19) = N + 2 P τ ,τ ,τ C ( τ , τ , τ ) S N ( λ ,τ , λ ,τ , λ ,τ ) (cid:16) λ ,τ λ , + (cid:17) N − P τ (cid:16) λ ,τ λ , + (cid:17) N − h V | , τ i , (112)where S N ( x, y, z ) = ( N yx − y − yx − ( y/x ) N (1 − y/x ) , for x = zz − N T N ( y,z ) − T N ( y,x ) z − x , for x = z , (113)with T N ( x, y ) = y N x/y − ( x/y ) N − x/y . In the limit N → ∞ the above complicated expression for R simplifies to an effectiveGaussian form R → N →∞ a N ξ p eff where ξ p eff ≡ X τ h , τ | , + i e − /ξ p ,τ − e − /ξ p ,τ (114)is an effective ”long chain” persistence length. This expression can be also obtained by using the simplified N, i → ∞ limiting form for h t i · t i + r i , Eq. (100), in the general formula for ( R/a ) , Eq. (108). It tends to a further limitingform when ξ p , ± ≫ ξ p eff → X τ h , τ | , + i ξ p ,τ ≃ (cid:26) ϕ U, ∞ ξ pU + ϕ B, ∞ ξ pB , T < T ∞ m ϕ U, ∞ [1 /ξ pU + 1 /ξ I ] − + ϕ B, ∞ ξ pB , T ∞ m < T < T ∞ (115)At T ∞ m this expression simplifies reduces to ξ p eff ≃ ( ξ pU + ξ pB ) / ξ I ≫ ξ pU , which is actually the case (Fig. 4).For T < T ∞ m , the longest persistence length dominates: ξ p eff ≃ ξ p , + ≃ ξ pU . Above T ∞ m , however, we see onceagain that there may be a competition between the persistence lengths, ξ p , ± , and the “transition probabilities”, h , ±| , + i , appearing in Eq. (115). This competition occurs now for T > T ∞ m , contrary to what was found for theshort distance behavior of the 2-point correlation function, Eq. (102), because the persistence lengths themselvesappear in Eq. (115), and not their inverses. For T ∞ m < T < T ∞ , depending on the weights, the smaller persistencelength, ξ p , − , may actually be dominant over the larger one, ξ p , + ; if so, ξ p eff ≃ ξ p , − ≃ ξ pB , which is actually the casewhen we consider standard parameter values for dsDNA and ssDNA (section VII), as shown in Fig. 4.0 VII. APPLICATION TO SYNTHETIC DNA THERMAL DENATURATION
Melting or thermal denaturation profiles are experimentally obtained by following the UV absorbance of a DNAsolution while slowly increasing the sample temperature. This method allows one to follow the temperature evolution ofthe fraction of base-pairs that have been disrupted, ϕ B ( T ). A typical profile has a sigmoid shape possibly with bumpsthat could appear depending on the DNA sequence. Different Ising-type models have been proposed [4, 5, 13, 14] formodeling denaturation curves by focusing on the influence of the base pair sequence, but they do not attempt to takeinto account properly the fluctuations of the DNA chains themselves. Yet, chain fluctuations increase with T andplay a crucial role in determining melting profiles. Moreover, these fluctuations concern both stiff helical segmentsand flexible coils with different bending rigidities.In this section, we compare the model developed above, whose key element is to account for internal state fluctuationson an equal footing with those of the chain, with a set of experimental data. We focus on the evolution of ϕ B ( T ) for asynthetic homopolynucleotide polydA-polydT. Six independent parameters appear in the theory: the polymerizationindex N , the three Ising parameters K , J and µ defined in Eq. (1) and Fig. 1, and bending moduli κ U for dsDNA and κ B for ssDNA. Note that we have also introduced a bending rigidity κ UB for domain walls. However κ UB appears inthe theory only in the effective cooperativity parameter J . Thus changing κ UB is equivalent to varying the bare J ,i.e. the energetic penalty to create a wall. Without any lost of generality, we choose to fix κ UB = κ U . Moreover, wechoose free boundary conditions for the end monomers, which is valid unless their state is fixed by the experimentalconditions [37] (although any type of boundary conditions can be treated using our model). Of the six parameters,three are determined experimentally: N , κ U , and κ B . Moreover, there is evidence that stacking interactions in dsDNAand ssDNA are of the same magnitude which justifies the choice of ˜ K = 0 adopted below [38].Figure 5 shows ϕ B ( T ) for a polydA-polydT of molecular weight M w = 1180 kDa in a solution of 0.1 SSC (0.015 MNaCl + 0.0015 M sodium citrate, pH 7.0) taken from [4]. FIG. 5: Fraction of broken base-pairs for a polydA-polydT vs. temperature (solution of 0.1 SSC, N = 1815). The solid linerepresents the theoretical law for µ = 1 . k B T m and J = 3 . k B T m where T m = 326 . N → ∞ for the sameparameter values corresponds to the broken line. In order to compare the data with our model predictions, we choose the experimental values persistence lengths, ℓ pds ≃
50 nm and ℓ pss ≃ κ U = ℓ pds /a = 147 and κ B = 2 ℓ pss /a = 5 .
54 at T = T m (taking a = 0 .
34 nm for one base-pair size and a factor of 2 for two flexible segments in parallel per coil segment). The tworemaining parameters ˜ µ and ˜ J are determined by fitting the experimental data. The solid line in Fig. 5 correspondsto ˜ µ = 1 . k B T m ≃ .
46 kJ/mol and ˜ J = 3 . k B T m ≃ .
13 kJ/mol leading to T m = 326 . µ , wefind a value very close to the experimental value of 10.5 kJ/mol [2]. Although the value of ˜ J is more difficult tointerpret, our result ˜ J ∼ µ is consistent with the idea that stacking interactions make the dominant contributionto DNA stability [5]. Chain fluctuations do not only renormalize the effective free energy, 2 ˜ L , required to break aninterior base-pair, but also the cooperativity parameter ˜ J : the latter varies almost linearly with T following Eq. (15)contrary to previous theories where ˜ J was taken as constant and supposed to be purely enthalpic in character [4]. Wehave for the total cooperativity parameter ˜ J = 4 . k B T m at T = T m , which shows that the bending contribution isroughly 25% (remembering that we have chosen κ UB = κ U ). The model fit thus leads to parameter values in accordwith experiment. Our model predictions for experimentally accessible A-T pair quantities are also in agreement with1accepted values [3, 39]: i) the loop initiation, factor, σ LI ≡ e − J ≃ − at T m ; and at physiological temperature, T ph , ii) an interior single base-pair opening probability ϕ B ( T ph ) ≃ − with a bubble initiation barrier of 17 k B T ,and iii) a free energy of 0 . k B T for breaking an additional base-pair in an already existing bubble. In reality, thefitted values of ˜ µ and ˜ J implicitly compensate for effects like loop entropy explicitly left out of the model [4]. Asshown in Fig. 9 of [4], effective Ising models without loop entropy, like ours, can be considered to account implicitly(and approximately) for loop entropy, provided that one allows for loop entropy contributions to both J and K . Thisloop entropy renormalization of the Ising model parameters will depend on the value of the loop entropy exponent k and the chain length N and could allow for a simple approximate way of accounting for the influence of loop entropywithin the framework of an effective Ising model (cf. [40]). This renormalization probably explains why our the modelvalue for σ LI is at the low end of the accepted spectrum.In Fig. 5, the curve corresponding to the thermodynamic limit ( N → ∞ ) is shown for the same parameter values.In this case, ϕ B ( T ) is given by Eq. (52) and the value of the melting temperature is obtained analytically as a functionof T ∞ m using L ( T ∞ m ) = 0 which is given in the limit of low temperature (˜ κ B ≫ k B T ) by k B T ∞ m ≃ µ + ˜ K ln(˜ κ U / ˜ κ B ) (116)Hence, the melting temperature is reached when the enthalpy required to create a link is perfectly balanced bythe difference in (entropy dominated) free energy between the two types of semi-flexible chains (U or B). Anotherquantity which has an experimental relevance is the width of the transition. In the thermodynamic limit this widthis narrow, but nonzero, due to the large but finite cooperativity parameter: ∆ T ∞ m ∝ /ξ I ( T ∞ m ) ≃ − J ( T ∞ m )][see Eq. (53)]. Hence, the thermodynamic limit clarifies the role of the two free model parameters: in conjunctionwith the experimentally known bending rigidities, µ sets the melting temperature and J fixes the transition width.This is in contrast to previous Ising-like models, where three fitting parameters were used J , ∂L /∂T , and T ∞ m with L assumed to by a linear function of T [4].Within the scope of our model the measured transition width is indicative of a very long Ising correlation length, ξ I , near the transition temperature, much larger than the pure U and B persistence lengths; therefore typical helix(U) and bubble (B) domains (of size ∼ ξ I ) are flexible within a small temperature window near the transition.Included in the predictions of our theory are mechanical and structural features of the composed chain, such aspersistence length or mean square end-to-end radius, R . This differs from purely Ising-type models [4, 5] and non-linear microscopic models [20, 21] where only thermodynamical quantities related to base-pairing are available. Thevariation of the effective persistence length ξ p eff (and thus the radius of gyration for long chains) vs. T is shown inFig. 4. It varies from ξ pU for T < T ∞ m to ξ pB for T > T ∞ m . Since the transition is very abrupt, we suggest that thedenaturation transition can also be followed experimentally by measuring directly the radius of gyration, for instanceby tethered particle motion [41], light scattering, or viscosity experiments. For instance, since the relative viscosity isproportional to c DNA R (where c DNA is the DNA concentration), it should clearly exhibit an abrupt thermal transitionat a given c DNA and N . Such a transition has indeed been observed for the viscosity of synthetic homopolynucleotidesolutions [42], in qualitative agreement with Fig. 4.In fitting our model to experiment for chains of length N = 1815, we have found that finite size effects play animportant role (Fig. 5). In the following section, we investigate such effects in detail. VIII. FINITE SIZE EFFECTS
It has been shown experimentally that DNA thermal denaturation varies with chain length, N [10]. In this section,we carefully study the effect of chain ends on the denaturation transition. In Fig. 6 are shown denaturation profilesfor various chain lengths from N = 100 to N → ∞ and the fixed parameters values ˜ µ = 4 .
46 kJ/mol, ˜ J = 9 .
13 kJ/moland ˜ K = 0 used in the previous section to fit the melting data. Within the scope of our model one observes thati) the melting temperature T m ( N ) is a decreasing function of N , varying as [ T m ( N ) − T ∞ m ] /T ∞ m ≃ / ( N − T ∗ at which ϕ B ≃ .
03; iii) the transition width ∆ T m ( N ) is adecreasing function of N .Concerning points (i) and (ii), the observed behavior for the coupled system with the present parameter values isdirectly related to the model result that T ∗ < T ∞ m , which is radically different from the behavior found for the simpleIsing model (for which melting curves, ϕ B , are strictly decreasing functions of N , because formally, T ∗ = ∞ when κ U = κ B = κ UB ). The present behavior for melting maps is also very different from the predictions of older empiricalIsing-like models of denaturation and Helix-Coil like transitions [9], for which the chemical potential µ appearing inthe end vector | V i is incorrectly identified with L . This identification results in melting curves independent of N ,i.e., T m = T ∗ .2 FIG. 6: Linear-Log plot of melting curves for N = 100 (dashed dotted green line), 500 (dotted purple line), 1815 (solid blue line), ∞ (dashed red line) in decreasing order at low temperature, T <
326 K (parameter values ˜ µ = 4 .
46 kJ/mol, ˜ J = 9 .
13 kJ/moland ˜ K = 0). We note that all the curves intersect at T ∗ , as discussed in the text. T m is defined by ϕ B = 0 .
5. Inset: Log-Logplot of model results for the shift in transition width ∆ T m − ∆ T ∞ m vs. polymer length. Dots correspond to the model resultsand the solid line is a law in 1 /N . Concerning point (iii), the transition width roughly follows the law (∆ T m − ∆ T ∞ m ) ∼ / ( N − N ∼ , finite size effects are important. For very short chains, e.g., N < N = 10. This point is crucial, since it has beenobserved experimentally that for polydA-polydT inserts between more stable G-C rich domains, melting curves aremuch wider for very short DNA chains ( N ∼
10 bp) [37] with a width that decreases with increasing N (observed for60 < N <
140 in [10]). In such experiments, the nature of end monomers clearly becomes extremely important.For a given N and T the local site dependent bubble opening probability, or melting map, ϕ B,i = 1 − h σ i i , (117)can be obtained from h σ i i given in Eq. (61). In Fig. 7 ϕ B,i is plotted for six different temperatures using the samemodel parameter values employed in Fig. 5; we observe that below T ∗ the chain unwinds from the ends, whereasabove this temperature an interior bond has a higher probability of being open than an end one. Far enough below T ∗ the melting curve heals rapidly to a plateau value close to ϕ B, ∞ on a length scale on the order of ξ I ≪ N . Atphysiological temperature, 310 K, the interior bond opening probability is ≃ − in agreement with the experimentalvalue for long runs of A-T pairs (which is an order of magnitude lower than randomly placed A-T pairs) [3]. At thistemperature the end bonds have opening probabilities two orders of magnitude greater than the interior ones.At T ∗ the melting curve is perfectly flat. This result, which is independent of chain length N , indicates thateach Ising variable can be considered to fluctuate independently: h σ i i is constant independent of i , despite a 2-point correlation function that does not factorize [ h σ i σ i + r i 6 = h σ i ih σ i + r i , cf. Eq. (71)] and a large Ising correlationlength ( ξ I ≫
1, see Fig. 4). Indeed, the influence of the renormalized stacking energy ( ∼ ˜ K < ∼ ˜ J ), which suppresses bubbleformation, and therefore ϕ B,i ( T ∗ ) = [1 − tanh(˜ µ/ ( k B T ∗ )] / ≃ .
03, which results from taking N = 1 or taking K = J = 0 for arbitrary N . In some ways this compensation leads to an effective non-interacting Ising system with( ∂ h c i /∂N ) T ∗ = ( ∂ϕ B /∂N ) T ∗ = 0, which explains why the melting curves for different values of N cross at T ∗ inFig. 6. Below T ∗ the combined effects of the renormalized stacking energy and entropy gain favoring interior bubblesare too small to overcome the destacking energy cost associated with an extra domain wall and the chain ends unwindfirst. Since K becomes more negative with increasing T faster than J increases, a temperature T ∗ is reached wherethe renormalized stacking and destacking effects just compensate. Above T ∗ the situation is reversed and the openingprobability is higher in the chain interior. For arbitrary N and T , the thermodynamic chemical potential, defined byˆ µ = ( ∂F/∂N ) T can be related to h c i via Maxwell-type relations, leading to ( ∂ ˆ µ/∂ ˜ µ ) T,N = N ( ∂ h c i /∂N ) T − h c i . At T ∗ this general relation simplifies to ( ∂ ˆ µ/∂ ˜ µ ) T ∗ ,N = −h c i| T ∗ = ( ∂f /∂ ˜ µ ) T ∗ ,N , characteristic of a non-interacting system.Upon examination of Eqs. (61) and (86), we see that T ∗ is determined by the condition that h σ i i = h c i ∞ , whichis obtained when the end vector | V i is identical to the eigenket | , + i and orthogonal to | , −i : h V | , + i = 1 and3 FIG. 7: Average melting maps for different temperatures for N = 1815 and the same parameter values used in Fig. 5: plotof the fraction of broken bases, ϕ B,i as a function of the base position i for, in increasing order, 310 K, 0 . T ∗ = 322 .
87 K, T ∗ = 326 .
13 K, T ∞ m = 326 .
24 K, T m = 326 . . T m = 329 .
66 K. h V | , −i = 0. Physically, this means that the coupled Ising-chain system is in a pure state, | , + i , that mixes thecanonical states in a special way. The temperature T ∗ can be obtained by solving h V | , −i = 0. Using Eqs. (11)and (39)-(40) this translates into e µ = e − J { sinh( L ) + [sinh ( L ) + e − J ] / } . After some manipulation usingEq. (52), this can be shown to be identical to h c i ∞ ( T ∗ ) = h c i ( N = 1 , T ∗ ) = tanh[˜ µ/ ( k B T ∗ )]. Furthermore, Eqs. (62)and (94) show that the Ising and chain 2-point correlation functions, h σ i · σ i + r i and h t i · t i + r i , also get simplified at T ∗ : the approximate forms, Eqs. (71) and (100) valid in general only for N, i → ∞ , become exact for arbitrary N and i at this special temperature. It is clear that at T ∗ the coupled system behaves as if there are no end effects andfinite chains have the same behavior as an infinite one.To shed additional light on this mechanism and illustrate the important role of internal bubble entropy for longchains, we now study an infinite chain and compare ϕ B, int = lim N,i →∞ ϕ B,i with ϕ B, end = lim N →∞ ϕ B, . Equa-tion (67) shows that ξ I plays here the role of a healing length , over which end effects relax (see Fig. 7). The ratio ofmatrix elements appearing in Eq. (67) gets simplified in the following way for special values of T : R V = h V | , −ih V | , + i = − e − µ , T < T ∗ , T = T ∗ tanh( µ/ , T = T ∞ m e µ , T > T ∞ m (118)At T ∞ m , Eq. (67) simplifies to h σ i i ∞ ( T ∞ m ) = tanh (cid:18) ˜ µ k B T ∞ m (cid:19) exp[ − ( i − /ξ I ( T ∞ m )] , (119)which shows that for very long chains ( N ≫ ξ I ) at T ∞ m ( > T ∗ ) ϕ B, end = [1 − tanh[˜ µ/ (2 k B T ∞ m )] / ≃ . ≪ ϕ B, ∞ ( T ∞ m ) = 1 /
2, revealing an internal opening probability more than three times higher than an end one. InFigs. 7 and 8, however, we observe that for T = T m and T ∞ m , ξ I > N = 1815 (cf. Fig. 4), and therefore end effects donot get damped out near the center of this finite chain. Indeed, at T ∞ m the opening probability ϕ B,i near the middleof a chain of length N = 1815 is much less than the value of 1/2 holding for an infinite one. For N = 1815 we stillobserve, however, noticeable differences ( ∼
10 to 20%) between interior and end opening probabilities near T m .The one-sequence-approximation has been defined by Poland and Scheraga [14] and consists in neglecting the manysmall bubbles which eventually collapse and considering only one large thermally excited bubble. It is valid fortemperatures sufficiently far below T ∞ m . For N → ∞ , ϕ B, int and ϕ B, end can be estimated in this approximation bysumming over, respectively, all the interior and end bubbles containing the fixed site in question. In the interior casewe thus have ϕ B, int ≃ ∞ X n =1 n exp h − β ∆ G ( n )int i = e − J ∞ X n =1 n e − nL = e − J ( L ) . (120)4 FIG. 8: Zoom of the melting map shown in Fig. 7 for T ∞ m (lower curve) and T m (upper curve). The factor of n in the sum is entropic in nature and equal to the number of ways of placing a fixed interior site withinan interior n -bubble. In the end case, where there is no entropic factor, ϕ B, end ≃ ∞ X n =1 exp h − β ∆ G ( n )end i = e K − J ∞ X n =1 e − nL = e − (2 J + µ ) L ) . (121)It is important to note that for the model parameters employed, the additional dimensionless free energy for breakingan additional base-pair in an already existing bubble, 2 L ≃ .
18 at physiological temperature, is less than one andmuch smaller than the bubble initiation energy cost for an end bubble, ∼ J ≃ a fortiori for an interior one, ∼ J ). This implies that even at physiological temperature, where the probability of bond opening is very small,bubbles covering a wide range of sizes contribute to the sums in Eqs. (120)-(121): both h n i int = coth( L ) ≃ /L ≃ h n i end = e − L / [2 sinh( L )] ≃ / (2 | L | ) ≃ ξ I , which varies as 1 / (2 | L | )when | L | ≪
1. Because L decreases with T (going to 0 at T ∞ m ), both ϕ B, int and ϕ B, end are increasing functions oftemperature. For T < T ∗ the bubble initiation energy cost dominates and ϕ B, int < ϕ B, end . Thanks to the entropicfactor, however, ϕ B, int increases more rapidly than ϕ B, end and the two curves cross over at T ∗ , an estimate of whichcan be obtained by equating the above two one-sequence-approximations for both these quantities. As a check on thepreceding discussion, the one-sequence-approximations for infinite chains obtained above for end and interior openingprobabilities can be shown to be in exact agreement with what one gets from the definition of ϕ B,i , Eq. (117), andEq. (67) by using the low temperature approximations for h V | , −i / h V | , + i ≃ − e − µ [Eq. (118)] and h c i ∞ , Eq. (52),i.e. expanding to lowest order in e − J , assuming that e − J ≪ sinh ( L ), which is the formal criterion for the validityof the one-sequence-approximation.Using the one-sequence-approximation, it is now easy to see how the cost in loop entropy associated with internalbubbles will modify the above results. In the so-called loop entropy models [4, 14, 14] the internal bubbles formedby the single strands are visualized as one polymer loop, whose entropic cost has been estimated first by Zimm [15].Hence, although ϕ B, end does not change, the interior opening probability does, becoming ϕ LE B, int ≃ ∞ X n =1 n/ ( n + n ) k exp[ − β ∆ G ( n )int ] , (122)where we have adopted a common simplified form for the loop entropy factor, f LE ( n ) = ( n + n ) − k parametrized by aconstant n , which may be as large as 100 [10], and an exponent k , usually assumed to be in the range 3 / ≤ k ≤ . n would severelyreduce the importance of loop entropy for short chains and probably reflects the presence of strong bending rigidityeffect (an important open question concerns how to incorporate bending rigidity into f LE in a physically correctway). Loop entropy clearly lowers the probability of interior bubble opening and will lead to an increase in T ∗ .The situation is further complicated if strand sliding, which may be important for periodic DNA, is taken intoaccount. For homopolymeric DNA like polydA-polydT, strand sliding leads to a modified loop entropy exponent, k ′ = k − k appears to be greater than 2 when self-avoidance is fully taken into account [18].The combined effects of loop entropy and strand sliding will lead to increases in both T ∗ and T ∞ m . Neither the5theoretical nor the experimental situation concerning T m ( N ) is entirely clear for homopolymeric DNA with free endsand further careful experiments are clearly called for. If loop entropy and strand sliding are included in the coupledIsing-chain model, T ∗ might become higher than T ∞ m , which would imply that T m ( N ) would increase with N [14],unlike what what we find to occur when these two effects are neglected.In the future, we intend to examine these questions by incorporating loop entropy and strand sliding directly intoour DNA model. With all other factors being equal, adding loop entropy and strand sliding increases the stability ofthe closed state, resulting in sharper melting curves and higher values of T m ( N ) (see Figs. 9 and 10 of [4] and [40]).The importance of this loop entropy contribution would change with chain length N and may lead to an increase in T m ( N ) with increasing N [14] at least for a certain range of chain sizes . IX. SUMMARY AND CONCLUSION
This paper presents a novel theoretical model of DNA denaturation, already introduced in [28], which focuses onthe coupling between the base-pair link state (unbroken or broken) and the rotational degrees of freedom of thesemi-flexible chain. The Hamiltonian includes local chain bending rigidities whose values depend on neighboringbase-pair states: around 5 k B T for bubbles and 150 k B T for connected base-pair segments. Because of the rotationalsymmetry, the model can be rewritten in terms of an effective Ising Hamiltonian by integrating out the rotationaldegrees of freedom of the chain. Hence, our model yields considerable insight into the empirical temperature-dependentparameters used in previous Ising-like models [4]. In particular, the melting temperature T m is no longer a fittingparameter, but emerges naturally as a function of: (i) experimentally known bending rigidities, κ U and κ B ; (ii) thebare energy required to open a base-pair, 2˜ µ ; (iii) the bare energy of a domain wall, or destacking, 2 ˜ J ; (iv) thedifference in bare stacking energy between ss and dsDNA, 2 ˜ K ; and (v) the polymerization index, N . Moreover, ourmodel allows structural features of the DNA chain, such as the mean size R , to be calculated as a function of T . Anabrupt transition for R is found at T m and explains, at least qualitatively, the thermal transition observed in viscositymeasurements [42].From an experimental perspective, our results obtained from exactly solving the coupled model can be summarizedas follows. First of all, we propose formulæ for chain free-boundary conditions, this information being encoded inthe end vector | V i given in Eq. (11). However, any other boundary condition can be treated following the sameroute, even though we shall not detail the calculations here. For example, a polydA-polydT sequence of length N sandwiched between more stable G-C sequences [37] can be seen near its melting transition as a DNA of length N ,with fixed boundary conditions, resulting in an end vector | V i = | U i .Once boundary conditions are set, our model predicts melting profiles, as measured for example from UV absorbanceexperiments. Even if the model relies upon six microscopic parameters, as discussed in section VII, most of them areknown experimentally, including the strand length N , and only two of them must be extracted from melting profiles:˜ µ , the bare half-energy required to break a base pair (which can also be estimated experimentally [2]), and ˜ J , thecooperativity parameter that indicates the cost of creating a domain wall between unbroken and broken base pairs(recent experiments on single DNA molecules aim at determining this quantity, see [44]). Melting profiles, giving thefraction of broken base pairs, ϕ B , as a function of the temperature T , are determined through the average of the Isingstate variable, h c i , because ϕ B ( N, T ) = [1 − h c i ( N, T )] / infinite chains , ( N → ∞ ), the expression for h c i ∞ ( T ) is rather simple (Eq. (52)): h c i ∞ ( T ) = sinh( L )[sinh 2( L ) + e − J ] / , (123)where the renormalized parameters J , K , and L = µ + K are given in Eqs. (15), (16), and (21). In theselatter equations, the function G ( κ ) has a simple algebraic expression (Eq. (13)) that reduces to a purely entropiccontribution, G ( κ ) ≃ ln(2 κ ), in the physically relevant low- T (spin wave) approximation. From the melting profile ϕ B, ∞ ( T ), the melting temperature, T ∞ m , is defined by ϕ B, ∞ = 1 /
2, in other words by L = 0. The finite transitionwidth is estimated by Eq. (53): ∆ T ∞ m ≃ k B [ T ∞ m ] exp[ − J ( T ∞ m )] / ˜ µ .For finite length strands ( N finite), the expression for h c i ( N, T ), Eq. (64), is more complex, because N has aninfluence on ϕ B ( N, T ) and thus T m . In addition, the interplay with mechanisms, such as loop entropy, not taken intoaccount in this study is nontrivial, as discussed in detail in section VIII. At the level tackled in the present paper, weobtain a simplified expression for h c i ( N, T ) when N is large, Eq. (69), which simplifies even further when N ≫ ξ I ≫ h c i ( N, T ) ≃ h c i ∞ + 2 R V ( ξ I /N ) (1 − h c i ∞ ) / , (124)6where R V , Eq. (65), which is a ratio of matrix elements pertaining to end effects, simplifies at certain special temper-atures, see Eq. (118). Note that finite-size effects are still important for sizes of several thousands of base pairs (seesections VII and VIII and Figs 7 and 8) and are not a purely academic debate.Three important correlation lengths can be calculated in the framework of our model. On the one hand, the Isingcorrelation length ξ I , gives access to the typical size of bubbles in the low temperature regime ( T < T m ), as well asto the typical size of unbound regions for T > T m . This quantity is calculated in Eq. (66) and assumes a simplifiedform at T ∞ m : ξ I ( T ∞ m ) ≃ exp[2 J ( T ∞ m )] / ≫
1, when J ( T ∞ m ) ≫
1. On the other hand, one effective chain persistencelength, ξ p eff , CF ≃ [ ϕ U, ∞ /ξ pU + ϕ B, ∞ /ξ pB ] − , provides information on the short distant behavior of the chain tangent-tangent correlation function, Eq. (102), and the other, ξ p eff , provides information on the typical chain conformations,in particular its mean-square-radius: h R i ≃ a N ξ p eff ≃ a N ( ϕ U, ∞ ξ pU + ϕ B, ∞ ξ pB ) , (125)where a is the monomer length (0.34 nm) and the approximation is valid for very long chains (large N ). Knowing thisquantity is of primary importance when interpreting data from Tethered Particle or Tweezer experiments [41, 45, 46],Atomic Force Microscopy [47] or viscosity measurements [42]. The results that we have obtained here for the chaintangent-tangent correlation function and mean-square-radius are very different from what one obtains by solvinga quenched random rigidity model, where the local joint rigidity can take on one of two values, κ and κ , withprobability ϕ and ϕ = 1 − ϕ . In this case there is only one effective correlation length: h t i · t i + r i = e − r/ξ p ran and h R i ≃ a N ξ p ran , where ξ p ran ≡ − / ln[ ϕ e − /ξ p ( κ ) + ϕ e − /ξ p ( κ ) ] and ξ p ( κ ) is given by Eq. (32).The foregoing analysis makes allowance for neither solvent entropic (hydrophobic), nor electrostatic effects, whichmight not only change the actual value of the bare Ising parameters ˜ J , ˜ K and ˜ µ , but also lead to additional entropiccontributions. To go further concerning solvent entropic effects would require molecular dynamics simulations, whichis beyond the scope of the present approach. The electrostatic effects in DNA melting are two-fold: an entropiccontribution arising from the difference in electrostatic energy between states U and B and an an enthalpic contributionarising from counterion release. At equilibrium the two effects partially compensate and the remaining contributionis, using a simple thermodynamic approach in the low salt limit [48, 49], approximately equal to ∆ G el = k B T ℓ b ( τ B − τ U ) ln( I/I ), where I is the ionic strength ( I = 1 M), ℓ b = e / (4 πǫk B T ) the Bjerrum length ( e is the elementarycharge and ǫ the water dielectric permittivity) and τ U and τ B are the linear charge densities of the pure U and Bchains, respectively. Although this is an approximate result and the determination of τ B − τ U remains controversial,these two contributions should be included in a more refined model. A natural extension emerging from this studyconcerns the ionic strength dependence of DNA melting profiles and effective persistence lengths.The current rapid development of force experiments in magnetic or optical tweezer traps [45], “Tethered ParticleMotion” experiments [41] or even more recently atomic force microscopy ones [47], presents a formidable opportunityto investigate directly the elastic properties of DNA strands as a function of temperature, salt concentration andlength N . Indeed, despite the pioneering work by Blake and Delcourt [10], a systematic experimental study of theeffect of N while controlling the nature of chain ends, is lacking, especially for homopolymers. This would be a wayto discriminate between the different models and also shed light on the role of loop entropy, which is neglected in thepresent model.The approach developed here can be extended to bubble dynamics. Very recent work [39, 50, 51] studied thegrowth of already nucleated bubbles using the Fokker-Planck equation applied to the Poland-Scheraga model (i.e.,an effective Ising model including loop entropy). The agreement with experimental results obtained by fluorescencecorrelation spectroscopy is remarkably good [39]. The issue of bubble nucleation is, however, not solved and will becontinued to be explored in the near future. In addition the mutual influence of thermally excited bubbles and chainflexibility should play an important role in determining global DNA conformations and strongly influence loopingdynamics [46, 52]. [1] J. SantaLucia Jr., Proc. Nat. Acad. Sci. USA [2] P. Pincet, E. Perez, G. Bryant, L. Lebeau, and C. Mioskowski, Phys. Rev. Lett. , 3091 (2006).[4] R.M. Wartell and E.W. Montroll, Adv. Chem. Phys.
129 (1972).[5] O. Gotoh, Adv. Biophys.
127 (2004).[8] M.Th. Record, Jr., Ch.F. Anderson, and Lohman T.M., Q. Rev. Biophys.
103 (1978).[9] P. Nelson,
Biological Physics. Energy, Information, Life. , W.H. Freeman and Compagny, New York, 2004, Section 9.[10] R. D. Blake and S. G. Delcourt, Biopolymers. J.
935 (2007).[12] C.J. Benham, J. Mol. Biol. , 425 (1996); R.M. Fye and C.J. Benham, Phys. Rev. E , 3408 (1999).[13] R.M. Wartell and A.S. Benight, Physics Rep.
67 (1985).[14] D. Poland and H.R. Scheraga,
Theory of Helix Coil Transition in Biopolymers , Academic Press, New York, 1970.[15] B.H. Zimm and Bragg J.K., J. Chem. Phys.
526 (1959); B.H. Zimm, J. Chem. Phys. et al. , Bioinformatics , 370 (1999).[17] R. D. Blake and S. G. Delcourt, Nucleic Acids Res. , 4988 (2000).[19] E. Carlon et al. , Phys. Rev. Lett. , 198101 (2002).[20] M. Peyrard and A.R. Bishop., Phys. Rev. Lett. , 2755 (1989).[21] T. Dauxois, M. Peyrard, and A.R. Bishop., Phys. Rev. E , 684 (1993).[22] D. Cule and T. Hwa, Phys. Rev. Lett. , 2375 (1997).[23] M. Peyrard, Nonlinearity , R1 (2004).[24] Y. Gao, K.V. Devi-Prasad and E.W. Prohovsky, J. Chem. Phys. , 6291 (1984).[25] J.-H. Jeon, W. Sung and F.H. Ree, J. Chem. Phys. , 164905 (2006).[26] T.E. Cheatham and P.A. Kollman, J. Mol. Biol.
434 (1996).[27] T. Garel, C. Monthus, and H. Orland, Europhys. Lett.
132 (2001).[28] J. Palmeri, M. Manghi, and N. Destainville, Phys. Rev. Lett. , 088103 (2007).[29] J. Yan and J.F. Marko, Phys. Rev. Lett. Dynamical Phenomena at Interfaces, Surfaces and Membranes , Eds. D. Beysens, N. Boccaraand G. Forgacs, Nova Science Publishers, Inc., New York, 1993, 323.[31] C. Storm and P.C. Nelson, Europhys. Lett.
760 (2003).[32] M.E. Fisher, Am. J. Phys.
343 (1964).[33] G.S. Joyce, Phys. Rev.
478 (1967).[34] M. Grifoni and P. H¨anggi, Phys. Rep.
229 (1998).[35] D. Chandler,
Introduction to Modern Statistical Mechanics , Oxford University Press, New York, 1987, Section 5.8.[36] C. Cohen-Tannoudji, B. Diu and F. Lalo¨e,
M´ecanique Quantique , Hermann, Paris, 1973, vol. 2.[37] G. Altan-Bonnet, A. Libchaber, and O. Krichevsky, Phys. Rev. Lett. , 3212 (2003).[39] H. C. Fogedby and R. Metzler, Phys. Rev. Lett. , 070601 (2007).[40] R. Blossey and E. Carlon, Phys. Rev. E , 061911 (2003)[41] N. Pouget et al. , Nucleic Acids Res. , e73 (2004).[42] R.B. Inman and R.L. Baldwin, J. Mol. Biol. , 452 (1964).[43] T. Garel and H. Orland, Biopolymers , 453 (2004).[44] C. Ke, M. Humeniuk, H. S. Gracz, and P. E. Marszalek, Phys. Rev. Lett. , 018302 (2007).[45] S.B. Smith, L. Finzi, and C. Bustamante, Science et al. , Nucleic Acids Res. , 4313 (2006).[47] P.A. Wiggins et al. , Nature Nanotechnology
137 (2006).[48] G.S. Manning, Biopolymers
937 (1972).[49] N. Korolev, A.P. Lyubartsev, and L. Nordenski¨old, Biophys. J. ,128105 (2006).[51] A. Bar, Y. Kafri and D. Mukamel, Phys. Rev. Lett. , 038103 (2007).[52] L. Finzi and J. Gelles, Science , 378 (1995). Symbol Quantity Mathematical definition Reference T m melting temperature ϕ U ( T m ) = ϕ B ( T m ) = 1 / N chain length Sec. IIU / B unbroken/broken base pair Sec. II ϕ U,B fractions of U’s and B’s Eq. (6) σ i internal degree of freedom (Ising variable) ± t i chain unit tangent vector (Heisenberg variable) k t i k = 1 Sec. II˜ κ i,i +1 local chain bending rigidity (coupling) κ U , κ B or κ UB Sec. II˜ J half-energy of a domain wall Sec. II˜ K difference in stacking energy between ds and ssDNA Sec. II˜ µ half-energy required to open a base-pair Sec. II κ, J, K, µ, . . . adimensional energies (in units of k B T ) κ = β ˜ κ , . . . Sec. II G ( κ ) free energy of a single joint of rigidity κ κ − ln[sinh( κ ) /κ ] Sec. II J renormalized (effective) Ising parameter J J − [ G ( κ U ) + G ( κ B ) − G ( κ UB )] Eq. (15) K renormalized (effective) Ising parameters K K − [ G ( κ U ) − G ( κ B )] Eq. (16) L effective chemical potential to open an interior base pair µ + K Sec. II h c i ∞ infinite size Ising “magnetization” sinh( L ) / [sinh ( L ) + e − J ] / Eq. (52) ϕ U, ∞ , ϕ B, ∞ fraction of unbroken (U) and broken (B) bonds ( N → ∞ ) [1 ± h c i ∞ ] / ϕ B,i site i bonding opening probability (melting map) [1 − h σ i i ] / T ∞ m infinite size ( N → ∞ ) melting temperature L ( T ∞ m ) = 0 Sec. II∆ T ∞ m transition width ( N → ∞ ) 2 | ∂ h c i ∞ /∂T | − at T ∞ m & Figs. 5 and 6 Eq. (53) h σ i + r σ i i Ising correlation function Sec. IV ξ I Ising correlation length e J / T ∞ m & Fig. 4 Sec. IV h t i + r · t i i chain correlation function Fig. 3(b) Sec. V ξ p , + , ξ p , − persistence lengths Fig. 4 Eq. (100) ξ p eff effective persistence length Fig. 4 Sec. VI R chain mean square (or gyration) radius h R i / Secs. II,VI T ∗ crossover temperature for finite chains ∂ϕ B /∂N | T ∗∗