Table Of ContentStructure factor and dynamics of the helix-coil transition
Herv´e Kunz1, Roberto Livi2, and Andra´s Su¨t˝o3
1 Institute of Theoretical Physics,
Ecole Polytechnique F´ed´erale de Lausanne,
7 CH-1015 Lausanne, Switzerland
0 2 Dipartimento di Fisica, Universit´a di Firenze,
0 INFM and INFN, Firenze,
2 via Sansone 1 I-50019 Sesto Fiorentino, Italy
n 3 Research Institute for Solid State Physics and Optics,
a Hungarian Academy of Sciences,
J P. O. Box 49, H-1525 Budapest 114, Hungary
7
1 (Dated: February 6, 2008)
Thermodynamical properties of the helix-coil transition were successfully described in the past
]
h by the model of Lifson, Poland and Sheraga. Here we compute the corresponding structure factor
c and show that it possesses a universal scaling behavior near the transition point, even when the
e transition is of first order. Moreover, we introduce a dynamical version of this model, that we
m
solve numerically. A Langevin equation is also proposed to describe the dynamics of the density
- of hydrogen bonds. Analytical solution of this equation shows dynamical scaling near the critical
t temperatureand predicts a gelation phenomenon above thecritical temperature. In the case when
a
t comparison of the two dynamical approaches is possible, the predictions of our phenomenological
s theory agree with theresults of theMonte Carlo simulations.
.
t
a PACS:87.10.Gg, 64.60.Fr, 64.60.Ht
m
-
d
I. INTRODUCTION thestructurefactorandfindthatnearthecriticalpointit
n
o exhibits universal scaling behavior, which depends only
c Ever since the discovery of the helicoidal structure of on α. Moreover, we also look at the transition from a
[ dynamical point of view. We do it in two ways: first we
DNA under the usual conditions, the phenomenon of its
describe this dynamics by a master equation such that
1 transformation to a coil under heating has been studied
v and described as an equilibrium phase transition. Basic the equilibrium state is chosen to be the one of the LPS
5 worksbyLifsonandsubsequentlybyPolandandSheraga model. Secondly,weintroduceaLangevintype equation
0 for the orderparameterwhichis the density of hydrogen
led to a model (LPS) considered today standard [1, 2].
4 bonds. The analysis of the solution of this equation re-
The model is one-dimensional and shows a phase tran-
1
vealsa criticalslowing-down,the relaxationtime diverg-
0 sition only because the interactions are long-ranged. A
ing at the criticalpoint with a power law still depending
7 more realistic description of the two strands moving and
onα. Intheabsenceofnoise,wefindthatabovethecrit-
0 interactinginthree dimensionalspacehasalsobeen pro-
/ posed and studied numerically [4, 5, 6? ] . A model of ical temperature a gelation phenomenon occurs, namely
t
a LPS–typewasshowntobeabletoconsistentlyfitallrel- the density of H-bonds approaches zero in a finite time.
m With the addition of noise we obtain also the analytic
evantsimulationresults[7]. Morerecently,someauthors
expression for the power–spectrum of the density of H-
- have proposed a new mechanism which yields an abrupt
d bonds. Finally,wecomparethetwodescriptionsandfind
first–order phase transition [8, 9], at variance with the
n that the dynamics of the order parameter derived from
previous theoretical considerations in favor of a continu-
o
the master equation and from the Langevin equations
c ous one. On the experimental side the situation is still
essentially agree.
: far from conclusive. Previous measurements indicated
v
i the presence of a sharp jump in the fraction of bound The following three sections are devoted to a study of
X basepairs[10, 11]. Morerecently,the resultsreportedin equilibrium properties. In Section II we summarize the
r [12] have been interpreted as an indication of a weakly basics of the LPS model. Formulas for the free energy
a
continuous phase transition. andthe pressurearegivenhere,aswellasscalingformu-
It seems, however, that little is known about the dy- las for the density inthe vicinity of the phase transition.
namics of this denaturation transition. We can mention SectionIIIcontainsthecomputationoftheone-andtwo-
thecontributionsbasedonamechanicalapproach[13,14] point correlation functions for finite systems and in the
and on a Langevin description of the two DNA–strands thermodynamic limit. Our results for the structure fac-
as polymers in continuous space [15]. torarepresentedinSectionIV.Results onthe dynamics
In this paper we investigate the original LPS model are collected in Sections V and VI. In Section V we de-
while considering the exponent α as a free parameter. scribe the microscopic dynamics realized numerically in
Forwhatconcernstheequilibriumpropertieswecompute this work. The equation for the order parameter and its
2
analysis both with and without noise is given in Section and
VI,wherewealsocomparethe resultsofthetwodynam-
βf(ρ)=sup{ǫρ−βp(ǫ)}≡ǫ ρ−βp(ǫ ) . (2.7)
ical approaches. The paper ends with a Summary. ρ ρ
ǫ
Here ρ is the density of 1 (H-bonds), and ρ and ǫ are
ǫ ρ
the respective values at which the suprema are attained.
II. THE LPS MODEL
The partition functions depend on the boundary con-
ditions imposed on the system. Typical choices are the
Let 0 and 1 denote respectively a broken or existing
free and the periodic conditions, or fixed 0 or 1 on the
H-bond at a given site. Then, a configuration of the
boundary sites. For example, with the choice 0...1 or
double-strand of DNA of length N is represented by a
1...0,
sequence n={n }N where n =0 or 1. In n there are
x x=1 x
alternating maximal connected subsequences (intervals)
Q = A (n)B (N −n) (2.8)
N,n k k
of 1 and 0, respectively. A weight a(j)>0 or b(j)>0 is
k≥1
associatedtoanintervaloflengthjof1or0,respectively. X
Ifthe lengthsofthe intervalsofnarel andm thenthe where
i j
weight of n is
k
A (n)= ··· a(l )δ (2.9)
e−V(n) = a(l ) b(m ) . (2.1) k i n,Pli
i j lX1≥1 lXk≥1iY=1
i j
Y Y and
As usual, thermal equilibrium can be discussed in the
k
canonical or in the grand-canonical ensemble. In the B (n)= ··· b(m )δ . (2.10)
k i n,Pmi
first case the description is restricted to configurations
mX1≥1 mXk≥1iY=1
with the total number of H-bonds, n= n , fixed. In
x x
the second case all configurations compatible with the Ak(n)=Bk(n)=0ifn≤0ork >n,thereforeQN,n =0
boundary condition are considered butPan activity eβµ for n ≥ N. One can show that in the thermodynamic
is introduced. The chemical potential µ is the energy limit the free energy and the pressure do not depend on
neededtobreakaH-bond. Thecanonicalpartitionfunc- the boundary condition. Below we give these quantities
tion is in the case of the simplest model (the LPS model in a
strict sense) obtained by choosing
Q = e−V(n) (2.2)
N,n a(l)=a b(l)=bull−α (2.11)
n:PXxnx=n
where u,α>1 and a,b>0 . With the notations c=ab,
and the grand-canonicalpartition function is
∞
Q (βµ)= Q enβµ = e−U(n) (2.3) ζα(x)= l−αe−lx (x≥0) (2.12)
N N,n
l=1
n n X
X X
and
where U(n) = V(n)−βµ n . Alternately, e−U(n) is
obtained from e−V(n) by repxlaxcing each factor a(l) by ǫc =lnu−ln[1+cζα(0)], (2.13)
P
a(l)elβµ. The free energy f and the pressure p are re- one obtains
spectively given by
βp=lnu+Θ(ǫ−ǫ )x(ǫ) (2.14)
c
1
βf =− Nl→im∞ N lnQN,n (2.4) where Θ is the Heaviside function and, for ǫ>ǫc, x(ǫ) is
the unique solution of the equation
n/N →ρ
uex−ǫ =1+cζ (x) . (2.15)
α
and
Thus, x(ǫ ) = 0 and x(ǫ) is monotone increasing with
c
1 ǫ. For ǫ < ǫ , βp takes on its minimum value, lnu.
βp= lim lnQ (βµ). (2.5) c
N→∞N N Comparison with (2.4), (2.5) and (2.11) shows that the
corresponding density and free energy are ρ ≡ 0 and
ǫ
Thetemperaturedoesnotappearexplicitlyinthecanon- f(ρ ) ≡ −lnu. Legendre transformation yields f(ρ) for
ǫ
ical partition function, therefore βf is independent of β all ρ. With ρ ≡ρ ,
c ǫc
and βp depends on it through ǫ = βµ. The free energy
and the pressure are Legendre transforms of each other, −lnu+ǫcρ, ρ≤ρc
βf(ρ)=
−lnu−y+ρ y+ln u , ρ≥ρ
βp(ǫ)=sup{ǫρ−βf(ρ)}≡ǫρ −βf(ρ ) (2.6) ( 1+cζα(y) c
ǫ ǫ (cid:16) (cid:17) (2.16)
ρ
3
where y =y(ρ) is the solution of the equation x and y we have 0’s. It follows that
cζ (y) N−1 r−1
1+ α−1 =ρ−1. (2.17) Q11(N)ρ (r)=a Q10[1,x]eǫ(y−x−1)Q01[y,N]
1+cζ (y) N
α
y=r+1x=2
X X
(3.22)
In fact, y(ρ) = x(ǫ ); y(ρ) is defined for ρ ≤ ρ ≤ 1 and
ρ c
By introducing
is monotone increasing, y(ρ )=0 and y(1)=∞. There-
c
fore, f(ρ) is monotone increasing from f(0)=−β−1lnu
N
to f(1) = 0; as the Legendre transform of the convex g(N)= Q01(x)e−ǫx (3.23)
p(ǫ), f(ρ) is also convex. Substituting y = 0 into (2.17)
x=2
yields X
we can rewrite (3.22) as follows:
1+cζ (0)
α
ρc = 1+c[ζ (0)+ζ (0)]. (2.18) aeǫN[g(r−1)+1][g(N −r)+1]
α−1 α ρ (r)= . (3.24)
N Q11(N)
Thisshowsthatρ =0ifα≤2andρ >0ifα>2. The
c c
first case corresponds to a continuous phase transition For the two-point correlation function ρN(r,r′) we ob-
with the order parameter ρ vanishing continuously as ǫ serve that in a configuration either r and r′ are in the
ǫ
approachesǫ from above. In the secondcase the transi- same cluster of 1’s or they belong to distinct clusters of
c
tion is of first order: ρ decreases continuously to ρ as ǫ 1’s. By assuming r′ ≥r+1 we obtain
ǫ c
decreasestoǫ andjumpstozerobelowǫ . Computation
αof≤ρ′ǫ3=(wdhρiǫle/idctǫisfrfionmiteEfqo.r(α2.>173))s.hTowhesttyhcpaetoρf′ǫcd+iv0e=rg∞encief ρN(r,r′)= aeǫN[g(r−1Q)+111(N][g)(N −r′)+1] ×
can be extracted from the asymptotic forms r′−r
1+a Q00(l)e−ǫl(r′−r−l) (3.25)
ρ ∼(ǫ−ǫ )2−α if 1<α<2
ǫ c l=1
X
ρ −ρ ∼ (ǫ−ǫc)α−2 if 2<α<3 (2.19)
ǫ c |ln(ǫ−ǫ )|−1 if α=2 These expressions can be simplified by considering that
c
(cid:26) the following relations hold for N ≥2:
validwhenǫ decreasesto ǫ . Equation(2.19)is obtained
c Q01(N) = uN+1e−ǫZ −uNZ
from the asymptotic form of ζ (x) for small x and from N N−1
α
x(ǫ) ≈ ǫ−ǫ for ǫ−ǫ ≪ 1. It will be used in Section uN
c c Q00(N) = u2e−2ǫZ −2ue−ǫZ +Z
VI. Finally, we note that for ǫ < ǫc (T > Tc), when a N+1 N N−1
ρǫ = 0, the decay of the density with increasing N can Q11(N) = auN(cid:2)ZN−1 (3(cid:3).26)
be computed,
where
2
ρN ∼ (1−e−ǫc(1−Tc/T))N . (2.20) Z = 1 dw 1 , (3.27)
N 2πi wN+1 ue−ǫ−R(w)
IC
This formula is valid for N ≫(1−T /T)−1.
c C is some circle around the origin and
∞
III. CORRELATION FUNCTIONS R(w)=w+c l−αwl+1 (3.28)
l=1
X
We will compute the one- and two-point correlation
Upon these results we have
functionsofH-bondsinthegrand-canonicalensemblefor
the LPS model. Let us denote the partition function of g(N)+1=uN+1e−ǫ(N+1)Z (3.29)
N
the box [x,y] by
and finally
Qαβ[x,y]=Qαβ(y−x+1) (3.21)
Z Z
ρ (r) = ue−ǫ r−1 N−r (3.30)
N
where we have imposed the boundary conditions α at x ZN−1
Qan1d1(β1)a=t aye(ǫN).oWteitthhoauttQpαreβju=dicQeβoαf,gQen0e0r(a1l)ity=wbeucoamnd- ρN(r,r′) = u2e−2ǫZr−1ZZN−r′Zr′−r (3.31)
N−1
pute the correlation functions with the boundary condi-
tionsα=1andβ =1. Webeginwiththedensityρ (r). where r′ ≥r.
N
In a configuration the site r is in a cluster of 1’s begin- Inthe low–temperaturephase(ǫ>ǫ )wecanperform
c
ningatasitex+1,endingatasitey−1,sothatatsites the thermodynamic limit(N →∞) byfixing |r−r′|and
4
letting the distance of r to the boundaries also tend to In this expression we have reparametrized w∗ as
infinity. Making use of the asymptotic (N ≫1) relation
w∗ =e−δ. (4.41)
Z ∼λNq(1+e−κN) (3.32)
N
Accordingly, T = T corresponds to δ = 0. The ex-
c
with κ>0, pression(4.40)depends onthe details of the LPSmodel.
However,wemayexpectthatinthecriticalregion,char-
λ−1 = euǫ w∗ (3.33) aemcteerrgizeeadnbdyt|hke|<re<su1ltsanddep|eTn−TdcTco|n<ly<on1tuhneivpearrsaamlfeetaetrurαes.
If this is the case, we can be confident of the predictions
and
of the model. Let us focus on the analysis of the critical
dR(w) region. The parameterδ allowsusto defineacorrelation
q−1 =w∗ dw |w=w∗ , (3.34) length ξ through the relation
where w∗ is the positive solution of the equation ue−ǫ = ξ =δ−1. (4.42)
R(w∗), we obtain
Near T , we find
c
ρ(r)=ρ=ue−ǫq =∂βp/∂ǫ, (3.35)
ξ = ξ0(α)(Tc/T −1)−α−11 if 1<α<2 (4.43)
the expected density of H-bonds, and ξ (α)(T /T −1)−1 if 2<α<3.
(cid:26) 1 c
ρ(r,r′)=qu2e−2ǫλ−|r′−r|Z|r′−r|. (3.36) In order to obtain explicit expressions we need to study
the function ζ (δ), cf. Eq. (2.12), for small values of δ.
α
We can see that the truncated correlation function For this purpose we can use the integral representation
ρT(r,r′)=ρ(r,r′)−ρ(r)ρ(r′) decays exponentially
1 +∞ tα−1
ρT(r,r′)∼ρ2e−κ|r′−r| (3.37) ζα(δ)= Γ(α) et+δ−1dt (4.44)
Z0
for |r′ −r| → ∞. We note that translation invariance valid for Reδ >0. By this relation we can write
is lost in the high–temperature phase, because there re-
mainsanexplicitdependenceonr andr′ inρ(r,r′),even ζα(δ)−ζα(δ−ik) (4.45)
if r and r′ arefar fromthe boundaries. This can be con- =(e−ik−1) ζ (δ)+(1−e−ik)B (k,δ)
α−1 α
sidered as a weakness of the LPS model. Nonetheless,
we shall see that we can still attribute a clear meaning where (cid:2) (cid:3)
to the structure factor in this phase. 1 +∞ tα−1
B (k,δ)= dt.
α Γ(α) (et+δ−ik−1)(1−e−t−δ)2
Z0
(4.46)
IV. STRUCTURE FACTORS
For 1 < α < 3 and δ → 0 the above expression reduces
to
The structure factor is more accessible to experiments
thanthecorrespondingcorrelationfunction. Itisdefined B (k,δ)≈δα−3b (k) (4.47)
α α
as follows
with
S (k)= 1 ρT(r,r′)eik(r′−r). (4.38)
N N N Γ(2−α) (1−ik)α−1−1 Γ(2−α)
Xr,r′ bα(k)= α−1 k2 +i k .
(cid:20) (cid:21)
We aim at computing the limit S(k) of this quantity as (4.48)
N → ∞. The two cases T < T and T > T are treated In particular,
c c
separately.
Γ(3−α)
Below T — The correlation functions become transla- b (0)= . (4.49)
c α
tionalinvariantforN →∞,i.e. ρT(r,r′)→ρT(0,r′−r). 2
N
Moreover ρT(0,r) decays exponentially and we obtain Thecompressibilityχ=S(0)cannowbecomputed. For
1<α<2thephasetransitioniscontinuousandwehave
S(k)= ρT(0,r)eikr. (4.39)
ρ ≈ ρ ξα−2 (4.50)
r 0
X
χ ≈ ρ2(2−α)ξ2α−3. (4.51)
Making use of (3.27), (3.28) and (3.37) we have 0
For 2 < α < 3 the phase transition is of first order and,
S(k) eδ
=2Re −1. as ρ→ρc, one has
ρ 1+cζ (δ)−eik[1+cζ (δ−ik)]
(cid:20) α α (cid:21)
(4.40) χ≈χ ξ3−α. (4.52)
0
5
We can also find the scaling limit for S(k). It is defined is atthe originandthe singularityis integrable. If, close
by taking the limits k →0 and ξ →∞ with the product to the origin,
kξ fixed at some finite value. If 1<α<2 we find
F(x)∼ p (s)|x|rj, s= sgn(x) (4.63)
j
S(k) 2(α−1) 1
= Re . (4.53) j
S(0) 2−α (1−ikξ)α−1−1 X
(cid:20) (cid:21)
then
For small kξ this gives a Lorentzian behavior
a
S(k) 4 F(x)e−iNxdx∼ pj(s)Γ(1+rj)(Nei2πs)−(1+rj) ,
= (4.54)
S(0) (kξ)2+4 Z−a sj
X
(4.64)
while in the critical region one has where rj is a sequence of real numbers such that −1 <
r < r < ··· . Using Eqs. (4.44)–(4.48), this formula
1 2
S(k) 2(α−1)cosπ(α−1) yields
= 2 . (4.55)
S(0) (2−α)(kξ)α−1
Z ∼c/(δ¯2Nα) (4.65)
N
On the other hand, for 2<α<3 we have
where δ¯=v−v . Moreover,if k 6=0 then
c
S(k) 2Re 1−(1−ikξ)α−1
= . (4.56)
S(0) (α−1)(α−2)(kξ)2 C (k)∼g(k)+g2(−k)eik(N−1) (4.66)
(cid:2) (cid:3) N
For |kξ| << 1 one has again a Lorentzian, but in the
and
critical region (4.56) takes the form
D (k)∼g(k)+g(−k)eik(N−1). (4.67)
S(k) 2 cosπ(α+1) N
= 2 . (4.57)
S(0) (α−1)(α−2)(kξ)3−α
This shows that the compressibility
Above T — In this phase the properties of the
c χ=S(0)∼2vv /(Nδ¯3). (4.68)
LPS model are more dependent on N than in the low- c
temperature one. Basically, one has to provide a suit-
Accordingly, in the limit N →∞ we have
able description of the partition function Z , defined in
N
(3.27). Let us introduce the parameter v =ue−ǫ. It can
S(k) v δ¯
easily be shown that if v >vc then = |g(k)|2 (2/δ¯)Reg−1(k)−1 − . (4.69)
S(0) v v
c c
1 +π (cid:2) (cid:3)
ZN = 2π e−iNθg(θ)dθ (4.58) Sinceδ¯=δ¯0(T/Tc−1),alsointhiscasewecanobtainthe
Z−π scalinglimit. Withthedefinition(4.43)ofthecorrelation
where length ξ, for 1 < α < 2 the normalized structure factor
reads
g−1(θ)=v−eiθ[1+cζ (−iθ)]. (4.59)
α
S(k) 1+2(kξ)α−1sinπα
= 2 (4.70)
In this case the structure factor can be expressed in the S(0) 1+(kξ)2(α−1)+2(kξ)α−1sinπα
form 2
while, for 2<α<3,
NS (k)=2v2ReC (k)−v2|D (k)|2−vD (0) (4.60)
N N N N
S(k) 1
where = . (4.71)
S(0) 1+(kξ)2
+π dθ
C (k)=Z−1 e−i(N−1)θg2(θ)g(θ+k) (4.61)
N N−1 2π We stress that the main result of these computations is
Z−π the existence of a universal scaling limit for the normal-
and ized structure factor S(k)/S(0), which depends only on
the parameterα. This suggestsalso that the LPSmodel
+π dθ belongs to a new universality class, in the language of
D (k)=Z−1 e−i(N−1)θg(θ)g(θ+k). (4.62)
N N−1 2π critical phenomena. Accordingly, its predictions about
Z−π
theprocessofDNAdenaturationcanbetakenwithsome
In order to find the asymptotic behavior (N → ∞) of confidence. Moreover, we expect that this result may
theseexpressionsweusethefollowingformulae. LetF(x) help to clarify experimentally the nature of the phase
beafunctionwhoseonlysingularityintheinterval[−a,a] transition, either continuous or first–order.
6
V. MICROSCOPIC DYNAMICS If γ is an interval (a sequence of consecutive integers
between 1 and N) then by definition the distance of x
We will describe the evolution of the configuration n to γ is d(x,γ) = min{|x−y| : y ∈ γ}. Because U(n) is
by adiscrete-time localdynamics. Althoughthe method additive in the contributions of the intervals of n,
is well-known, we present it in some detail.
In a time step, n can change only at a single site x. U(n)=− [lna(li)+liǫ]− lnb(mj) , (5.78)
The new configuration is Txn and, hence, (Txn)y = ny Xi Xj
if y 6= x, and (T n) = 1− n . Note that T2 equals
x x x x
the identity. The transition probability from n to m is z(n,x)dependsonlyontheintervalsofnandTxnwhose
pn,m. More precisely, it is the conditional probability distance to x is 0 or 1. If in n there is a single interval
that, if the configuration is n in time t, it will be m in with this property, there will be three in Txn, and vice
time t+1. It then follows that pn,m = 0 unless m = n versa; if there are two in n, there will also be two in
or m=Txn for some x. The irreducibility of the matrix Txn. We obtain altogether three different functions of
pn,m and, hence, the ergodicity of the dynamics and the the lengths of intervals, and their reciprocals. These are
uniquenessoftheequilibriumstateisguaranteedifpn,Txn listed in Table 1.
is indeed positive for all x. (Note, however, that both
n and Txn have to satisfy the boundary conditions. If (1)l11(1)l2 → (1)l10(1)l2 R1
thesearegivenbyfixingn1 andnN then2≤x≤N−1.) (0)l11(1)l2 → (0)l10(1)l2 R2
In a numerical experiment the transition probability is (1)l11(0)l2 → (1)l10(0)l2 R3
decomposed as (0)l11(0)l2 → (0)l10(0)l2 R4
(1)l10(1)l2 → (1)l11(1)l2 R−1
pn,Txn =pxW(Txn|n,x) (5.72) (0)l10(1)l2 → (0)l11(1)l2 R1−1
2
where px is the probability to choose x for an attempt (1)l10(0)l2 → (1)l11(1)l2 R3−1
of change, and W(Txn|n,x) is the conditional probabil- (0)l10(0)l2 → (0)l11(0)l2 R−1
4
ity that once n and x have been selected, the change is
executed. Clearly, xpx = 1 must hold; for periodic or TABLEI: Transition ratesz(n,x). Thenotations(X)l1 and
freeboundaryconditionsthetypicalchoiceforpx is1/N. (X)l2 indicate a left-cluster of length l1 and a right-cluster
The probability toPstay in n if the change is attempted of length l2 made of X-symbols, respectively. R1 = abv,
in x is R2 = v(l1l+11)α, R3 = v(l2l+21)α, R4 = avb(l1+l1ll22+1)α where
v=ue−ǫ.
W(n|n,x)=1−W(T n|n,x) (5.73)
x
Typical choices of W are
and, thus,
N N z(n,x)
W(T n|n,x)=min{1,z(n,x)} or . (5.79)
pn,n =1− pn,Txn = pxW(n|n,x) . (5.74) x 1+z(n,x)
x=1 x=1
X X
They satisfy (5.77) because z(T n,x) = 1/z(n,x). We
x
This is nonzero unless W(n|n,x)=0 for all x.
shall work with the first one.
Let p (n), t = 0,1,2,..., be the probability of n at
t
time t. The master equation of the evolution is
VI. THE ORDER PARAMETER EQUATION
N
pt+1(n)=pt(n)pn,n+ pt(Txn)pTxn,n . (5.75)
x=1 The microscopic dynamics of denaturation as de-
X
scribed by the master equation is quite complex and un-
The local transition probability W is to be chosen so
tractable analytically. Most of the time, however, one
that the grand-canonical distribution π be the unique
is interested in the evolution of the order parameter,
equilibrium solutionof(5.75). Imposing the conditionof
namely the density ρ of 1 (H-bonds). This quantity, ob-
detailed balance
tained as an average over the sample, should vary much
more slowly than the other degrees of freedom to which
π(n)W(T n|n,x)=π(T n)W(n|T n,x) , (5.76)
x x x
itiscoupled. Thiscouplinghastwoeffects. (i)Itcreates
we still have an infinity of choices: any W satisfying an effective thermodynamic force exerted on the order
parameter and depending nonlinearly on it. This force
W(T n|n,x) shouldvanishifthedensity,andviaitthepressure,reach
x =e−U(Txn)+U(n) =z(n,x) (5.77)
W(n|T n,x) their equilibrium value, corresponding to the prescribed
x
value of ǫ. We choose it therefore in the form
and 0<W(T n|n,x)≤1 for all x and all n of length N
x
with the right boundary conditions will do. −[(βf)′(ρ)−βµ]
7
where f(ρ) is the equilibrium free energy for each value 2. ρ > ρ . Then for 0 ≤ t ≤ t , ρ evolves according to
0 c c
of ρ, cf. Eqs.(2.16). So the deterministic part ofthe dy- ρ˙ = −γ(ǫ −ǫ), where t is defined by ρ(t ) = ρ . For
ρ c c c
namics realizes the searchof the supremum in Eq. (2.6). t>t the force is constant as in 1., hence
c
One can recognize here also the mean-field approach to
the problem of the critical slowingdown in phase transi- ρ(t)=ρc−γǫc(1−Tc/T)(t−tc) for tc ≤t≤tf (6.86)
tions. (ii) Due to the very large number of the coupled
degrees of freedom, the coupling creates a random force where tf =tc+ρc/γǫc(1−Tc/T) and ρ(tf)=0.
of zero average, responsible for the fluctuations of the What we see here is the phenomenon of gelation: The
density. We choose the random force as a white noise system reaches equilibrium in a finite time [16].
η(t) with a correlation (ii) T < T (ǫ > ǫ ). Then lim ρ(t) = ρ > ρ .
c c t→∞ ǫ c
For large t, by linearizing ǫ about ǫ we find ǫ −ǫ ≈
hη(t)η(t′)i=2γβ−1δ(t−t′). (6.80) (ρ−ρ )/ρ′ [recall that ρ′ =ρdρ /dǫ] and ρ
ǫ ǫ ǫ ǫ
Thus, the time evolution of ρ(t) is governed by the
Langevin equation ρ(t)−ρǫ ∼Ae−[γ/ρ′ǫ]t. (6.87)
ρ˙ =−γ[(βf)′(ρ)−βµ]+ β η(t) (6.81) When T is close to Tc, from Eq. (2.19) we deduce ρ′ǫ ∼
(α−2)ǫ (T /T −1)α−3 and
c c
where γ > 0. We expect this equpation to describe cor-
rectlythebehavioroftheorderparameterinthe vicinity ρ(t)−ρǫ ∼exp(−t/τ) (6.88)
ofthe transitionpoint. The comparisonoftheanalytical
predictions based on Eq. (6.81) with the numerical solu- with
tion of the master equation can be a first check of this
τ =ρ′/γ ∼γ−1(α−2)(T /T −1)α−3.
hypothesis. From now on we consider only the LPS case ǫ c
(2.11). From Eq. (2.16) we deduce
CASE 2: 1<α<2, ρ =0.
c
(βf)′(ρ)=ǫρ, (6.82) (i) T > Tc. Now starting from ρ0 = ρ(0) > 0, for all
times ǫ ≥ ǫ > ǫ, therefore ρ˙ ≤ −γǫ (1 −T /T) and
ρ c c c
therefore (6.81) reads
ρ(t) attains zero in a finite time t < ρ /γǫ (1−T /T).
f 0 c c
Again, we find gelation.
ρ˙(t)=−γ[ǫ −ǫ]+ β η(t) (6.83)
ρ(t)
(ii) If T < T , for large times we find the same result as
c
As we have seen in the former secption, the microscopic in (6.87) and (6.88) with
dynamics depends only on α, c = ab and v. It is im-
portant to realize that ǫρ−ǫ also depends only on these τ =ρ′(ǫ)/γ ∼γ−1(2−α)(Tc/T −1)1−α.
parameters. Indeed, for given ρ we solve Eq. (2.17) for
y, replace x by y(ρ) in Eq. (2.15) and extract ǫ . Thus, CASE 3: α = 2. For T > T we obtain gelation, for
ρ c
T <T an exponential relaxationto ρ =ρ with decay
c ∞ ǫ
ǫρ−ǫ=y(ρ)−ln{1+cζα[y(ρ)]}+lnv. (6.84) time
On the other hand, we expect to obtain the value of γ τ ∼[γ(T /T −1)ln2(T /T −1)]−1.
c c
from the microscopic dynamics.
Note that ǫ ≡ ǫ if ρ ≤ ρ and tends to ∞ as ρ
ρ c c
As in the case of the static structure factor, we can
approaches1. Invertingρ showsthatǫ startswithzero
ǫ ρ
obtain a scaling limit for the time-dependent order pa-
derivative if α ≤ 3, while ǫ′ > 0 if α > 3. We will
ρc+0 rameter ρ(t). It can be determined by taking ρ → 0,
confine ourselves to the case α≤3.
t → ∞, T → T . We restrict this analysis to the case
c
α>3/2.
(i) T >T . The gelation time t diverges as
c f
A. Evolution without noise
2−α
t ∼(T/T −1)−(1−ν) where ν = .
CASE 1: 2 < α < 3. In this case ρ > 0. The evo- f c α−1
c
lution is naturally different above and below the critical
temperature. (ii) T <Tc. Then
(i) T > T (ǫ < ǫ ). Let ρ = ρ(0). We distinguish
c c 0
ρ(t) γǫ
between two subcases.
=z (1−T/T )t (6.89)
c
ρ ρ
1. ρ0 <ρc. Then ρ˙ =−γǫc(1−Tc/T) and ǫ (cid:18) ǫ (cid:19)
ρ(t)=ρ −γǫ (1−T /T)t for 0≤t≤t (6.85) where z(t) is the unique solution of
0 c c f
where t =ρ /γǫ (1−T /T) and ρ(t )=0. z˙ =−(z1/ν −1) (6.90)
f 0 c c f
8
tending to 1 as t tends to infinity. Therefore, if t(1− Making use of (6.95) and (6.96) one obtains
T/T )1−ν ≫1 then
c
C (t)=χe−|t|/τ (6.99)
ρ
γǫ
ρ(t)∼ρ 1+exp − (1−T/T )t . (6.91)
ǫ νρ c and
(cid:20) (cid:18) ǫ (cid:19)(cid:21)
On the other hand, in the critical region χτ 1
P (ω)= (6.100)
ρ π 1+(ωτ)2
1≪t≪|T/T −1|−(1−ν)
c
where τ =χ/γ. For T րT we obtain
c
one has
−ν/(1−ν) γ−1ξ2α−3 if 1<α<2
νγǫt τ ∼ (6.101)
ρ(t)∼ 1−ν . (6.92) ( γ−1ξ3−α if 2<α<3
(cid:20) (cid:21)
with ξ defined in Section IV. In the language of critical
B. The effect of noise phenomena Eq. (6.101) can be interpreted as a scaling
law valid for ωτ not too large. The dynamical critical
exponent z is then given by
In order to compute the time–dependent correlation
functions of the density ρ, we have added a white noise
term η(t) to the deterministic order parameter equation z = 2α−3 if 1<α<2 (6.102)
(6.81). On the other hand, this choice may lead to a ( 3−α if 2<α<3 .
consistency problem. In fact, it is well known that in
this case the stationary distribution P(ρ) of the density Near to but above Tc, one has instead
is of the form
τ ∼γ−1ξ2/N (6.103)
1
P(ρ)∼exp − [βf(ρ)−ǫρ] . (6.93)
2 since χ=2vvc/Nδ¯3, see (4.68) .
(cid:26) (cid:27)
The averageof ρ takenwith P(ρ) is, in general,different
from the equilibrium value ρǫ, given by the equation C. Comparison with the numerical experiment
(βf)′(ρ )=ǫ, (6.94)
ǫ Wehaveperformednumericalsimulationsaccordingto
cf. Eq. (2.6). This is a direct consequence of having the microscopic MonteCarlo dynamics described in Sec-
interpreted f(ρ) in (6.81) as the equilibrium free energy. tion V, in order to check its agreement with the order–
Therefore, we will only investigate the effect of noise by parameterequation. Wereportherethe resultsobtained
considering the linearized version of Eq. (6.81) close to for the two cases α = 2.5 (first–order phase transition)
equilibrium, i.e. and α = 1.5 (second–order phase transition). The mi-
croscopic dynamics is found to depend on the control
ρ˙ =−γ(βf)′′(ρ )(ρ−ρ )+ βη(t). (6.95) parameter v =ue−ǫ. Since ǫ=βµ, v can be used at the
ǫ ǫ
place of the temperature T ∼β−1 for describing numer-
Below Tc, f′′(ρ) is related to the copmpressibility, ical results. Without prejudice of generality one can fix
ab = 1. According to Eq. (2.13) the critical value v is
∂2f c
=k T/χ (6.96) given by the relation
∂ρ2 B
v =1+ζ (0) (6.104)
whichisawelldefinedquantityinthelimitN →∞. This c α
isnotthecaseaboveT ,wherethedeterministicequation
c Mostoftheresultsdiscussedinthissectionhavebeenob-
predicts the occurrence of gelation. In order to circum- tainedforlatticesizeN =5×103. Thetimeevolutionof
vent this difficulty, we assume that (6.95) and (6.96) are
ρ(t) has been extended over time spans ranging between
still valid above T , while we maintain the explicit de-
c 100τ and 600τ. In this subsection τ = N denotes the
pendence on N of both ρ and χ. The time–dependent
ǫ natural time unit of lattice updates (not to be confused
correlation function C (t) is defined as
ρ with the τ used in Eq.(6.99) ) . Fluctuations have been
C (t)=hρ(t+t′)ρ(t′)i−hρ(t+t′)ihρ(t′)i (6.97) smoothened by averaging over a large number of initial
ρ conditions (typically 104). The initial conditions have
where the bracket denotes the time average. The power been sampled among high density equilibrium states ob-
spectrum of ρ is given by tained for v =0.1 (ρ(0)≈0.9). For the special situation
concerningthecaseα=2.5andρ(0)<ρ theinitialcon-
c
1 +∞ ditionshavebeensampledbyattributingtoeachsitethe
P (ω)= eiωtC (t)dt. (6.98)
ρ 2π ρ value ”1” with probabilityρ(0) and”0” with probability
Z−∞
9
0.9
(1−ρ(0)) . Lacking an equilibrium state of finite den-
sity due tothe first–ordernature ofthe phasetransition, 0.8
this is a natural choice for sampling states of fixed ini-
0.7
tialdensityρ(0). Wehavealsoverifiedthatotherrecipes
can modify the duration of the transient evolution, but 0.6
ρ
we have not observedany difference in the long time be-
0.5
havior. Finally, we have also verified that the results do
notdependonthechoiceofboundaryconditions. Inpar- 0.4
ticular, the results reported in this subsection have been
0.3
obtained for fixed boundary conditions, where the first
(last) lattice site is put in contact with a fixed ”0” state 0.20 100 2τ00 300 400
on its right (left).
1. The case α=2.5 FIG. 2: The density of H-bonds ρ as a function of τ for
Thetimeevolutionofρ(t)hasbeenanalyzedforseveral α = 2.5 and v = 2.4, 2.45, 2.5, 2.53 (from top to bottom).
valuesofv inthe range[2.0,3.0],startingfromhighden- Notethat a linear decrease in timeis eventually approached,
sity equilibrium states. Close to the theoretical critical thus yielding denaturation after a finite time (gelation phe-
value, v = 2.341486..., the dynamics has been sampled nomenon).
c
over a time lapse up to 600τ, in order to obtain a re-
liable identification of the critical point. Despite finite
not unexpected: finite-size and memory effects cannot
size effects are expected to affect significantly MC dy-
allow for the identification of a single transition point
namics in models with long–range interactions like the
between the two dynamical regimes. This analysis could
LPS model, we have been able to obtain the numerical
be slightly improved by considering much larger lattices
estimatev =2.34±0.02,whichagreesverywellthethe-
c
and longer simulations, while maintaining at least the
oretical expectation. Moreover, for v > v (i.e. T > T )
c c
same number of averages over initial conditions. We did
we have also verified that ρ(t) initially evolves according
notproceedalongthisline,becausethelimitofourcom-
to the dynamics ρ˙ =−γ[ǫ −ǫ] (see Fig.1).
ρ
putationalcapabilitieshasalreadybeenreachedwiththe
0 values adopted in the reported simulations.
Thisscenariodoesnotchangealsoforv <v (T <T ).
c c
We still find a qualitative agreement with the theoreti-
−0.5
cal predictions based on the order parameter dynamics.
Specifically, for small times the dynamics is again ruled
∆ρ/∆τ
−1
by ρ˙ =−γ[ǫ −ǫ] (see Fig.3) .
ρ
−1.5 0
−−20.06 −0.05 −0.04 −0.03 −0.02 −0.01 0 −0.5
[ε − ε(ρ)]
∆ρ/∆τ
FIG. 1: ρ˙ ≡ ∆ρ/∆τ versus ǫ−ǫ(ρ) for α = 2.5 and v = −1
2.4, 2.45, 2.5, 2.53. Thedashedlineisthebest–fitoftheslope
(−γ ≈30), common to all values of v . It has been drawn to −1.5
guidetheeyes for appreciating theextension of thetransient
evolution, duringwhich the four lines overlap.
−−20.06 −0.05 −0.04 −0.03 −0.02 −0.01 0
[ε − ε(ρ)]
For larger values of time, ρ(t) decreases linearly in
time, as predicted by Eq.(6.86). This confirms the pres- FIG. 3: ρ˙ ≡ ∆ρ/∆τ versus ǫ−ǫ(ρ) for α = 2.5 and v =
ence of the phenomenon of gelation, according to which 2.2, 2.25, 2.3, 2.33. Thedashedlineisthebest–fitoftheslope
the systemreachesadenutaredstateinafinite time (see (−γ ≈30), common to all values of v . It has been drawn to
Fig.2). guidetheeyesfor appreciating theextension of thetransient
Wehavenotbeenabletocheckthequantitativeagree- evolution, duringwhich thefour lines overlap.
ment with the theoretical predictions reported in sub-
section A. In particular, we have not found an effective Afterwards, one observes the decay to a finite asymp-
criterion for locating the crossover between the initial totic value ρ > ρ (see Fig.4) . According to Eq.(6.87)
ǫ c
and the asymptotic dynamics of ρ(t). It is expected to sucha decayis predicted to be exponential, with the de-
occur atthe criticaltime t (ρ(t )=ρ ), but inMC sim- cay rate τ given in Eq.(6.88) .
c c c
ulations the crossover region extends over a long time Also in this case a clear quantitative verification is
lapse, which typically begins before and extends far be- prevented by the long lapse of time characterizing the
yond t . We want to point out that such a scenario is crossover between the transient and the asymptotic dy-
c
10
0.9
well with the theoretical expectation v = 3.6060508....
c
On the other hand, the crossover between the transient
0.8
and the asymptotic dynamics is found to extend much
more than in the case α = 2.5. For instance, even
0.7 very close to v the transient behavior ruled by the law
ρ c
ρ˙ =−γ[ǫ −ǫ]lastsoverafewunitsofτ (seeFig.6). Such
ρ
0.6 a linear region reduces the more v is far from v
c
0.5
−0.5
0.4
0 20 40 τ 60 80 100
−1
∆ρ/∆τ
FIG.4: Thedensity of H-bondsρas a function of τ for α=
2.5 and v = 2.2, 2.25, 2.3, 2.33 (from top to bottom). Note −1.5
thatthedecaytoafinitevalueρǫ slows downsignificantly as
v approaches vc =2.34±0.02. −2
−0.07 −0.06 −0.05 −0.04 −0.03 −0.02 −0.01 0
[ε −ε(ρ)]
namics.
Further peculiar behaviors of the microscopic dynam-
ics emerge for v >v andρ(0)<ρ . In this case anyini- FIG. 6: ρ˙ ≡ ∆ρ/∆τ versus ǫ−ǫ(ρ) for α = 1.5 and v =
c c
3.59, 3.60, 3.61. A linear dependence can be attributed to
tialconditioncannotcorrespondto anequilibriumstate.
thefirst few units of τ.
Accordingly,hystereticeffectsdetermineagrowthofρ(t)
from ρ(0), until a value close to ρ is approached. Then,
c
However, the long–time dynamics allows to identify
ρ(t) starts to decrease until a linear–in–time decay is
the occurrence of the gelation phenomenon for v > v
eventually approached, as predicted by Eq.(6.85) (as an c
(T >T ). In particular, ρ eventually exhibits a decay in
example, see Fig.5) . c
time which is faster than linear. In order to exemplify
0.8 this behavior in Fig.7 we show ρ(τ) for some values of v
muchlargerthanv . Theasymptoticbehaviorisruledby
c
a decay faster than a power-law,but also slowerthan an
0.7
exponential. This finding confirms the qualitative agree-
mentwiththepredicitonsoftheorderparameterdynam-
ρ0.6 ics.
0.5
0,1
0.4
0 5 1τ0 15 20
ρ
FIG. 5: The density of H-bonds ρ as a function of τ for 0,01
α = 2.5 and v = 2.6. The evolution has been obtained by
averaging over randomly seeded initial conditions of initial
densityρ(0)=0.8and0.4 . Inordertobetterappreciatethe
hysteretic effect associated with the first–order nature of the 0,001
1 τ 10
phasetransition(lowercurve),wehavereportedtheevolution
over a relatively short time scale. Note also that the initial
conditions affect only the transient evolution: both curves FIG. 7: The density of H-bonds ρ as a function of τ for
convergetothesameasymptoticdynamics,which eventually α = 1.5 and v = 3.8, 4.5, 6.0 (from top to bottom) . The
turnsouttoalinear–in-timedecay,yieldingthegelationphe- double–logarithmic scale shows an aymptotic decay in time
nomenon. which is faster than linear.
The phenomenon of gelation is recovered also in this Wewanttopointoutthatwehavereportedtheresults
case, as predicted on the basis of the order parameter forlargevaluesforv,becausethischoiceallowstoreduce
equation. Nonetheless, a quantitative analysis is pre- the duration of the crossover from the transient regime.
vented for the above mentioned reasons. In fact, for the same values of v reported in Fig.6 the
2. The case α=1.5 long–time behavior sets in for much larger values of τ.
Numericalanalysispredicts a criticalvalue ofthe con- Moreover, note that also in these cases finite size effects
trol parameter v =3.60±0.02. This result agrees quite donotallowtoreachacompletelydenaturatedstate(ρ=
c
