J. Mol. Bid. (1991) 217, 413-419
Computer
Simulation
of DNA
Supercoiling
Konstantin V. Klenin2, Alexander V. Vologodskii’, Vadim V. Anshelevich’ Alexander M. Dykhne2 and Maxim D. Frank-Kamenetskii’ ‘Institute of Molecular Genetics USSR Academy of Sciences, Moscow 123182, USSR ‘Branch
of I. V. Kurchatov Institute of Atomic Energy Troitslc, Moscow Region 142092, USSR
(Received 4 April
1990; accepted l? September 1990)
We treat supercoiled DNA within a wormlike model with excluded volume. A modified Monte Carlo approach has been used, which allowed computer statistical-mechanical simulations of moderately and highly supercoiled DNA molecules. Even highly supercoiled molecules do not have a regular shape, though with an increase in writhing the chains look more and more like branched interwound helixes. The averaged writhing ( Wr) z @7 ALk. The superhelical free energy F is calculated as a function of the linking number, Lk. The calculations have shown that the generally accepted quadratic dependence of F on Lk is valid for a variety of conditions, though it is by no means universal. Significant deviations from the quadratic dependence are expected at high superhelical density under ionic conditions where the effective diameter of DNA is small. The results are compared with the available experimental data.
1. Introduction
or even moderately supercoiled molecules, because the probability of their occurrence due to thermal motion is negligible. Several attempts have been undertaken to overcome this problem and calculate the parameters for highly supercoiled DNA molecules (Tan & Harvey, 1989; Hao & Olson, 1989). Here, we extend our previous Monte Carlo calculations to make it possible to generate supercoiled DNA molecules with arbitrary supercoiling. This enables one to visualize supercoiled molecules theoretically and to calculate their various characteristics. Most significantly, we could calculate, using this approach, the dependence of superhelix energy on the DNA linking number for a wide range of superhelical densities. We have studied this dependence for different duplex diameter values and have come to the conclusion that the functional relation between superhelical energy and linking number should strongly depend on ionic strength.
Numerous attempts have been made to supply a theoretical description of DNA supercoiling. To this end, a purely mechanical model, which treated DNA as an elastic rod, has been applied (Benham, 1977, 1983; Le Bret: 1979, 1984). However, due to serious mathematical difficulties, people either completely ignored the entropy effects or treated only very short chains. The more adequate approach of Shimada & Yamakawa (1985), who treated DNA within the wormlike model, is also applicable to very short chains only. Computer simulation of supercoiled DNA by the Monte Carlo method is a much more fruitful approach. It was first proposed by Benham (1978) and was realized by Vologodskii et al. (1979). Whereas in the first calculations the DNA molecule was modelled as a freely jointed, infinitely thin chain (Vologodskii et al., 1979; Le Bret, 1980; Chen, 1981), later on a much more adequate model, which treated DNA as a wormlike chain with excluded volume, was developed (Frank-Kamenetskii et al., 1985; Shimada & Yamakawa, 1988; Klenin et al., 1989). This allowed a quantitative treatment of the experimental data on equilibrium thermal distribution of topoisomers, yielding the torsional and bending rigidities of duplex (Klenin et al., 1989). In its traditional form, however, the Monte Carlo approach does not permit the simulation of highly
2. Methods of Calculation We rely on the procedure of Metropolis et al. (1953), which was used earlier to generate an equilibrium distribution of closed wormlike molecules over the writhing
number (Frank-Kamenetskii et al., 1985; Shimada & Yamakawa, 1988; Klenin et al., 1989). A closed DNA comprising n Kuhn statistical lengths is modelled as a spatial curve that consists of kn equal straight elementary 413
0022-2836/91/030413X)7
$03.00/O
0 1991 Academic Press Limited
K. V. Klenin
414 segments. Such a conformation
has the energy:
(1) where the summation is done over all the hinges connecting the elementary segments, k, is the Boltzmann constant, T is the absolute temperature, 8, is the angle (in radians) between 2 adjacent segments, and c( is the bending rigidity constant. The latter parameter was defined in such a way that the Kuhn statistical length corresponded to k segments (see Frank-Knmenetskii et al., 1985). An elementary step in the Metropolis procedure consists in the rotation of an arbitrary number of adjacent segments by a, random angle within the interval (-&,, f&J. The &, value is chosen in such a way that about half of the steps would be “successful”. Within t.he frame(Frankwork of the above standard procedure Kamenetskii et ~2.: 1985; Shimada & Yamakawa, 1988; Klenin et nl., 1989) one obtains a set of chains, where chains with large absolute values of writhing (Wr) occur very seldom. Such a sample makes it possible to determine the equilibrium distribution function over writhing; P(Wr), only in a comparatively narrow range of Wr (several standard deviations Jm). The P(Wr) function is connected with the free energy F(Wr) and the energy of a given microstate, E, in the following way: P( Wr) cc exp (-F(
Wr)/&T)
cc fd{r) &WY-
Wr({rj)) exp (-E((r))lW).
(2)
Here, the integration embraces all the possible conformations of the chain {r}, Wr({r}) is a function of the chain conformation, UJr is just a number. The E((r}) value is defined by eqn (I). To calculate the P( Wr) function for an arbitrary range of the Wr value, we have used the following trick. Let US use instead of energy in the Metropolis procedure the following value: E&@-j) = E((r))+k,Tg(Wr((rj)),
et al, different Q values, we used the fact that the sampies for different Q overlapped. The F(Wr) functions for 2 values of p were shifted until the overlapping regions match. The writhing calculations were performed as described (Vologodskii et al., 1979; Le Bret, 1980). The excluded volume effects were allowed for by defining the DXA effective diameter, d. Within the framework of the XiIetropolis procedure, this meant that a chain was assumed to have infinite energy if 2 non-adja.cent segments in it were closer t,ogether than d. The knotted chains were discriminated in a similar way. The analysis of chains with respect to knots was done as described by Frank-Kamenetskii & Vologodskii (1981). Our calculations were performed for values of n of 5 and 9 and for values of d of 0 to Cl (in terms of Kuhn statistical length, b). The parameter lc deserves special consideration. The k value should be sufficiently large to prevent any change of the results with its further increase. We chose k = 10 for n = 5 and k = 5 for n = 9. Our previous experience (Klenin et aE., 1989) indicated that, this choice guarantees a reliable estimation of the writhing variance, (( Wr)‘). However, this does not, necessarily mean a low statistical error for large Wr values. Unfortunately, since the computer time sharply increased with k, we could not check that our data for the F(Wr) function corresponded to saturation with respect to this parameter for all the WY values under study. We performed only the following control computer experiment. En the case of n = 5, we repeated all calculaLions for k = 6. We observed syst,ematic deviations from the results presented below, which correspond to k = 10; namely that the F value is overestimated at high writhing, the effect being especially significant for small d. However, in any case the deviation does nob exceed 10%. We therefore believe that our actual error is not higher than this. Table 1 summarizes the CJvalues that we used. Five elementa,ry steps in the Metropolis procedure required about, 1 s on a 386 microprocessor. For given values of parameters, from 3 x lo4 to 1.5 X lo5 conformations were generated.
(3) 3. Results
where g( Wr) is an arbitrary function (see below). Substituting E&(r)) instead of E({r}) in eqn (2), one obtains the distribution function: P,(Wr) cc jd(r)
x6(Wr-Wr({r))) xexp(-E((r})lk,T-g(Wr((rj))).
(4)
Due to the delta-function under the integration symbol, one can put the Wr number instead of the Wr((r}) function in the exponent and remove the exp (g( Wr)) function from under the integration symbol. As a result one arrives at: P,( Wr) cc exp ( - F( JYr)/k,T - g( VT)).
(5)
(a) S’hupe of mpercoibeci
molecules
Figure 1 shows typical conformations of chains with different writhing numbers. As expected, even highly supercoiled molecules do not have a regular shape. However, one can clearly see that with increasing Wr the molecules acquire the features of the branched int,erwound helix. We never observed a conformation that would resemble a chain wound around a torus. Note that sharp angles in the chains in Figure 1 are due to an insufficiently high number
Hence, the expression for the free energy is: F( Wr) = - k,T(Jn [P.( Wr)] +g( Wr)) + const.
Table 1
(6)
The maximum of the P,( Wr) function thus obtained is shifted as compared with the maximum of the equilibrium distribution function P( Wr). The range of the Wr values for which the P,, Wr) function can be determined reliably is shifted accordingly. Throughout this paper we use the linear function g( Wr) = -q Wr: where p is a calculation parameter. The shift is regulated by changing the (I parameter. Through ca.lculations for different p values one can determine the shape of the F( WT) function over different ranges of the Wr value. To splice the resu1t.s obtained for
The q valu,es used in this study ,
a
n = 5, k = 10 (a = 24030)
d/b 000
002
0.06
0.10
PO 2.3
0.0 4.5
O-0 5.0 85
O-0 7.5 13.0
n = 9, k = 5 (a = 1.!381) O-00 0.02
0.06
0.0 2.4
0.0 3.2 6.0 8.0
0.0 2.5 50
0.10 0.0 35 7.0 103 lxl
Computer Simulation of straight segments in the Kuhn statistical length. However, further increase in the subdivision of t.he chain would lead to a substantial increase of computational time without significant increase in the accuracy of our calculations (see Methods of Calculation). (b) Superhelical
ALk = ATw+
-wr=o.93
9
;.-.,::1;
-Wr=24
415
Supercoiling -wr=5.31
-Wr = 6.10
-Wr = 6.34
-wr=7.01
energy
We are starting with the well-known relation between the linking difference, ALk, the change of number, Wr twisting, ATw, and the writhing (White, 1969; Fuller, 1971): As before (Vologodskii
qf DNA
et al.,
Wr. 1979; Klenin
(7) et al.,
-wr = I.66
/c ; -Wr=2.84
Figure 1. Typical conformations of chains obtained in the course of computer simulation. The number of Kuhn statistical lengths n = 9 (corresponds to 2650 bp). The values of writhing are indicated; d/b = 0.02.
1989), we rely on the simplest assumption about the independence of bending and torsional deformations of the duplex in the open DNA forms. The agreement between theory and experiment achieved by Klenin et al. (1989) confirms this assumption. Equation (7) leads to the following expression for the equilibrium distribution function over the linking difference: P(AL@ = Jl P(ATw)P( Wr) 6(ALk:-ATwWr)d(ATw)d( WY),
-Wr=3.56
(8)
where P(ATw) and P( Wr) are the distribution functions of mutually independent variables, ATw and Wr, for non-closed DNA just before closing. Note that experimentally the equilibrium distribution over the linking difference is achieved by treating closed circular DNA with type I topoisomerases, which introduce and remove single-stranded breaks, or by ligating nicked DNA circles (Depew & Wang, 1975; Pulleyblank et al., 1975; Horowitz & Wang, 1984). The ATw value is the sum of small independent deviations from the equilibrium value of the winding angle between adjacent base-pairs. As a result, the distribution of the ATw value is Gaussian : P(ATw) cc exp
-
(ATw)’ ~((ATw)~)
1 (9) ’
where the variance, ((ATw)‘), depends on the DNA length, L (in cm), and on the value of the torsional
416
K. V. Klenin 40
et al.
r i/
d/b = 0.02
d/b =O.OO
30-
20-
2 ;: a
0
< I
d/b = 0.06
IO I
i/,
i 0
0
1234012345
, , , I/,, 2
4
6
80
, ,i 2
4
6
0
Wr
Wr
Figure 2. The dependence of the superhelix free energy on writhing at different values of the effective diameter of the duplex, d, for chains with n = 5. Different points correspond to different Q values (see Table 1).
Figure 3. The dependence of the super-helix free energy on writhing at different values of the effective diameter of the duplex, d. for chains with n = 9. Different points correspond to different q values [see Table 1).
rigidity 1985);
numerically. The resulting helical energy on the linking in Figures 4 and 5.
constant
C
(Frank-Kamenetskii
et al.,
The C value has been reliably estimated earlier as 3 x IO-l9 erg cm. Thus, the P(ATw) function is known. The distribution function over writhing, P( Wr), and the corresponding free energy, F( Wr), were calculated as described in Methods of Calculation. We simulated DNA molecules contained five and nine Kuhn statistical lengths (1470 and 2650 base-pairs, respectively: the Kuhn statistical length being equal to 100 nm (Hagerman, 1981; Klenin et al., 1989)). We have also varied the effective diameter of DNA, d; which strongly depends on ionic strength. Note that the Kuhn statistical length is virtually ionic strength-independent within a wide range of ionic strengths (Hagerman, 1981). Figures 2 and 3 show the results presented in the form of F( Wr). The superhelical free energy F(ALk) was calculated from equation (8) and: F(ALk) The
integration
= -k,T in
equation
In P(ALk). (8) was
1’) performed
dependence of superdifference is presented
4. Discussiora We have eIabora.ted a Monte Carlo approach that makes it possible to simulate moderately and highly supercoiled molecules. This opens new possibilities for understanding and analysing different aspects of DNA supercoiling. First of all, we now have an idea as to what DNA molecules with a physiological superhelical density look like. Figure 1 demonstrates that, although their conformations significantly vary, all of them may be classified as branched interwound helixes. Their shape resembles electron microscopy patterns of supercoiled DNA molecules (Boles et al., 1990). Our method of simulation makes it possible to calculate different physical characteristics of closed circular DNA with the fixed number of supercoils (AL,k), such as gyration radius, scattering function, etc. To this end, one has to calculate the chosen physical characteristic, x, a,s a function of Wr. Then the x( WY) value should be averaged over the writhing number, provided that the ALk value is
Computer
Simulation
of DNA
417
Supercoiling loo:
80 -
d/b=O.OO
!’
,
:
,:
d/b =0*02
/ :
60 -
a < 8o-
d/b=O-06
20
--I__ 2
0
4
6
0
2
4
6
A 8
Figure 4. The dependence of the superhelix free energy on the linking difference for different values of the effective diameter of the duplex, d, for chains with n = 5. The broken line corresponds to quadratic approximation (see eqn (15)).
fixed,
which
xWW
3
6
9
12
0
3
6
9
12
Figure 5. The dependence of the superhelix free energy on the linking difference for different values of the effective diameter of the duplex, d, for chains with n = 9. The broken line corresponds to quadratic approximation (see eqn (15)).
leads to the expression:
11s = $,(ALk)
x( Wr)P(ATw)P(
h(ALk-ATw-
Wr)
Wr)d(ATw)d( WY). (12)
We have used the above procedure for the case of x( Wr) = Wr to calculate the (Wr) value as a function of ALk. Figure 6 shows the results for two values of the DNA chain length and for different values of the DNA diameter. These data show that, for long enough DNA molecules, the (WY) value is close to @7 ALk for a very wide range of ALk values. Boles et aZ. (1990) studied this problem experimentally and arrived at a remarkably similar quantitative result. Note that equation (12), with the aid of equations (2) and (9), can be rewritten as: x(ALJ-9 = ~.fx(Ir))exp -Wk-
~ 0
where l/Z is the normalizing factor. The equation (13) is a function of the chain tion. Thus, to calculate the averaged value given ALk value one should use the
D
c
0.41
(-E(irl)lW
W~(~r}))2/(2((ATw)2)))~~r~,
0.6.
0
(13) 1 value in conformaof x for a standard
2
4
6
8
IO
12
14
Figure 6. The results of calculation of averaged writhing as a function of ALk for n = 5 (top) and n = 9 diameter, d, are: (bottom). The values of DNA C (d,‘b = O-06); B (d/b = @02); A (d/b = @OO); D (d/b = 0.10).
K. V. Klenin
418 Metropolis E&(r))
procedure
= E({r))+k,T(AL~-
with
the energy: Wr((r)))2/ (X@~W)~)).
(14)
Equation (14) has a very simple physical meaning. It is the sum of bending and torsional energy of DNA with given values of ALk and Wr (Fuller, 1971; Hao & Olson, 1989). Note that Hao 6t Olson (1989) used a Metropolis-Monte Carlo procedure similar to ours to calculate the equilibrium conformation of supercoiled DNA. They focused, however, only on the shape of supercoiled DNA in its global elastic energy minimum. By contrast, our approach makes it. possible to calculate real measurable properties of highly supercoiled molecules, such as free energy, gyration radius, diffusion coefficient, etc. Here, we have concentrated mainly on calculations of the free energy of supercoiling, F(ALk). This significant characteristic of DNA supercoiling has been studied before, both theoretieahy (Vologodskii et al., 1979; Shimada & Yamakawa, 1985; Klenin et aE., 1989) and experimenta,lly (Depew & Wang, 1975; Pulleyblank et al., 1975; Horowitz & Wang, 1984), chiefly for low absolute values of superhelical studies density. Only the very first experimental 1970; Hsieh & Wang, 1975) (Bauer & Vinograd, measured superhelical energy at a, high negative superhelical density by studying binding with intercalating dyes. Unfortunately, t.hese studies were done at extremely high ionic strengths and the accuracy of measurement was not very high. The detailed shape of the F(ALk) function at high ALk values is essential to the study of conformational transitions in closed circular DNA induced by For theoretical interpretation negative supercoiling. of these data, people use the standard quadratic expression:
F(ALk) = ksT(ALk)2/(2((Ar,k)2)),
(15)
which was obtained from the data on thermal distribution of supercoiled DNA over ALk at low ALk values, making its validity for highly supercoiled molecules questionable. The results of our calculations presented in Figures 4 and 5 show that, under a variety of conditions, equa$ion (15) is valid for high superhelica1 densities. However, this is by no means the general case. For small values of the DNA effective diameter, significant deviations from equation (15) are predicted at high values of superhelical density. A small DNA diameter means the screening of the Coulomb repulsion potential and corresponds to a higher ionic strength. Therefore, one can expect that, at a high concentration of a monovalent salt or in the presence of bivalent cations, the dependence of F on ALk at high superhelical densities would be much weaker than predicted by equation (15). Tf this were really the case, one could expect a serious revision of the quantitative estimations of energy parameters for conformational transitions in supercoiled D1\‘A induced by supercoiling, which were obtained from equation (15).
et al.
However, the most dramatic effect of nonquadratic energy dependence should be expected on the width of structural transitions induced by negative supercoiling. Indeed, the width of transit,ion 6Lk depends on the second derivative of superhelix energy over ALk taken at the transition point (Kozyavkin et al., 1990):
6Lk = (4k,T/6Tw)(d2F(ALk)/d(ALk)2)-1,
(16)
where 6Tw is the change of the Tw value in unconstrained DNA that accompanies the transition. Our numerical calculat.ions show that with the changing DNA diameter, the d’F(ALk)/d(ALk)’ value may change more than fivefold. This agrees reasonably with data obtained by Koz.yavkin et at. (199(l), who reported a dramatic broadening of the transition zone. We hope that our theoretical perdictions will stimulate extensive experimental studies of different. properties of highly supercoiled DNA molecules. We thank N. R. Cozzarelli for fruitful
discussion.
References Baser, W. & Vinograd, ,J. (1970). J. No/. Biol. 47. 419-435. Bedlam, C. J. ( 1977). Proc. Nat. Acad. Sci., U.S. A 74, 2397-2401. Benham, C. J. (1978). J. h’ol. Biol. 123, 361-370. Benham, C. J. (1983). Biopolymers, 22, 2477-2495. 5oles, C. D., White, J. & Cozzarelli, N. R. (1990). J. :!&ot. Biol. 213, 931-951. Chen, Y. (1981). J. Chem. Phys. 75, 2447-2453. Depew , R. E. & Wang, J. C. ( 1975). Proc. ivat. Aead. Sci., U.S.A. 72, 4275-4279. Frank-Kamenetskii, M. D. & Vologodskii, A. V. (1981). Sov. Phys.-Usp. 24, 679-696. Frank-Kamenetskii, M. D., tukashin, A. V., Anshelevich, V. V. & Vologodskii, A. V. (1985). J. Biomol. Struct. Dynum. 2, 1005-1012. Fuller, F. B. (1971). Proc. Nat. Acud. Sci., U.S.A. 68. 815-819. Hagerman, P. J. (1981). Biopolymers, 20; 1583-1535. Hao, M.-H. & Olson, W. K. (1989). Macromolecules, 22; 3292-3303. Horowitz, D. S. Br.Wang, J. C. (1984). J. iLIo2. Biol. 173, 75-91. Hsieh, T.-S. & Wang, J. C. (1975). Biochemistry, 14: 527-535. Klennin, K. V.; Vologodskii, A. V., Anshelevich, V. V. Klisko, V. Y., Dykhne, A. M. & Frank-Kamenetskii, M. D. (1989). J. Biomol. Struct. Dynam. 6, 707-714. Kozyavkin, S. A., Slesarev. A. I., Malkhosyan, S. R. & Pa.nyutin, I. G. (1990). Eur. J. Bioch.em. 191, 105-l 13. Le Bret,, M. (1979). Biopolymers, 18, 1709-1725. Le Bret,. &l. (1980). Biopolymers, 19, 619-637. Le &et, WI. (1984). Biopolymers, 23, 1835-1867. Metropolis, N., Rosenbluth, A. W., Rosenbluth, ,M. N., Teller. A.. H. & Teller, E. (1953). J. Chem. Phys. 21, 5087-1092. Pulleyblank, D. E. & Shure, M., Tang, ID., Vinograd, J. & Vosberg; H.-P. (1975). Proc. Nat. Acad. Sci., U.S.A. 72, 4280-4284. Shimada, J. & Yamakawa, H. (1985). J. Mol. Biol. 184, 319-329.
Computer Simulation Shimada, J. & Yamakawa, H. (1988). Biopolymers, 27, 657-673. Tan, R. K. Z. & Harvey, S. C. (1989). J. Mol. Biol. 205, 573-591.
of DNA Supercoiling
419
Vologodskii, A. V., Anshelevich, V. V., Lukashin, A. V. & Frank-Kamenetskii, M. D. (1979). Nature (London), 280, 294-298. White, (1969). Amer. J. Math. 91, 693-728.
Edited by P. won Hippel
Computer
Simulation
of DNA
Supercoiling
Konstantin V. Klenin2, Alexander V. Vologodskii’, Vadim V. Anshelevich’ Alexander M. Dykhne2 and Maxim D. Frank-Kamenetskii’ ‘Institute of Molecular Genetics USSR Academy of Sciences, Moscow 123182, USSR ‘Branch
of I. V. Kurchatov Institute of Atomic Energy Troitslc, Moscow Region 142092, USSR
(Received 4 April
1990; accepted l? September 1990)
We treat supercoiled DNA within a wormlike model with excluded volume. A modified Monte Carlo approach has been used, which allowed computer statistical-mechanical simulations of moderately and highly supercoiled DNA molecules. Even highly supercoiled molecules do not have a regular shape, though with an increase in writhing the chains look more and more like branched interwound helixes. The averaged writhing ( Wr) z @7 ALk. The superhelical free energy F is calculated as a function of the linking number, Lk. The calculations have shown that the generally accepted quadratic dependence of F on Lk is valid for a variety of conditions, though it is by no means universal. Significant deviations from the quadratic dependence are expected at high superhelical density under ionic conditions where the effective diameter of DNA is small. The results are compared with the available experimental data.
1. Introduction
or even moderately supercoiled molecules, because the probability of their occurrence due to thermal motion is negligible. Several attempts have been undertaken to overcome this problem and calculate the parameters for highly supercoiled DNA molecules (Tan & Harvey, 1989; Hao & Olson, 1989). Here, we extend our previous Monte Carlo calculations to make it possible to generate supercoiled DNA molecules with arbitrary supercoiling. This enables one to visualize supercoiled molecules theoretically and to calculate their various characteristics. Most significantly, we could calculate, using this approach, the dependence of superhelix energy on the DNA linking number for a wide range of superhelical densities. We have studied this dependence for different duplex diameter values and have come to the conclusion that the functional relation between superhelical energy and linking number should strongly depend on ionic strength.
Numerous attempts have been made to supply a theoretical description of DNA supercoiling. To this end, a purely mechanical model, which treated DNA as an elastic rod, has been applied (Benham, 1977, 1983; Le Bret: 1979, 1984). However, due to serious mathematical difficulties, people either completely ignored the entropy effects or treated only very short chains. The more adequate approach of Shimada & Yamakawa (1985), who treated DNA within the wormlike model, is also applicable to very short chains only. Computer simulation of supercoiled DNA by the Monte Carlo method is a much more fruitful approach. It was first proposed by Benham (1978) and was realized by Vologodskii et al. (1979). Whereas in the first calculations the DNA molecule was modelled as a freely jointed, infinitely thin chain (Vologodskii et al., 1979; Le Bret, 1980; Chen, 1981), later on a much more adequate model, which treated DNA as a wormlike chain with excluded volume, was developed (Frank-Kamenetskii et al., 1985; Shimada & Yamakawa, 1988; Klenin et al., 1989). This allowed a quantitative treatment of the experimental data on equilibrium thermal distribution of topoisomers, yielding the torsional and bending rigidities of duplex (Klenin et al., 1989). In its traditional form, however, the Monte Carlo approach does not permit the simulation of highly
2. Methods of Calculation We rely on the procedure of Metropolis et al. (1953), which was used earlier to generate an equilibrium distribution of closed wormlike molecules over the writhing
number (Frank-Kamenetskii et al., 1985; Shimada & Yamakawa, 1988; Klenin et al., 1989). A closed DNA comprising n Kuhn statistical lengths is modelled as a spatial curve that consists of kn equal straight elementary 413
0022-2836/91/030413X)7
$03.00/O
0 1991 Academic Press Limited
K. V. Klenin
414 segments. Such a conformation
has the energy:
(1) where the summation is done over all the hinges connecting the elementary segments, k, is the Boltzmann constant, T is the absolute temperature, 8, is the angle (in radians) between 2 adjacent segments, and c( is the bending rigidity constant. The latter parameter was defined in such a way that the Kuhn statistical length corresponded to k segments (see Frank-Knmenetskii et al., 1985). An elementary step in the Metropolis procedure consists in the rotation of an arbitrary number of adjacent segments by a, random angle within the interval (-&,, f&J. The &, value is chosen in such a way that about half of the steps would be “successful”. Within t.he frame(Frankwork of the above standard procedure Kamenetskii et ~2.: 1985; Shimada & Yamakawa, 1988; Klenin et nl., 1989) one obtains a set of chains, where chains with large absolute values of writhing (Wr) occur very seldom. Such a sample makes it possible to determine the equilibrium distribution function over writhing; P(Wr), only in a comparatively narrow range of Wr (several standard deviations Jm). The P(Wr) function is connected with the free energy F(Wr) and the energy of a given microstate, E, in the following way: P( Wr) cc exp (-F(
Wr)/&T)
cc fd{r) &WY-
Wr({rj)) exp (-E((r))lW).
(2)
Here, the integration embraces all the possible conformations of the chain {r}, Wr({r}) is a function of the chain conformation, UJr is just a number. The E((r}) value is defined by eqn (I). To calculate the P( Wr) function for an arbitrary range of the Wr value, we have used the following trick. Let US use instead of energy in the Metropolis procedure the following value: E&@-j) = E((r))+k,Tg(Wr((rj)),
et al, different Q values, we used the fact that the sampies for different Q overlapped. The F(Wr) functions for 2 values of p were shifted until the overlapping regions match. The writhing calculations were performed as described (Vologodskii et al., 1979; Le Bret, 1980). The excluded volume effects were allowed for by defining the DXA effective diameter, d. Within the framework of the XiIetropolis procedure, this meant that a chain was assumed to have infinite energy if 2 non-adja.cent segments in it were closer t,ogether than d. The knotted chains were discriminated in a similar way. The analysis of chains with respect to knots was done as described by Frank-Kamenetskii & Vologodskii (1981). Our calculations were performed for values of n of 5 and 9 and for values of d of 0 to Cl (in terms of Kuhn statistical length, b). The parameter lc deserves special consideration. The k value should be sufficiently large to prevent any change of the results with its further increase. We chose k = 10 for n = 5 and k = 5 for n = 9. Our previous experience (Klenin et aE., 1989) indicated that, this choice guarantees a reliable estimation of the writhing variance, (( Wr)‘). However, this does not, necessarily mean a low statistical error for large Wr values. Unfortunately, since the computer time sharply increased with k, we could not check that our data for the F(Wr) function corresponded to saturation with respect to this parameter for all the WY values under study. We performed only the following control computer experiment. En the case of n = 5, we repeated all calculaLions for k = 6. We observed syst,ematic deviations from the results presented below, which correspond to k = 10; namely that the F value is overestimated at high writhing, the effect being especially significant for small d. However, in any case the deviation does nob exceed 10%. We therefore believe that our actual error is not higher than this. Table 1 summarizes the CJvalues that we used. Five elementa,ry steps in the Metropolis procedure required about, 1 s on a 386 microprocessor. For given values of parameters, from 3 x lo4 to 1.5 X lo5 conformations were generated.
(3) 3. Results
where g( Wr) is an arbitrary function (see below). Substituting E&(r)) instead of E({r}) in eqn (2), one obtains the distribution function: P,(Wr) cc jd(r)
x6(Wr-Wr({r))) xexp(-E((r})lk,T-g(Wr((rj))).
(4)
Due to the delta-function under the integration symbol, one can put the Wr number instead of the Wr((r}) function in the exponent and remove the exp (g( Wr)) function from under the integration symbol. As a result one arrives at: P,( Wr) cc exp ( - F( JYr)/k,T - g( VT)).
(5)
(a) S’hupe of mpercoibeci
molecules
Figure 1 shows typical conformations of chains with different writhing numbers. As expected, even highly supercoiled molecules do not have a regular shape. However, one can clearly see that with increasing Wr the molecules acquire the features of the branched int,erwound helix. We never observed a conformation that would resemble a chain wound around a torus. Note that sharp angles in the chains in Figure 1 are due to an insufficiently high number
Hence, the expression for the free energy is: F( Wr) = - k,T(Jn [P.( Wr)] +g( Wr)) + const.
Table 1
(6)
The maximum of the P,( Wr) function thus obtained is shifted as compared with the maximum of the equilibrium distribution function P( Wr). The range of the Wr values for which the P,, Wr) function can be determined reliably is shifted accordingly. Throughout this paper we use the linear function g( Wr) = -q Wr: where p is a calculation parameter. The shift is regulated by changing the (I parameter. Through ca.lculations for different p values one can determine the shape of the F( WT) function over different ranges of the Wr value. To splice the resu1t.s obtained for
The q valu,es used in this study ,
a
n = 5, k = 10 (a = 24030)
d/b 000
002
0.06
0.10
PO 2.3
0.0 4.5
O-0 5.0 85
O-0 7.5 13.0
n = 9, k = 5 (a = 1.!381) O-00 0.02
0.06
0.0 2.4
0.0 3.2 6.0 8.0
0.0 2.5 50
0.10 0.0 35 7.0 103 lxl
Computer Simulation of straight segments in the Kuhn statistical length. However, further increase in the subdivision of t.he chain would lead to a substantial increase of computational time without significant increase in the accuracy of our calculations (see Methods of Calculation). (b) Superhelical
ALk = ATw+
-wr=o.93
9
;.-.,::1;
-Wr=24
415
Supercoiling -wr=5.31
-Wr = 6.10
-Wr = 6.34
-wr=7.01
energy
We are starting with the well-known relation between the linking difference, ALk, the change of number, Wr twisting, ATw, and the writhing (White, 1969; Fuller, 1971): As before (Vologodskii
qf DNA
et al.,
Wr. 1979; Klenin
(7) et al.,
-wr = I.66
/c ; -Wr=2.84
Figure 1. Typical conformations of chains obtained in the course of computer simulation. The number of Kuhn statistical lengths n = 9 (corresponds to 2650 bp). The values of writhing are indicated; d/b = 0.02.
1989), we rely on the simplest assumption about the independence of bending and torsional deformations of the duplex in the open DNA forms. The agreement between theory and experiment achieved by Klenin et al. (1989) confirms this assumption. Equation (7) leads to the following expression for the equilibrium distribution function over the linking difference: P(AL@ = Jl P(ATw)P( Wr) 6(ALk:-ATwWr)d(ATw)d( WY),
-Wr=3.56
(8)
where P(ATw) and P( Wr) are the distribution functions of mutually independent variables, ATw and Wr, for non-closed DNA just before closing. Note that experimentally the equilibrium distribution over the linking difference is achieved by treating closed circular DNA with type I topoisomerases, which introduce and remove single-stranded breaks, or by ligating nicked DNA circles (Depew & Wang, 1975; Pulleyblank et al., 1975; Horowitz & Wang, 1984). The ATw value is the sum of small independent deviations from the equilibrium value of the winding angle between adjacent base-pairs. As a result, the distribution of the ATw value is Gaussian : P(ATw) cc exp
-
(ATw)’ ~((ATw)~)
1 (9) ’
where the variance, ((ATw)‘), depends on the DNA length, L (in cm), and on the value of the torsional
416
K. V. Klenin 40
et al.
r i/
d/b = 0.02
d/b =O.OO
30-
20-
2 ;: a
0
< I
d/b = 0.06
IO I
i/,
i 0
0
1234012345
, , , I/,, 2
4
6
80
, ,i 2
4
6
0
Wr
Wr
Figure 2. The dependence of the superhelix free energy on writhing at different values of the effective diameter of the duplex, d, for chains with n = 5. Different points correspond to different Q values (see Table 1).
Figure 3. The dependence of the super-helix free energy on writhing at different values of the effective diameter of the duplex, d. for chains with n = 9. Different points correspond to different q values [see Table 1).
rigidity 1985);
numerically. The resulting helical energy on the linking in Figures 4 and 5.
constant
C
(Frank-Kamenetskii
et al.,
The C value has been reliably estimated earlier as 3 x IO-l9 erg cm. Thus, the P(ATw) function is known. The distribution function over writhing, P( Wr), and the corresponding free energy, F( Wr), were calculated as described in Methods of Calculation. We simulated DNA molecules contained five and nine Kuhn statistical lengths (1470 and 2650 base-pairs, respectively: the Kuhn statistical length being equal to 100 nm (Hagerman, 1981; Klenin et al., 1989)). We have also varied the effective diameter of DNA, d; which strongly depends on ionic strength. Note that the Kuhn statistical length is virtually ionic strength-independent within a wide range of ionic strengths (Hagerman, 1981). Figures 2 and 3 show the results presented in the form of F( Wr). The superhelical free energy F(ALk) was calculated from equation (8) and: F(ALk) The
integration
= -k,T in
equation
In P(ALk). (8) was
1’) performed
dependence of superdifference is presented
4. Discussiora We have eIabora.ted a Monte Carlo approach that makes it possible to simulate moderately and highly supercoiled molecules. This opens new possibilities for understanding and analysing different aspects of DNA supercoiling. First of all, we now have an idea as to what DNA molecules with a physiological superhelical density look like. Figure 1 demonstrates that, although their conformations significantly vary, all of them may be classified as branched interwound helixes. Their shape resembles electron microscopy patterns of supercoiled DNA molecules (Boles et al., 1990). Our method of simulation makes it possible to calculate different physical characteristics of closed circular DNA with the fixed number of supercoils (AL,k), such as gyration radius, scattering function, etc. To this end, one has to calculate the chosen physical characteristic, x, a,s a function of Wr. Then the x( WY) value should be averaged over the writhing number, provided that the ALk value is
Computer
Simulation
of DNA
417
Supercoiling loo:
80 -
d/b=O.OO
!’
,
:
,:
d/b =0*02
/ :
60 -
a < 8o-
d/b=O-06
20
--I__ 2
0
4
6
0
2
4
6
A 8
Figure 4. The dependence of the superhelix free energy on the linking difference for different values of the effective diameter of the duplex, d, for chains with n = 5. The broken line corresponds to quadratic approximation (see eqn (15)).
fixed,
which
xWW
3
6
9
12
0
3
6
9
12
Figure 5. The dependence of the superhelix free energy on the linking difference for different values of the effective diameter of the duplex, d, for chains with n = 9. The broken line corresponds to quadratic approximation (see eqn (15)).
leads to the expression:
11s = $,(ALk)
x( Wr)P(ATw)P(
h(ALk-ATw-
Wr)
Wr)d(ATw)d( WY). (12)
We have used the above procedure for the case of x( Wr) = Wr to calculate the (Wr) value as a function of ALk. Figure 6 shows the results for two values of the DNA chain length and for different values of the DNA diameter. These data show that, for long enough DNA molecules, the (WY) value is close to @7 ALk for a very wide range of ALk values. Boles et aZ. (1990) studied this problem experimentally and arrived at a remarkably similar quantitative result. Note that equation (12), with the aid of equations (2) and (9), can be rewritten as: x(ALJ-9 = ~.fx(Ir))exp -Wk-
~ 0
where l/Z is the normalizing factor. The equation (13) is a function of the chain tion. Thus, to calculate the averaged value given ALk value one should use the
D
c
0.41
(-E(irl)lW
W~(~r}))2/(2((ATw)2)))~~r~,
0.6.
0
(13) 1 value in conformaof x for a standard
2
4
6
8
IO
12
14
Figure 6. The results of calculation of averaged writhing as a function of ALk for n = 5 (top) and n = 9 diameter, d, are: (bottom). The values of DNA C (d,‘b = O-06); B (d/b = @02); A (d/b = @OO); D (d/b = 0.10).
K. V. Klenin
418 Metropolis E&(r))
procedure
= E({r))+k,T(AL~-
with
the energy: Wr((r)))2/ (X@~W)~)).
(14)
Equation (14) has a very simple physical meaning. It is the sum of bending and torsional energy of DNA with given values of ALk and Wr (Fuller, 1971; Hao & Olson, 1989). Note that Hao 6t Olson (1989) used a Metropolis-Monte Carlo procedure similar to ours to calculate the equilibrium conformation of supercoiled DNA. They focused, however, only on the shape of supercoiled DNA in its global elastic energy minimum. By contrast, our approach makes it. possible to calculate real measurable properties of highly supercoiled molecules, such as free energy, gyration radius, diffusion coefficient, etc. Here, we have concentrated mainly on calculations of the free energy of supercoiling, F(ALk). This significant characteristic of DNA supercoiling has been studied before, both theoretieahy (Vologodskii et al., 1979; Shimada & Yamakawa, 1985; Klenin et aE., 1989) and experimenta,lly (Depew & Wang, 1975; Pulleyblank et al., 1975; Horowitz & Wang, 1984), chiefly for low absolute values of superhelical studies density. Only the very first experimental 1970; Hsieh & Wang, 1975) (Bauer & Vinograd, measured superhelical energy at a, high negative superhelical density by studying binding with intercalating dyes. Unfortunately, t.hese studies were done at extremely high ionic strengths and the accuracy of measurement was not very high. The detailed shape of the F(ALk) function at high ALk values is essential to the study of conformational transitions in closed circular DNA induced by For theoretical interpretation negative supercoiling. of these data, people use the standard quadratic expression:
F(ALk) = ksT(ALk)2/(2((Ar,k)2)),
(15)
which was obtained from the data on thermal distribution of supercoiled DNA over ALk at low ALk values, making its validity for highly supercoiled molecules questionable. The results of our calculations presented in Figures 4 and 5 show that, under a variety of conditions, equa$ion (15) is valid for high superhelica1 densities. However, this is by no means the general case. For small values of the DNA effective diameter, significant deviations from equation (15) are predicted at high values of superhelical density. A small DNA diameter means the screening of the Coulomb repulsion potential and corresponds to a higher ionic strength. Therefore, one can expect that, at a high concentration of a monovalent salt or in the presence of bivalent cations, the dependence of F on ALk at high superhelical densities would be much weaker than predicted by equation (15). Tf this were really the case, one could expect a serious revision of the quantitative estimations of energy parameters for conformational transitions in supercoiled D1\‘A induced by supercoiling, which were obtained from equation (15).
et al.
However, the most dramatic effect of nonquadratic energy dependence should be expected on the width of structural transitions induced by negative supercoiling. Indeed, the width of transit,ion 6Lk depends on the second derivative of superhelix energy over ALk taken at the transition point (Kozyavkin et al., 1990):
6Lk = (4k,T/6Tw)(d2F(ALk)/d(ALk)2)-1,
(16)
where 6Tw is the change of the Tw value in unconstrained DNA that accompanies the transition. Our numerical calculat.ions show that with the changing DNA diameter, the d’F(ALk)/d(ALk)’ value may change more than fivefold. This agrees reasonably with data obtained by Koz.yavkin et at. (199(l), who reported a dramatic broadening of the transition zone. We hope that our theoretical perdictions will stimulate extensive experimental studies of different. properties of highly supercoiled DNA molecules. We thank N. R. Cozzarelli for fruitful
discussion.
References Baser, W. & Vinograd, ,J. (1970). J. No/. Biol. 47. 419-435. Bedlam, C. J. ( 1977). Proc. Nat. Acad. Sci., U.S. A 74, 2397-2401. Benham, C. J. (1978). J. h’ol. Biol. 123, 361-370. Benham, C. J. (1983). Biopolymers, 22, 2477-2495. 5oles, C. D., White, J. & Cozzarelli, N. R. (1990). J. :!&ot. Biol. 213, 931-951. Chen, Y. (1981). J. Chem. Phys. 75, 2447-2453. Depew , R. E. & Wang, J. C. ( 1975). Proc. ivat. Aead. Sci., U.S.A. 72, 4275-4279. Frank-Kamenetskii, M. D. & Vologodskii, A. V. (1981). Sov. Phys.-Usp. 24, 679-696. Frank-Kamenetskii, M. D., tukashin, A. V., Anshelevich, V. V. & Vologodskii, A. V. (1985). J. Biomol. Struct. Dynum. 2, 1005-1012. Fuller, F. B. (1971). Proc. Nat. Acud. Sci., U.S.A. 68. 815-819. Hagerman, P. J. (1981). Biopolymers, 20; 1583-1535. Hao, M.-H. & Olson, W. K. (1989). Macromolecules, 22; 3292-3303. Horowitz, D. S. Br.Wang, J. C. (1984). J. iLIo2. Biol. 173, 75-91. Hsieh, T.-S. & Wang, J. C. (1975). Biochemistry, 14: 527-535. Klennin, K. V.; Vologodskii, A. V., Anshelevich, V. V. Klisko, V. Y., Dykhne, A. M. & Frank-Kamenetskii, M. D. (1989). J. Biomol. Struct. Dynam. 6, 707-714. Kozyavkin, S. A., Slesarev. A. I., Malkhosyan, S. R. & Pa.nyutin, I. G. (1990). Eur. J. Bioch.em. 191, 105-l 13. Le Bret,, M. (1979). Biopolymers, 18, 1709-1725. Le Bret,. &l. (1980). Biopolymers, 19, 619-637. Le &et, WI. (1984). Biopolymers, 23, 1835-1867. Metropolis, N., Rosenbluth, A. W., Rosenbluth, ,M. N., Teller. A.. H. & Teller, E. (1953). J. Chem. Phys. 21, 5087-1092. Pulleyblank, D. E. & Shure, M., Tang, ID., Vinograd, J. & Vosberg; H.-P. (1975). Proc. Nat. Acad. Sci., U.S.A. 72, 4280-4284. Shimada, J. & Yamakawa, H. (1985). J. Mol. Biol. 184, 319-329.
Computer Simulation Shimada, J. & Yamakawa, H. (1988). Biopolymers, 27, 657-673. Tan, R. K. Z. & Harvey, S. C. (1989). J. Mol. Biol. 205, 573-591.
of DNA Supercoiling
419
Vologodskii, A. V., Anshelevich, V. V., Lukashin, A. V. & Frank-Kamenetskii, M. D. (1979). Nature (London), 280, 294-298. White, (1969). Amer. J. Math. 91, 693-728.
Edited by P. won Hippel
