- Email: [email protected]
Monte Carlo simulation of the kinetics of the helix—coil transition in a synthetic DNA
20 December 1985
CHEMICAL PHYSICS LETTERS
Volume 122, number 5
MONTE
CARLO
OF THE
KlNETICS
SIMULATION
OF THE
HELIX-COIL
TRANSITION
IN A SYNTHETIC
DNA
Values of the rue a~ which an alkrnaring copolymcric nucleic acid mclw ar~er a hzmpcn~urejump have been oblaincd wilh a Monle Carlo simulation. Using a Glauber-king scheme 10 model the dynamics. the helix-random coil transilion shows a pure “monodispersive” rslaxalion in accord with previous annlglical and experimenlal prsdiclions.
1. Introduction Since the early 1950’s when Pauling, and Watson and Crick, respectively, proposed the hydrogen-bonded structures for polypeptides and for DNA, the orderdisorder transition between helical and randomly coiled conformations of certain biological macromolecules in solution has been studied by many authors. The main efforts have been devoted to analyze a wide variety of equilibrium aspects of this reaction [l-6] _ Simultaneously, although less frequently, both theoretical [7-l I] and experimental [12--141 studies on the kinetic features of nccleic acid reactions have appeared in the iiterature. We are interested in these dynamic questions. We focus here on the relaxation rate of a helix-coil transition in a nucleic acid. The reaction is induced in our case by a thermal shift of the equilibrium. Theoretical studies have shown that this process can be understood as an evolution from the perturbed chemical equilibrium among particles confined to a onedimensional lattice, so that only two states are available for each particle [ 1,4,8,15] _These particles are actually the hydrogen bonds connecting the two strands of the double helix, and the two states avail0 009-2614/85/$ (North-Holland
0330 0 Elsevier Science Publishers B.V. Physics Publishing Division)
able stand for intact or broken bonds. The confmement of the particles to one dimension, and the fact that a perceptible advance in the kinetics involves many elementary breaking processes, lead to fundamental differences between this chemical equilibrium and ordinary ones. In addition, the new.equilibrium approached after a temperature jump is usually reached very rapidly, so that relaxation techniques seem to be specially appropriate to follow kinetically the behaviour of the system [7-g]_ A priori a whole distribution of relaxation times may contribute to a different extent in the “polydispersive” relaxation spectrum [9,16] _However, under appropriate conditions concerning the size of the system and the intensity of the applied perturbation, a singleexponential decrease can be identified as describing the whole relaxation [8,9]. This will be the case here, in accord with the experimental results available [ 12]_ Part of the aforementioned differences between the helix-coil equilibrium and ordinary ones lies in the need to introduce some degree of cooperativity, reflecting the influence of the neighbors on the state of a given unit. Melting seems mostly to be a cooperative phenomenon, the breaking of the bonds connecting one base pair facilitating the breaking of those of 459
Volume 122, number 5
CHEMICAL PHYSICS LEITERS
neighboring base pairs. Several generalizations of the one-dimensional Ising model with nearest-neighbor correlations, originally conceived to deal with the ‘one-dimensionaI ferromagnet, have been successfully used in this context [8,9,17]. In order to study the kinetics of the helix-coil transition we have adopted here the dynamical version of the Ising model proposed by Glauber [8]. A previous formulation of Glauber’s dynamics applied to nucleic acids is due to Gael [S] _He restricted his analysis to the very vicinity of the melting temperature where half of the total number of bonds are broken. Several approximations were then used to decouple the hierarchy of equations resulting from Glauber’s treatment. In order to avoid the simplifications involved in Goel’s scheme, we performed similar calculations, but via a Monte Carlo simulation [ 191 of Glauber’s dynamics_ The results of such a simulation are reported in this paper. They confirm the predicted analytical behaviour and fit the corresponding experimental results satisfactorily [12]. To retain what is the most essential in the thermal helix-coil transition, a particularly simple system consisting of a DNA with only one kind of base pairs was considered in ref. [8] _III this way, the difficulties arising from the intramolecular heterogeneity of natural DNA [20] could be overridden. In our approach here we have also restricted ourselves to a synthetic DNA. More concretely, we have chosen the alternating nucleic acid dAT. Unfortunately some difficulties appear when dealing with such copolymers. We refer to the appearance of some degree of intrachain folding: a single strand with alternative A and T nucleotides can bend back on itself, forming intrachain double-stranded helices [3,21]. This effect is more pronounced at high salt concentrations, and seems to be the origin of a distinct behaviour observed experimentally in the denaturation rate of dAT under temperature jumps starting very close to the melting point [12]. In order to be conveniently reflected in our simulation, this additional aspect would require some modifications on the Glauber-Ising scheme used here. Considering those situations for which denaturation is unaffected by intrachain folding, the most important aspect of the helix-coil transition in nucleic acids, i.e. the breaking of the interchain bonds, turns out to be appropriately described in our simulation here.
460
20 December 1935
2: The model and its Monte Carlo simulation In this section we begin by summarizing the Glauber-Ising model as it is used for a one-component nucleic acid. The DNA molecule may be schematically represented by a ladder whose rungs are identified as bonding complexes between a pair of bases. We introduce a parameter p for each complex, so that p takes the value +I if the complex is broken, and -1 if it is intact. Let the state for a given contiguration at time r, {p(t)), be denoted by the set of values /.rl (t), _._, am, where Nis the total number of units. The total number of broken complexes is then
This is a statistical-mechanical quantity subjected to an appropriate average with P( {p]; t), the probability distribution for a complexion {/.r} at time t, 92 0) = {X1 n(iP])P(CP];
r) -
(2)
Accordin:to (2), we must solve for P( {p); t) in order to describe properly the kinetic behaviour of% (t)_ If we let-w,-(e) be the probability per unit time that the ith cornFlex flips from the value 3 to -b, while the others remain momentarily fued, then we may write the time derivative of P( {p); t) as N
(3) The set of transition probabilities, {w), can be determined according to Glauber. In terms of a constant r. describing the time scale on which transitions take place, and to be fmed later on, we have WJo
= + TO1 11 - PL$ + + 703 -*)&1
fl=tanh.T,
y=tanh2U,
+ L$+dl
1
(4)
where the two basic dimensionless parameters of the model, namely J and U, appear for the fust time. The meaning of these two quantities can be understood by writing the Ising Hamiltonian of the system, with pe-
Volume 122, number 5
CHEMICAL PHYSICS LETTJXS
riodic boundary conditions,
pN+l = pl, as
U accounts for the nearest-neighbor stacking energy, and measures the cooperativity of the melting process [l-4,8,9]_ On the other hand, Dstands for the energy of an intact complex relative to a broken one when units are independent. U is usually assumed temperature-independent [ 1,4], and semi-empirically evaluated from thermodynamic quantities. In taking U in this way, certain aspects not included explicitly in the model, as for example end melting [3], or configurational entropy [ 1.31, have been compensated in an ad hoc fashion. The parameter J is, on the other hand, a temperature-dependent field. As in previous approaches [ 1,8] , we shall follow the assumption that it depends on temperature through J=a(T-
T,),
(6)
T, being the melting temperature_ This relation holds more precisely when close to T, _ The proportionality constant a is a system-dependent parameter to be fitted later by comparison with experimental results for the copolymeric nucleic acid investigated here. From an analytical point of view we would proceed further by writing the evolution equation satisfied for the average (b(r)>_ This equation, which follows from (3) and (4), contains terms proportional to pair correlation functions C+(t) am), which likewise depend on higher correlation functions. An exact solution of this hierarchy is rarely attained_ Therefore we should resort to several simplifications. By assuming that the system is in the vicinity of T, , Goel [8] obtained a first-order approximation for t+(t)), in terms of the zero&order values of the correlation functions calculated for an infinite system at the melting point. In this way he found for the average fraction of broken bonds f(t) = 7Z (t)/N,
f(r) = : {I + ewp& + Cw
I-0
ianh J
- -19h-J
1-
(7)
The static equilibrium value given by the Ising model calculated in the limit N-t 00
f=i
{I +sinhJ/[exp(~u)+sinh2J11/2~
(8)
20 Demmber 1985
can easily be recovered from the long-time Emit, t + m, of (7) within the J+ 0 approximation invoked to obtain (7). If one introduces a relaxation function defmed as usually by #(t) = v(t) -f(-)]/v(O) f(m)]. it is readily seen that Q(t) shows a typical exponential behaviour given by Q(r) = exp[-(1
- r)
t/7ol -
69
In order to go beyond Goel’s simplifications, we have performed a computer Monte Carlo simulation of this dynamical model, using the set of transition probabilities specified in (4). For the moment, r. was chosen unity for “one Monte Carlo step/particle” [ 19]_ The real value of this time scale will be found later on by comparison with the experimental results. The Monte Carlo technique applied to our system considers the site variable% as a set of bits in the memory of a computer taking values 0 or 1_ Once an initial condition has been specified, the system evolves
through the random selection of the individual sites. When one of these sites is selected, the computer compares its transition probability as given by (4) with a random number z, 0 < z < 1, which is generated. If w > z the value of the site variable is reversed, and otherwise the system is not modified. By repeating this procedure many times one obtains successive states of the system. Finally the averages of the timedependent quantities of interest are carried out over the complexions generated in this way.
3. Results of Monte CarIo simulations As mentioned in section 1, we have chosen an alternating copolymeric nucleic acid, dAT, with N= lo4 units, corresponding to a molecular weight of a few million. This value falls entirely within the range of observed sizes [3,7], and in particular it agrees with those experimentally investigated by Spatz and Baldwin [ 121. Moreover, this relatively large choice of N is justified since otherwise end effects would probably become significant and the statistical-mechanical treatment itself could be highly influenced by fluctuations. For the cooperative parameter we have chosen U= 2. It is the same value used by Goel, and at the melting temperature (T, = 320 K) it yields a stacking energy AC (AG = 4kB nr) = 5 kcal/mol. Only es461
Volume 122, number 5
timates of the order of magnitude of tbis quantity have been obtained experimentally [1,4], but our value falls within the range provided by other models [2--41 and indicated by. previous experiments [22] _ Finally .T, or more precisely the proportionality factor a, is the most involved parameter to be fmed. .We must think of it as being intimately related with the conditions under which the melting process takes place. Since we try to compare our simulated relaxations with those experimentally investigated by Spatz and Baldwin, we have chosen a in a way that adjusts our results to the experimental ones over the entire range analyzed. This was indeed possible with u = 0.024 K-l _In addition, a range of values between a = 0.010 and n = 0.028, depending on the system, can be thermodynamically argued according to the estimated values of AFf and AS involved in the helix-coil transition [l-3] - Our estimate falls within this range. The temperature jumps investigated cover the most significant region of the melting, i.e. the portion of the sigmoid cooperative helix-coil curves (see fig. 1) centered on T = T,, where small changes in temperature cause pronounced modifications in the number of broken bonds. Every initial configuration, corresponding usually to temperatures below T,, was selected through the appropriate run started at T, and lengthened enough to attain thermal equilibrium. Once a specific equilibrium configuration for a given initial temperature was so chosen, it was perturbed with a temperature jump, and we’ observed the kinetics of
T - T. Fig. l_ Equilibrium average fraction of broken bonds plotted versus temperature. The solid line corresponds to the analytical results for N+ - as given by (8). The tierisks are results of our simulation.
462
20 December 1985
CHEMICAL PHYSICS LmERS
* r_i_
0,-l
-1
-?I
n (A6
(t)
3bO
)
600
900
1200
t (MC51 Fig. 2. Rela?ration of the fraction of broken bonds after a temperature jump from TO = T, - 1 to Tf = Tm + 1, as observed in our simulation. ‘Ihe slope of the solid line. obtained from a linear regression, yields the corresponding relaxation rate. The time t is in Monte Carlo steps (MCS).
melting at the final temperature by processing the outputs for the fraction of broken bonds. A typical relaxation is shown in fig. 2, where we have plotted semilogarithmically the change in the fraction of broken bonds Af(r) E:f(m) -f(f) versus time. As one can see, the effect of the temperature jump is large and easily followed. Except in very few patiological cases to be discussed later, the observed kinetic curves are linear over a considerable extent of the reaction_ This proves that, at least within the range analyzed here, the relaxations are largely pure exponentials or “monodispersive”. The unique significant discrepancies encountered, arise, as expected, in the long-time limit when the unavoidable equilibrium fluctuations become more and more evident. In every run we allowed the system to reach the final equilibrium configurations from which we evaluated f(m)- These values compare satisfactorily with the predicted analytical ones given by (8). as shown in fig. l- On the other hand, by avoiding the aforementioned long-time tails, unimportant in respect to the kinetic behaviour of the denaturation, the slope of the linear part of such curves can easily be evaluated. This slope is defined usually as the melting rate. Ftially, regarding the few patholo8ical cases mentioned previously, we can realize that they correspond to relaxations comprised entirely belo,w T,, for example those -2orT, - 1. starting with To = T, - 3andTf=Tm In these cases the small values off(f) are highly distorted by the intrinSic numerical errors of the simulations.
Volume 122, number 5
CHEMICAL PHYSICS
LlXCERS
20 December
1985
T,- i,
Tr- To Fig. 3. Linear variation of log(l/T) versus Tf - To corresponding to our simulation results for: A.-To = Tm - 3; n , TO = Tm - 2;Z,To = Tm -1; X.TO = Tm.Theexperimental dataof Spatz and Baldwin corresponding to the values of To closest to ours, are plotted for comparison: a, TO = T, - 2.7; q, To = T, - 1.7 _The broken lines fitting the experimental results show the same linear behavior. ‘Gne is given in Monte Carlo steps/particle.
Fig. 4. Scmilogarlthmic plot of the relaxation rates versus Tf - Tm. showing that our simulation results, with the same symbols as in fig. 3. are independent of the initial temperature, falling on a common line. The experimental results for a wide range of To (a, To = Tm - 1.7; a, TO = Tm - 2.7; q, To= Tm --3.7;~, TO = T, -5_7;m,To = Tm -6.6)may also be grouped along this line. The experimental points for TO = Tm - 0.15. a, display a sharply distincr behavior as a consequence of the intrachain fokting appearing when TO = Tm-
According pect
that
to Spatz
the melting
and Baldwin rate would
the final temperature_ in fig. 3 on a logarithmic for the melting
[12]
depend
we could strongly
exon
To this end, we have plotted scale our calculated
values
rates
versus the temperature jumps. The values observed by Spatz and Baldwin [ 121 are also plotted for comparison_ Those arising from our simulations have been transformed to adapt our time scale in the computer with the laboratory one. In this logarithmic
plot,
fects are not included in our model which purely deals with interchain ones, generically responsible for the essentials of the helix-coil transition.
such a transformation
only
amounts
by a global translation, and the set of nearly parallel lines found in the simulation and in’the experiment compare well. The melting rate rises rapidly, linearly in this logarithmic scale, as the jump size is increased. On the other hand, over the entire range of initial temperatures studied here, the melting rates depend onIy on the fmal temperature. This can easily be concluded from fig. 4, where we have plotted melting rates versus the final temperature. All the results now fall on a common line which is also identified from the experimental results as the one corresponding to a wide, range of To. When To is much closer to T, and specially for high salt concentrations, a secondary dependence of the rate on To seems to originate in the intra; chain folding. This feature is not exhibited in OUTresults_ But this can easily be explained. Intrachain ef-
References
111E.W.
MontiolJ and N.S. Gael, Biopolymers 4 (1966) 855, and references therein; N.S. Goel and E.W. Montioll. Biopolymers 6 (1968) 731. 121D-M. Crothers and N.R. Kallenbach, J. Chcm. Phys. 45 (1966) 917. PI D. Poland-and H-A. Scheram, Theory of helix-coil transitions in biopolymcrs (Academic Press, New York, 1970). [41 M. Leunp. F.C. Cboo and B.Y. Tong, EGopolymers 16 (1977) 1233. 151 I. Morgenstem, K. Binder and A. Baumglirter. J. Chem. Phys. 69 (1978) 253. Y.L. Lyubchenko [61 A-V. Vologodskii, B-R. Myan. and M.D. Frank-Kamenetskii, J. Biomol. Struct. Dynam. 2 (1984) 131. I73 D.M. Crothers, J. Mol. Biol. 9 (1964) 712. [aI N-S, Coel, Biopolymcrs 6 (1968) 55. t91 G. Schwarz and J. Engel, Angew. Chem. 84 (1972) 615. 1101T- Tanaka, M. Suzuki and A. Wada, J. Chem. Phys. 58 (1973) 5707.
463
Volume
122. number 5
[ll]
A. Suyama and A.Wada,
[12]
H.Ch. 213.
[13] [14]
Spaa
CI+EMICAL Biopalymers
23 (1984)
and R.L. Baldwin, J. Mol. Biol. 11
409.
PHYSICS
LErTERs
[17]
F-M. Pohl and T-M. Jovin, J. Mol. BioL 67 (1972) 375. Y. Yakabe. T. Sane. N. Kure. K. Mimkami and T. Yasunaga.,Biopolymer 21 (1982) 1703. 1151 H. Reiss. D.A. McQuarrie. J.P. McTague and E-R. Cohen, J. Chem. Phys. 44 (1966) 4567. [16] A. Baumgk-tner and K. Binder, J. Stat. Phys. 18 (1978) 423.
Camp. (1985),
be published.
[18] R. Glaubex. J. Math. Phys. 4 (1963) 294. [ 191 K. Bind?. Advan. Phys. 23 (1974) 917; Monte Carlo methods in statistical physics (Springer. Berlin, 1979). [20] V. Blomtield, 6.00th~~ and 1. Tmoco Jr., Physical chemistry of nucleic acids (Harper and Row, New [21] [22]
York, 1974). M. GG, J. Phys. Sot. Japan 23 (1967) 597.. W-M. Hung and P.O.P. Ts’o. J. Mol. Biol. 16 (1966) 523.
464
-20 D&ember-lb85
J.L. VallCs and F. Sagu&. Acta Cientika
to
(1965)
-.:
CHEMICAL PHYSICS LETTERS
Volume 122, number 5
MONTE
CARLO
OF THE
KlNETICS
SIMULATION
OF THE
HELIX-COIL
TRANSITION
IN A SYNTHETIC
DNA
Values of the rue a~ which an alkrnaring copolymcric nucleic acid mclw ar~er a hzmpcn~urejump have been oblaincd wilh a Monle Carlo simulation. Using a Glauber-king scheme 10 model the dynamics. the helix-random coil transilion shows a pure “monodispersive” rslaxalion in accord with previous annlglical and experimenlal prsdiclions.
1. Introduction Since the early 1950’s when Pauling, and Watson and Crick, respectively, proposed the hydrogen-bonded structures for polypeptides and for DNA, the orderdisorder transition between helical and randomly coiled conformations of certain biological macromolecules in solution has been studied by many authors. The main efforts have been devoted to analyze a wide variety of equilibrium aspects of this reaction [l-6] _ Simultaneously, although less frequently, both theoretical [7-l I] and experimental [12--141 studies on the kinetic features of nccleic acid reactions have appeared in the iiterature. We are interested in these dynamic questions. We focus here on the relaxation rate of a helix-coil transition in a nucleic acid. The reaction is induced in our case by a thermal shift of the equilibrium. Theoretical studies have shown that this process can be understood as an evolution from the perturbed chemical equilibrium among particles confined to a onedimensional lattice, so that only two states are available for each particle [ 1,4,8,15] _These particles are actually the hydrogen bonds connecting the two strands of the double helix, and the two states avail0 009-2614/85/$ (North-Holland
0330 0 Elsevier Science Publishers B.V. Physics Publishing Division)
able stand for intact or broken bonds. The confmement of the particles to one dimension, and the fact that a perceptible advance in the kinetics involves many elementary breaking processes, lead to fundamental differences between this chemical equilibrium and ordinary ones. In addition, the new.equilibrium approached after a temperature jump is usually reached very rapidly, so that relaxation techniques seem to be specially appropriate to follow kinetically the behaviour of the system [7-g]_ A priori a whole distribution of relaxation times may contribute to a different extent in the “polydispersive” relaxation spectrum [9,16] _However, under appropriate conditions concerning the size of the system and the intensity of the applied perturbation, a singleexponential decrease can be identified as describing the whole relaxation [8,9]. This will be the case here, in accord with the experimental results available [ 12]_ Part of the aforementioned differences between the helix-coil equilibrium and ordinary ones lies in the need to introduce some degree of cooperativity, reflecting the influence of the neighbors on the state of a given unit. Melting seems mostly to be a cooperative phenomenon, the breaking of the bonds connecting one base pair facilitating the breaking of those of 459
Volume 122, number 5
CHEMICAL PHYSICS LEITERS
neighboring base pairs. Several generalizations of the one-dimensional Ising model with nearest-neighbor correlations, originally conceived to deal with the ‘one-dimensionaI ferromagnet, have been successfully used in this context [8,9,17]. In order to study the kinetics of the helix-coil transition we have adopted here the dynamical version of the Ising model proposed by Glauber [8]. A previous formulation of Glauber’s dynamics applied to nucleic acids is due to Gael [S] _He restricted his analysis to the very vicinity of the melting temperature where half of the total number of bonds are broken. Several approximations were then used to decouple the hierarchy of equations resulting from Glauber’s treatment. In order to avoid the simplifications involved in Goel’s scheme, we performed similar calculations, but via a Monte Carlo simulation [ 191 of Glauber’s dynamics_ The results of such a simulation are reported in this paper. They confirm the predicted analytical behaviour and fit the corresponding experimental results satisfactorily [12]. To retain what is the most essential in the thermal helix-coil transition, a particularly simple system consisting of a DNA with only one kind of base pairs was considered in ref. [8] _III this way, the difficulties arising from the intramolecular heterogeneity of natural DNA [20] could be overridden. In our approach here we have also restricted ourselves to a synthetic DNA. More concretely, we have chosen the alternating nucleic acid dAT. Unfortunately some difficulties appear when dealing with such copolymers. We refer to the appearance of some degree of intrachain folding: a single strand with alternative A and T nucleotides can bend back on itself, forming intrachain double-stranded helices [3,21]. This effect is more pronounced at high salt concentrations, and seems to be the origin of a distinct behaviour observed experimentally in the denaturation rate of dAT under temperature jumps starting very close to the melting point [12]. In order to be conveniently reflected in our simulation, this additional aspect would require some modifications on the Glauber-Ising scheme used here. Considering those situations for which denaturation is unaffected by intrachain folding, the most important aspect of the helix-coil transition in nucleic acids, i.e. the breaking of the interchain bonds, turns out to be appropriately described in our simulation here.
460
20 December 1935
2: The model and its Monte Carlo simulation In this section we begin by summarizing the Glauber-Ising model as it is used for a one-component nucleic acid. The DNA molecule may be schematically represented by a ladder whose rungs are identified as bonding complexes between a pair of bases. We introduce a parameter p for each complex, so that p takes the value +I if the complex is broken, and -1 if it is intact. Let the state for a given contiguration at time r, {p(t)), be denoted by the set of values /.rl (t), _._, am, where Nis the total number of units. The total number of broken complexes is then
This is a statistical-mechanical quantity subjected to an appropriate average with P( {p]; t), the probability distribution for a complexion {/.r} at time t, 92 0) = {X1 n(iP])P(CP];
r) -
(2)
Accordin:to (2), we must solve for P( {p); t) in order to describe properly the kinetic behaviour of% (t)_ If we let-w,-(e) be the probability per unit time that the ith cornFlex flips from the value 3 to -b, while the others remain momentarily fued, then we may write the time derivative of P( {p); t) as N
(3) The set of transition probabilities, {w), can be determined according to Glauber. In terms of a constant r. describing the time scale on which transitions take place, and to be fmed later on, we have WJo
= + TO1 11 - PL$ + + 703 -*)&1
fl=tanh.T,
y=tanh2U,
+ L$+dl
1
(4)
where the two basic dimensionless parameters of the model, namely J and U, appear for the fust time. The meaning of these two quantities can be understood by writing the Ising Hamiltonian of the system, with pe-
Volume 122, number 5
CHEMICAL PHYSICS LETTJXS
riodic boundary conditions,
pN+l = pl, as
U accounts for the nearest-neighbor stacking energy, and measures the cooperativity of the melting process [l-4,8,9]_ On the other hand, Dstands for the energy of an intact complex relative to a broken one when units are independent. U is usually assumed temperature-independent [ 1,4], and semi-empirically evaluated from thermodynamic quantities. In taking U in this way, certain aspects not included explicitly in the model, as for example end melting [3], or configurational entropy [ 1.31, have been compensated in an ad hoc fashion. The parameter J is, on the other hand, a temperature-dependent field. As in previous approaches [ 1,8] , we shall follow the assumption that it depends on temperature through J=a(T-
T,),
(6)
T, being the melting temperature_ This relation holds more precisely when close to T, _ The proportionality constant a is a system-dependent parameter to be fitted later by comparison with experimental results for the copolymeric nucleic acid investigated here. From an analytical point of view we would proceed further by writing the evolution equation satisfied for the average (b(r)>_ This equation, which follows from (3) and (4), contains terms proportional to pair correlation functions C+(t) am), which likewise depend on higher correlation functions. An exact solution of this hierarchy is rarely attained_ Therefore we should resort to several simplifications. By assuming that the system is in the vicinity of T, , Goel [8] obtained a first-order approximation for t+(t)), in terms of the zero&order values of the correlation functions calculated for an infinite system at the melting point. In this way he found for the average fraction of broken bonds f(t) = 7Z (t)/N,
f(r) = : {I + ewp& + Cw
I-0
ianh J
- -19h-J
1-
(7)
The static equilibrium value given by the Ising model calculated in the limit N-t 00
f=i
{I +sinhJ/[exp(~u)+sinh2J11/2~
(8)
20 Demmber 1985
can easily be recovered from the long-time Emit, t + m, of (7) within the J+ 0 approximation invoked to obtain (7). If one introduces a relaxation function defmed as usually by #(t) = v(t) -f(-)]/v(O) f(m)]. it is readily seen that Q(t) shows a typical exponential behaviour given by Q(r) = exp[-(1
- r)
t/7ol -
69
In order to go beyond Goel’s simplifications, we have performed a computer Monte Carlo simulation of this dynamical model, using the set of transition probabilities specified in (4). For the moment, r. was chosen unity for “one Monte Carlo step/particle” [ 19]_ The real value of this time scale will be found later on by comparison with the experimental results. The Monte Carlo technique applied to our system considers the site variable% as a set of bits in the memory of a computer taking values 0 or 1_ Once an initial condition has been specified, the system evolves
through the random selection of the individual sites. When one of these sites is selected, the computer compares its transition probability as given by (4) with a random number z, 0 < z < 1, which is generated. If w > z the value of the site variable is reversed, and otherwise the system is not modified. By repeating this procedure many times one obtains successive states of the system. Finally the averages of the timedependent quantities of interest are carried out over the complexions generated in this way.
3. Results of Monte CarIo simulations As mentioned in section 1, we have chosen an alternating copolymeric nucleic acid, dAT, with N= lo4 units, corresponding to a molecular weight of a few million. This value falls entirely within the range of observed sizes [3,7], and in particular it agrees with those experimentally investigated by Spatz and Baldwin [ 121. Moreover, this relatively large choice of N is justified since otherwise end effects would probably become significant and the statistical-mechanical treatment itself could be highly influenced by fluctuations. For the cooperative parameter we have chosen U= 2. It is the same value used by Goel, and at the melting temperature (T, = 320 K) it yields a stacking energy AC (AG = 4kB nr) = 5 kcal/mol. Only es461
Volume 122, number 5
timates of the order of magnitude of tbis quantity have been obtained experimentally [1,4], but our value falls within the range provided by other models [2--41 and indicated by. previous experiments [22] _ Finally .T, or more precisely the proportionality factor a, is the most involved parameter to be fmed. .We must think of it as being intimately related with the conditions under which the melting process takes place. Since we try to compare our simulated relaxations with those experimentally investigated by Spatz and Baldwin, we have chosen a in a way that adjusts our results to the experimental ones over the entire range analyzed. This was indeed possible with u = 0.024 K-l _In addition, a range of values between a = 0.010 and n = 0.028, depending on the system, can be thermodynamically argued according to the estimated values of AFf and AS involved in the helix-coil transition [l-3] - Our estimate falls within this range. The temperature jumps investigated cover the most significant region of the melting, i.e. the portion of the sigmoid cooperative helix-coil curves (see fig. 1) centered on T = T,, where small changes in temperature cause pronounced modifications in the number of broken bonds. Every initial configuration, corresponding usually to temperatures below T,, was selected through the appropriate run started at T, and lengthened enough to attain thermal equilibrium. Once a specific equilibrium configuration for a given initial temperature was so chosen, it was perturbed with a temperature jump, and we’ observed the kinetics of
T - T. Fig. l_ Equilibrium average fraction of broken bonds plotted versus temperature. The solid line corresponds to the analytical results for N+ - as given by (8). The tierisks are results of our simulation.
462
20 December 1985
CHEMICAL PHYSICS LmERS
* r_i_
0,-l
-1
-?I
n (A6
(t)
3bO
)
600
900
1200
t (MC51 Fig. 2. Rela?ration of the fraction of broken bonds after a temperature jump from TO = T, - 1 to Tf = Tm + 1, as observed in our simulation. ‘Ihe slope of the solid line. obtained from a linear regression, yields the corresponding relaxation rate. The time t is in Monte Carlo steps (MCS).
melting at the final temperature by processing the outputs for the fraction of broken bonds. A typical relaxation is shown in fig. 2, where we have plotted semilogarithmically the change in the fraction of broken bonds Af(r) E:f(m) -f(f) versus time. As one can see, the effect of the temperature jump is large and easily followed. Except in very few patiological cases to be discussed later, the observed kinetic curves are linear over a considerable extent of the reaction_ This proves that, at least within the range analyzed here, the relaxations are largely pure exponentials or “monodispersive”. The unique significant discrepancies encountered, arise, as expected, in the long-time limit when the unavoidable equilibrium fluctuations become more and more evident. In every run we allowed the system to reach the final equilibrium configurations from which we evaluated f(m)- These values compare satisfactorily with the predicted analytical ones given by (8). as shown in fig. l- On the other hand, by avoiding the aforementioned long-time tails, unimportant in respect to the kinetic behaviour of the denaturation, the slope of the linear part of such curves can easily be evaluated. This slope is defined usually as the melting rate. Ftially, regarding the few patholo8ical cases mentioned previously, we can realize that they correspond to relaxations comprised entirely belo,w T,, for example those -2orT, - 1. starting with To = T, - 3andTf=Tm In these cases the small values off(f) are highly distorted by the intrinSic numerical errors of the simulations.
Volume 122, number 5
CHEMICAL PHYSICS
LlXCERS
20 December
1985
T,- i,
Tr- To Fig. 3. Linear variation of log(l/T) versus Tf - To corresponding to our simulation results for: A.-To = Tm - 3; n , TO = Tm - 2;Z,To = Tm -1; X.TO = Tm.Theexperimental dataof Spatz and Baldwin corresponding to the values of To closest to ours, are plotted for comparison: a, TO = T, - 2.7; q, To = T, - 1.7 _The broken lines fitting the experimental results show the same linear behavior. ‘Gne is given in Monte Carlo steps/particle.
Fig. 4. Scmilogarlthmic plot of the relaxation rates versus Tf - Tm. showing that our simulation results, with the same symbols as in fig. 3. are independent of the initial temperature, falling on a common line. The experimental results for a wide range of To (a, To = Tm - 1.7; a, TO = Tm - 2.7; q, To= Tm --3.7;~, TO = T, -5_7;m,To = Tm -6.6)may also be grouped along this line. The experimental points for TO = Tm - 0.15. a, display a sharply distincr behavior as a consequence of the intrachain fokting appearing when TO = Tm-
According pect
that
to Spatz
the melting
and Baldwin rate would
the final temperature_ in fig. 3 on a logarithmic for the melting
[12]
depend
we could strongly
exon
To this end, we have plotted scale our calculated
values
rates
versus the temperature jumps. The values observed by Spatz and Baldwin [ 121 are also plotted for comparison_ Those arising from our simulations have been transformed to adapt our time scale in the computer with the laboratory one. In this logarithmic
plot,
fects are not included in our model which purely deals with interchain ones, generically responsible for the essentials of the helix-coil transition.
such a transformation
only
amounts
by a global translation, and the set of nearly parallel lines found in the simulation and in’the experiment compare well. The melting rate rises rapidly, linearly in this logarithmic scale, as the jump size is increased. On the other hand, over the entire range of initial temperatures studied here, the melting rates depend onIy on the fmal temperature. This can easily be concluded from fig. 4, where we have plotted melting rates versus the final temperature. All the results now fall on a common line which is also identified from the experimental results as the one corresponding to a wide, range of To. When To is much closer to T, and specially for high salt concentrations, a secondary dependence of the rate on To seems to originate in the intra; chain folding. This feature is not exhibited in OUTresults_ But this can easily be explained. Intrachain ef-
References
111E.W.
MontiolJ and N.S. Gael, Biopolymers 4 (1966) 855, and references therein; N.S. Goel and E.W. Montioll. Biopolymers 6 (1968) 731. 121D-M. Crothers and N.R. Kallenbach, J. Chcm. Phys. 45 (1966) 917. PI D. Poland-and H-A. Scheram, Theory of helix-coil transitions in biopolymcrs (Academic Press, New York, 1970). [41 M. Leunp. F.C. Cboo and B.Y. Tong, EGopolymers 16 (1977) 1233. 151 I. Morgenstem, K. Binder and A. Baumglirter. J. Chem. Phys. 69 (1978) 253. Y.L. Lyubchenko [61 A-V. Vologodskii, B-R. Myan. and M.D. Frank-Kamenetskii, J. Biomol. Struct. Dynam. 2 (1984) 131. I73 D.M. Crothers, J. Mol. Biol. 9 (1964) 712. [aI N-S, Coel, Biopolymcrs 6 (1968) 55. t91 G. Schwarz and J. Engel, Angew. Chem. 84 (1972) 615. 1101T- Tanaka, M. Suzuki and A. Wada, J. Chem. Phys. 58 (1973) 5707.
463
Volume
122. number 5
[ll]
A. Suyama and A.Wada,
[12]
H.Ch. 213.
[13] [14]
Spaa
CI+EMICAL Biopalymers
23 (1984)
and R.L. Baldwin, J. Mol. Biol. 11
409.
PHYSICS
LErTERs
[17]
F-M. Pohl and T-M. Jovin, J. Mol. BioL 67 (1972) 375. Y. Yakabe. T. Sane. N. Kure. K. Mimkami and T. Yasunaga.,Biopolymer 21 (1982) 1703. 1151 H. Reiss. D.A. McQuarrie. J.P. McTague and E-R. Cohen, J. Chem. Phys. 44 (1966) 4567. [16] A. Baumgk-tner and K. Binder, J. Stat. Phys. 18 (1978) 423.
Camp. (1985),
be published.
[18] R. Glaubex. J. Math. Phys. 4 (1963) 294. [ 191 K. Bind?. Advan. Phys. 23 (1974) 917; Monte Carlo methods in statistical physics (Springer. Berlin, 1979). [20] V. Blomtield, 6.00th~~ and 1. Tmoco Jr., Physical chemistry of nucleic acids (Harper and Row, New [21] [22]
York, 1974). M. GG, J. Phys. Sot. Japan 23 (1967) 597.. W-M. Hung and P.O.P. Ts’o. J. Mol. Biol. 16 (1966) 523.
464
-20 D&ember-lb85
J.L. VallCs and F. Sagu&. Acta Cientika
to
(1965)
-.: