- Email: [email protected]
A simple model for DNA denaturation transition
Physica A 314 (2002) 607 – 612
www.elsevier.com/locate/physa
A simple model for DNA denaturation transition Maria Serena Causoa; ∗ , Barbara Coluzzib , Peter Grassbergerc a INFM-NEST
and Scuola Normale Superiore, I-56126 Pisa, Italy di Fisica, Universit!a di Roma La Sapienza, and INFN, Sez. di Roma I, I-00185 Roma, Italy c NIC - Forschungszentrum J, ulich, D-52425 J,ulich, Germany
b Dipartimento
Abstract We present a simple model of DNA in which many important aspects of the real system are disregarded, while few relevant interactions which drive the system towards a .rst-order phase transition, as experimentally observed, are taken into account. c 2002 Elsevier Science B.V. All rights reserved. PACS: 87.15.Aa; 64.60.Kw Keywords: DNA; Thermal denaturation transition
1. Introduction The nature of the DNA denaturation transition has been experimentally investigated over the years [1]. A multistep behaviour in light-absorption spectra as the double-stranded domains melt upon heating was observed already in the 1950s in heterogeneous DNA, containing both AT and GC pairs. This was a clear signal of a .rst-order thermal denaturation transition. On the other hand, it is well known that the collapse transition for homopolymers in dilute regime (the -transition) is second order. Therefore, a simple question can be raised: What is the relevant interaction driving the transition out of the critical -collapse behaviour and making it .rst order? Our aim is to look for a simple model of DNA in which the relevant interactions producing a .rst-order phase transition from single- to double-stranded phase are taken into account. In order to have a low-temperature double-stranded phase a certain degree of inhomogeneity in the interactions between base pairs of the two strands is needed. ∗
Corresponding author. E-mail address: [email protected] (M.S. Causo).
c 2002 Elsevier Science B.V. All rights reserved. 0378-4371/02/$ - see front matter PII: S 0 3 7 8 - 4 3 7 1 ( 0 2 ) 0 1 1 4 6 - 9
608
M.S. Causo et al. / Physica A 314 (2002) 607 – 612
The minimal inhomogeneity consists in modeling the system as a diblock copolymer in which monomers in diCerent blocks exert a short-range attractive interaction on each other, while monomers in the same block interact only via steric interactions. It was numerically shown that this model undergoes a second-order zipping transition [2]. Therefore, this minimal inhomogeneity is not enough to produce a .rst-order phase transition.
2. The model A higher degree of inhomogeneity can be introduced modeling the system as two interacting polymers in which each monomer corresponds to a base and has its complementary base at the same position in the other chain and considering a short-range attraction only between complementary bases [3]. The system can be therefore represented as a pair of N -step self-avoiding walks (SAWs) (!1 ; !2 ) starting at the same origin on a 3D cubic lattice. The Boltzmann weight associated with a generic con.guration is N −H=KT e (1) = (1 − !i1 ;!j2 )(1 − !i1 ;!j1 )(1 − !i2 ;!j2 ) exp !i1 ;!i2 ; i=0
i=j
where = − =kT ˆ , ˆ being the binding energy. The partition function can be written as ZN =
N
cN; n e n
(2)
n=0
with n = card{i | !i1 = !i2 ; i ¿ 0}, and cN; n represents the number of allowed con.gurations of two N -step SAWs with n contacts. There are two conFicting tendencies in the system. On the one hand, the system experiences an entropic gain in the molten phase when the two strands open up to form denaturated bubbles and freely moving ends, since the number of accessible con.gurations in phase space is larger than in the double-stranded phase. On the other hand, there is the tendency to build energetically favoured contacts between the two strands. In the thermodynamic limit N → ∞ the balance between these opposite tendencies leads when ∗ to a phase transition between the high-temperature swollen phase, and the low-temperature double-stranded phase. The crossover behaviour near the transition (i.e., for → ∗ , N → ∞ with N ( − ∗ ) .nite) can be described by ZN ( ) = ( )N N
∗
−1
(( − ∗ )N ) ;
(3)
where ( ) is the critical fugacity at , ∗ is the tricritical susceptibility exponent, (x) is a scaling function which is .nite at x = 0, and the crossover exponent. From the basic assumption that for ¿ ∗ and N → ∞ the system behaves as an N -step SAW,
M.S. Causo et al. / Physica A 314 (2002) 607 – 612
609
and for ¡ ∗ as a 2N -step SAW, it follows that: (x) ∼ |x|(−
∗
)=
;
(4)
where is the susceptibility exponent for self-avoiding walks in 3D [4]. Morover, for ( − ∗ ) positive and small, we have 1 1 − ˙ ( − ∗ )1= : ∗ ( ) ( )
(5)
The scaling of the number of contacts between the two chains, which is proportional to the energy, can be obtained by taking the derivative of the logarithm of ZN ( ) with respect to : ∗ ( − )1=−1 N for ¿ ∗ ; (6) nN ( ) ∼ N for = ∗ ; 1=( ∗ − ) for ¡ ∗ : Taking the fraction of contacts n=N as the order paramenter of the system, we see that a .rst-order transition is obtained for = 1. Taking once more the derivative with respect to , one .nds that the speci.c heat for large N displays a peak which scales as N 2−1 and is located at = ∗ + const=N . 3. Numerical results In order to determine the nature of the transition we studied our model by MonteCarlo simulations [3] using the pruned-enriched Rosenbluth method (PERM) [5] with Markovian anticipation. In the present case the algorithm was implemented in such a way that the two chains grow simultaneously. Following the PERM strategy, the whole system grows according to the Rosenbluth method, while con.gurations with very large=small weight are cloned=pruned. We studied the scaling behaviour of the probability distribution of contacts as well as the scaling behaviour of energy, speci.c heat and geometric quantities such as the end-to-end distance of the two chains. The probability distribution of contacts behaves at a critical point as PN (n) = N − p(n=N ) ;
(7)
with p(x) scaling function. At a typical .rst-order phase transition, instead, no scaling is expected and the distribution function has a two-peak structure, which becomes sharper as N increases. For our system we .nd that the probability distribution of contacts scales with exponent = 1 at the estimated critical value ∗ = 1:3413(4). The scaling form is shown in Fig. 1. A peak at n = 0 and a shoulder at n=N ≈ 0:5 can be observed. The valley connecting them does not increase in depth with N , but considering values of slightly below and slightly above ∗ one sees the maximum of P(n) jumping discontinuously from zero to the shoulder of the distribution at ∗ .
610
M.S. Causo et al. / Physica A 314 (2002) 607 – 612
NP (n)
10
N=500 N=1000 N=1500 N=2000 N=2500 N=3000
1
0.1
0
0.1
0.2
0.3
0.4 n/N
0.5
0.6
0.7
0.8
Fig. 1. Probability distribution of the contact numbers at the estimated critical value.
12 10
N=500 N=1000 N=1500 N=2000 N=2500 N=3000
3 2.5 NP (n)
8 NP (n)
3.5
N=500 N=1000 N=1500 N=2000 N=2500 N=3000
6 4
2 1.5 1
2
0.5 0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
n/N
n/N
0.6 N=500 0.55 N=1000 0.5 N=1500 N=2000 0.45 N=2500 N=3000 0.4 0.35 0.3 0.25 0.2 0.15 0.1 0.05 -10
160 140 120 100 C
/N
Fig. 2. Probability distribution of the contact numbers at = 0:999 ∗ and = 1:001 ∗ .
N=500 N=1000 N=1500 N=2000 N=2500 N=3000
80 60 40 20
-5
0 (ε-ε*) N
5
10
0 1.32 1.325 1.33 1.335 1.34 1.345 1.35 1.355 1.36 ε
Fig. 3. The average contact number per monomer plotted against ( − ∗ )N with = 1:3413 and = 1 (left), and the speci.c heat per base pair (right).
M.S. Causo et al. / Physica A 314 (2002) 607 – 612
611
Therefore, we are in the presence of a .rst-order phase transition with discontinuous order parameter which, nevertheless, displays a tricritical scaling (see Fig. 2). The left panel in Fig. 3 shows that the average number of contacts scales too, but with a crossover exponent compatible with = 1 (and therefore = 1), which is the signal of a .rst-order phase transition. Also the scaling of the peak of the speci.c heat and its location are in agreement with = 1. 4. What can be the reason of the scaling at = ∗ ? One can look at the system from an extended scaling point of view, de.ning two correlation lengths 1 and 2 . The .rst, which we call the geometrical correlation length, is identi.ed as the Flory radius of any of the two polymers, 1 = (!N1 − !01 )2 1=2 =
(!N2 − !02 )2 1=2 . It follows the scaling law 1 ∼ N ! both in the denaturated as in the double-stranded phase, with ! being the Flory exponent. The second, 2 , is the thermal correlation length, which is .nite in the double-stranded phase and grows as the denaturation transition is approached. It is de.ned as the mean diameter of the molten ‘bubbles’ or the end-to-end distance Rend = (!N1 − !N2 ) 1=2 . De.ning a thermal correlation length exponent !T with 2 ∼ ( ∗ − )−!T one can show that ! : (8) = !T Numerical results for the thermal correlation length are shown in Fig. 4. For .xed
¿ ∗ , Rend goes to a constant for N → ∞, but as ∗ the observable does not saturate, but grows as N ! . This shows that 2 de.nes indeed a divergent thermal correlation length (see Fig. 4 (left)). The estimate of the thermal exponent !T is obtained from the behaviour of the quantity limN →∞ R2end for ∗ shown in Fig. 4 (right). Our estimate for !T agrees with the Flory exponent ! and therefore Eq. (8) is satis.ed with = 1. The absence of a .nite correlation length at the phase transition is due to q = 3.3 q = 3.65 q = 3.77 q = 3.815 q = 3.825 q = 3.835 q = 3.85 q = 3.88
Rend2
1000
100 const / (ε-ε* )0.59
limN→∞ Rend
10000
q = 4.2 q = 4.7 q = 5.5
100
10
10
1
1
10
100 N
1000
1 0.001
0.01
0.1
1
ε-ε∗
Fig. 4. Average squared end-to-end distance R2end for various values of q = e , plotted against N on a (left), and thermal correlation length 2 =limN →∞ Rend plotted against − ∗ on a double logarithmic scale (right).
612
M.S. Causo et al. / Physica A 314 (2002) 607 – 612
the fact that the system is topologically one dimensional. The energy cost of melting together two adjacent denaturated bubbles does not increase with N , i.e., the surface tension which is responsible for the .nite correlation length at a .rst-order phase transition vanishes. Applying results for general polymer networks to extreme con.gurations allowed for the system, further evidence for a .rst-order phase transition was found in Ref. [6]. References [1] [2] [3] [4] [5]
R.M. Wartell, A.S. Benight, Phys. Rep. 126 (1985) 67. E. Orlandini, F. Seno, A.L. Stella, Phys. Rev. Lett. 84 (2000) 294. M.S. Causo, B. Coluzzi, P. Grassberger, Phys. Rev. E 62 (2000) 3958. S. Caracciolo, M.S. Causo, A. Pelissetto, Phys. Rev. E 57 (1998) R1215. P. Grassberger, Phys. Rev. E 56 (1997) 3682; H. Frauenkron, M.S. Causo, P. Grassberger, Phys. Rev. E 59 (1999) R16; S. Caracciolo, M.S. Causo, P. Grassberger, A. Pelissetto, J. Phys. A 32 (1999) 2931. [6] Y. Kafri, D. Mukamel, L. Peliti, Phys. Rev. Lett. 85 (2000) 4988.
www.elsevier.com/locate/physa
A simple model for DNA denaturation transition Maria Serena Causoa; ∗ , Barbara Coluzzib , Peter Grassbergerc a INFM-NEST
and Scuola Normale Superiore, I-56126 Pisa, Italy di Fisica, Universit!a di Roma La Sapienza, and INFN, Sez. di Roma I, I-00185 Roma, Italy c NIC - Forschungszentrum J, ulich, D-52425 J,ulich, Germany
b Dipartimento
Abstract We present a simple model of DNA in which many important aspects of the real system are disregarded, while few relevant interactions which drive the system towards a .rst-order phase transition, as experimentally observed, are taken into account. c 2002 Elsevier Science B.V. All rights reserved. PACS: 87.15.Aa; 64.60.Kw Keywords: DNA; Thermal denaturation transition
1. Introduction The nature of the DNA denaturation transition has been experimentally investigated over the years [1]. A multistep behaviour in light-absorption spectra as the double-stranded domains melt upon heating was observed already in the 1950s in heterogeneous DNA, containing both AT and GC pairs. This was a clear signal of a .rst-order thermal denaturation transition. On the other hand, it is well known that the collapse transition for homopolymers in dilute regime (the -transition) is second order. Therefore, a simple question can be raised: What is the relevant interaction driving the transition out of the critical -collapse behaviour and making it .rst order? Our aim is to look for a simple model of DNA in which the relevant interactions producing a .rst-order phase transition from single- to double-stranded phase are taken into account. In order to have a low-temperature double-stranded phase a certain degree of inhomogeneity in the interactions between base pairs of the two strands is needed. ∗
Corresponding author. E-mail address: [email protected] (M.S. Causo).
c 2002 Elsevier Science B.V. All rights reserved. 0378-4371/02/$ - see front matter PII: S 0 3 7 8 - 4 3 7 1 ( 0 2 ) 0 1 1 4 6 - 9
608
M.S. Causo et al. / Physica A 314 (2002) 607 – 612
The minimal inhomogeneity consists in modeling the system as a diblock copolymer in which monomers in diCerent blocks exert a short-range attractive interaction on each other, while monomers in the same block interact only via steric interactions. It was numerically shown that this model undergoes a second-order zipping transition [2]. Therefore, this minimal inhomogeneity is not enough to produce a .rst-order phase transition.
2. The model A higher degree of inhomogeneity can be introduced modeling the system as two interacting polymers in which each monomer corresponds to a base and has its complementary base at the same position in the other chain and considering a short-range attraction only between complementary bases [3]. The system can be therefore represented as a pair of N -step self-avoiding walks (SAWs) (!1 ; !2 ) starting at the same origin on a 3D cubic lattice. The Boltzmann weight associated with a generic con.guration is N −H=KT e (1) = (1 − !i1 ;!j2 )(1 − !i1 ;!j1 )(1 − !i2 ;!j2 ) exp !i1 ;!i2 ; i=0
i=j
where = − =kT ˆ , ˆ being the binding energy. The partition function can be written as ZN =
N
cN; n e n
(2)
n=0
with n = card{i | !i1 = !i2 ; i ¿ 0}, and cN; n represents the number of allowed con.gurations of two N -step SAWs with n contacts. There are two conFicting tendencies in the system. On the one hand, the system experiences an entropic gain in the molten phase when the two strands open up to form denaturated bubbles and freely moving ends, since the number of accessible con.gurations in phase space is larger than in the double-stranded phase. On the other hand, there is the tendency to build energetically favoured contacts between the two strands. In the thermodynamic limit N → ∞ the balance between these opposite tendencies leads when ∗ to a phase transition between the high-temperature swollen phase, and the low-temperature double-stranded phase. The crossover behaviour near the transition (i.e., for → ∗ , N → ∞ with N ( − ∗ ) .nite) can be described by ZN ( ) = ( )N N
∗
−1
(( − ∗ )N ) ;
(3)
where ( ) is the critical fugacity at , ∗ is the tricritical susceptibility exponent, (x) is a scaling function which is .nite at x = 0, and the crossover exponent. From the basic assumption that for ¿ ∗ and N → ∞ the system behaves as an N -step SAW,
M.S. Causo et al. / Physica A 314 (2002) 607 – 612
609
and for ¡ ∗ as a 2N -step SAW, it follows that: (x) ∼ |x|(−
∗
)=
;
(4)
where is the susceptibility exponent for self-avoiding walks in 3D [4]. Morover, for ( − ∗ ) positive and small, we have 1 1 − ˙ ( − ∗ )1= : ∗ ( ) ( )
(5)
The scaling of the number of contacts between the two chains, which is proportional to the energy, can be obtained by taking the derivative of the logarithm of ZN ( ) with respect to : ∗ ( − )1=−1 N for ¿ ∗ ; (6) nN ( ) ∼ N for = ∗ ; 1=( ∗ − ) for ¡ ∗ : Taking the fraction of contacts n=N as the order paramenter of the system, we see that a .rst-order transition is obtained for = 1. Taking once more the derivative with respect to , one .nds that the speci.c heat for large N displays a peak which scales as N 2−1 and is located at = ∗ + const=N . 3. Numerical results In order to determine the nature of the transition we studied our model by MonteCarlo simulations [3] using the pruned-enriched Rosenbluth method (PERM) [5] with Markovian anticipation. In the present case the algorithm was implemented in such a way that the two chains grow simultaneously. Following the PERM strategy, the whole system grows according to the Rosenbluth method, while con.gurations with very large=small weight are cloned=pruned. We studied the scaling behaviour of the probability distribution of contacts as well as the scaling behaviour of energy, speci.c heat and geometric quantities such as the end-to-end distance of the two chains. The probability distribution of contacts behaves at a critical point as PN (n) = N − p(n=N ) ;
(7)
with p(x) scaling function. At a typical .rst-order phase transition, instead, no scaling is expected and the distribution function has a two-peak structure, which becomes sharper as N increases. For our system we .nd that the probability distribution of contacts scales with exponent = 1 at the estimated critical value ∗ = 1:3413(4). The scaling form is shown in Fig. 1. A peak at n = 0 and a shoulder at n=N ≈ 0:5 can be observed. The valley connecting them does not increase in depth with N , but considering values of slightly below and slightly above ∗ one sees the maximum of P(n) jumping discontinuously from zero to the shoulder of the distribution at ∗ .
610
M.S. Causo et al. / Physica A 314 (2002) 607 – 612
NP (n)
10
N=500 N=1000 N=1500 N=2000 N=2500 N=3000
1
0.1
0
0.1
0.2
0.3
0.4 n/N
0.5
0.6
0.7
0.8
Fig. 1. Probability distribution of the contact numbers at the estimated critical value.
12 10
N=500 N=1000 N=1500 N=2000 N=2500 N=3000
3 2.5 NP (n)
8 NP (n)
3.5
N=500 N=1000 N=1500 N=2000 N=2500 N=3000
6 4
2 1.5 1
2
0.5 0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
n/N
n/N
0.6 N=500 0.55 N=1000 0.5 N=1500 N=2000 0.45 N=2500 N=3000 0.4 0.35 0.3 0.25 0.2 0.15 0.1 0.05 -10
160 140 120 100 C
Fig. 2. Probability distribution of the contact numbers at = 0:999 ∗ and = 1:001 ∗ .
N=500 N=1000 N=1500 N=2000 N=2500 N=3000
80 60 40 20
-5
0 (ε-ε*) N
5
10
0 1.32 1.325 1.33 1.335 1.34 1.345 1.35 1.355 1.36 ε
Fig. 3. The average contact number per monomer plotted against ( − ∗ )N with = 1:3413 and = 1 (left), and the speci.c heat per base pair (right).
M.S. Causo et al. / Physica A 314 (2002) 607 – 612
611
Therefore, we are in the presence of a .rst-order phase transition with discontinuous order parameter which, nevertheless, displays a tricritical scaling (see Fig. 2). The left panel in Fig. 3 shows that the average number of contacts scales too, but with a crossover exponent compatible with = 1 (and therefore = 1), which is the signal of a .rst-order phase transition. Also the scaling of the peak of the speci.c heat and its location are in agreement with = 1. 4. What can be the reason of the scaling at = ∗ ? One can look at the system from an extended scaling point of view, de.ning two correlation lengths 1 and 2 . The .rst, which we call the geometrical correlation length, is identi.ed as the Flory radius of any of the two polymers, 1 = (!N1 − !01 )2 1=2 =
(!N2 − !02 )2 1=2 . It follows the scaling law 1 ∼ N ! both in the denaturated as in the double-stranded phase, with ! being the Flory exponent. The second, 2 , is the thermal correlation length, which is .nite in the double-stranded phase and grows as the denaturation transition is approached. It is de.ned as the mean diameter of the molten ‘bubbles’ or the end-to-end distance Rend = (!N1 − !N2 ) 1=2 . De.ning a thermal correlation length exponent !T with 2 ∼ ( ∗ − )−!T one can show that ! : (8) = !T Numerical results for the thermal correlation length are shown in Fig. 4. For .xed
¿ ∗ , Rend goes to a constant for N → ∞, but as ∗ the observable does not saturate, but grows as N ! . This shows that 2 de.nes indeed a divergent thermal correlation length (see Fig. 4 (left)). The estimate of the thermal exponent !T is obtained from the behaviour of the quantity limN →∞ R2end for ∗ shown in Fig. 4 (right). Our estimate for !T agrees with the Flory exponent ! and therefore Eq. (8) is satis.ed with = 1. The absence of a .nite correlation length at the phase transition is due to q = 3.3 q = 3.65 q = 3.77 q = 3.815 q = 3.825 q = 3.835 q = 3.85 q = 3.88
Rend2
1000
100 const / (ε-ε* )0.59
limN→∞ Rend
10000
q = 4.2 q = 4.7 q = 5.5
100
10
10
1
1
10
100 N
1000
1 0.001
0.01
0.1
1
ε-ε∗
Fig. 4. Average squared end-to-end distance R2end for various values of q = e , plotted against N on a (left), and thermal correlation length 2 =limN →∞ Rend plotted against − ∗ on a double logarithmic scale (right).
612
M.S. Causo et al. / Physica A 314 (2002) 607 – 612
the fact that the system is topologically one dimensional. The energy cost of melting together two adjacent denaturated bubbles does not increase with N , i.e., the surface tension which is responsible for the .nite correlation length at a .rst-order phase transition vanishes. Applying results for general polymer networks to extreme con.gurations allowed for the system, further evidence for a .rst-order phase transition was found in Ref. [6]. References [1] [2] [3] [4] [5]
R.M. Wartell, A.S. Benight, Phys. Rep. 126 (1985) 67. E. Orlandini, F. Seno, A.L. Stella, Phys. Rev. Lett. 84 (2000) 294. M.S. Causo, B. Coluzzi, P. Grassberger, Phys. Rev. E 62 (2000) 3958. S. Caracciolo, M.S. Causo, A. Pelissetto, Phys. Rev. E 57 (1998) R1215. P. Grassberger, Phys. Rev. E 56 (1997) 3682; H. Frauenkron, M.S. Causo, P. Grassberger, Phys. Rev. E 59 (1999) R16; S. Caracciolo, M.S. Causo, P. Grassberger, A. Pelissetto, J. Phys. A 32 (1999) 2931. [6] Y. Kafri, D. Mukamel, L. Peliti, Phys. Rev. Lett. 85 (2000) 4988.