The Electrostatic Contribution to DNA Base-Stacking Interactions RICHARD A. FRIEDMAN and BARRY H O N K *

Department of Biochemistry and Molecular Biophysics, Columbia University, 630 West 168th Street, New York, New York 10032

SYNOPSIS

Base-stacking and phosphate-phosphate interactions in B-DNA are studied using the finite difference Poisson-Boltzmann equation. Interaction energies and dielectric constants are calculated and compared to the predictions of simple dielectric models. No extant simple dielectric model adequately describes phosphate-phosphate interactions. Electrostatic effects contribute negligibly to the sequence and conformational dependence of base-stacking interactions. Electrostatic base-stacking interactions can be adequately modeled using the Hingerty screening function. The repulsive and dispersive Lennard-Jones interactions dominate the dependence of the stacking interactions on roll, tilt, twist, and propellor. The Lennard-Jones stacking energy in ideal B-DNA is found to be essentially independent of sequence.

INTRODUCTION The nature and cause of the sequence dependence of the conformation and flexibility of DNA is one of the most basic and important problems in molecular biophysics.lv2This is because the interaction of both beneficial and harmful ligands with DNA3 and the role of DNA in supramolecular structures4 depend on DNA sequence. Attempts to understand the sequence dependence of DNA conformation have usually focused on interactions between base pairs rather than on interactions involving the sugarphosphate backbone.'^^-'^ To understand the sequence dependence of DNA conformation, it is necessary to ascertain the relative extent to which different kinds of interactions determine the sequence and conformational dependence of the stacking energy. Some workers interpret DNA structure in terms of steric (repulsive) intera c t i o n ~while , ~ others find that both steric and dispersive (attractive ) interactions are important 1~13,94 and that electrostatic effects are negligible.13yg4 However, most workers find that electrostatic inBiopolymers, Vol. 32, 145-159 (1992) 0 1992 John Wiley & Sons, Inc.

CCC 0006-3525/92/020143-l5$04.00

* To whom correspondence should be addressed.

teractions play a significant or even predominant ro1e.6,9-12,14-25 M ost notably, two recent comprehensive, predictive, and somewhat successful base-only statistical mechanical theories of DNA bending have appeared: Sarai and co-workers find that electrostatic interactions between base pairs are the primary determinant of DNA structure, and play a major role in DNA fle~ibility’-’~;Olson and co-workers find that steric interactions are the primary determinant, but that electrostatic and dispersive interactions play a significant role.6-8 Any treatment of electrostatic interactions consists of a model of the charge distributions on the interacting molecules, and a dielectric model according to which the charge distributions interact.26 Various semiempirical, ab i n i t i ~ , ~and ~,” even experimental3’ methods of describing the charge distribution on DNA bases have been used. However, modeling the interactions between the charge distributions has been less developed It is usual practice to use a dielectric constant that is either a constant,6 or some simple function of distance.9,28,31-35 Such models do not take into account the effect on the electrostatic energy of the irregular boundary between the low dielectric solute and the high dielectric solvent that may contain Our laboratory has developed a method of calculat14~15~19,22927928

145

146

FRIEDMAN AND HONIG

ing electrostatic energies using the finite difference Poisson-Boltzmann ( FDPB ) equation that takes these effects into account.36This method, which is rigorous at the level of continuum electrostatics, has proven successful in predicting the pKas3' and redox potentials of proteins, 3 7 ~ 3 9and the solvation energies of ions4' and small molecule^.^^^^^ We have recently applied the method to DNA, focusing on phosphatephosphate interaction^.^^ The present paper extends this work to base-stacking interaction^.^^ Our goals are to test simple dielectric models, and to assess the relative role of electrostatic and Lennard-Jones ( L J ) interactions in determining the sequence and conformational dependence of stacking energies in DNA. The present paper is limited to consideration of base-base and phosphate-phosphate interactions. Interactions with solvent, which will form the subject of a subsequent paper, are only mentioned here where qualitative conclusions are affected.45 Dielectric saturation is not included in our calculations. Although water ordering has been observed in DNA there is no evidence at present that it will change the local dielectric constant of the first shell of water around DNA in sol ~ t i o n .Indeed, ~ ~ , ~the ~ nmr experiments of Halle and Piculell on dilute aqueous synthetic polyelectrolyte solutions48as further interpreted by Leyte4' demonstrate that the water reorientational anisotropy and hence the dielectric saturation in the vicinity of a polyion is minimal. As DNA is much wider and has a charge density only 40% higher than the polylons in their study, their result pertains to DNA as well. The success of our method in reproducing experimental and computer simulation results in small molecules and proteins further supports the neglect of dielectric saturation in treating In particular, Jayaram et al. have found that in the saturation of univalent ions, dielectric saturation is compensated by electrostriction, so that bulk dielectric constants remain ~ a l i d . ~ '

and co-workers, except that the only nonzero parameters are rise = 3.38 A and twist = 36". Hence roll, tilt, and propellor = 0". Unless otherwise specified, the parameters are based on 1980 fiber structure of B-DNA of Arnott and c o - w o r k e r ~ . ~ All ~,~',~~ bases in the step energy calculations have ideal BDNA parameters except as noted. Steps characterized by a nonzero propellor have the propellor imparted to both base pairs. Steps characterized by nonzero roll or tilt, or by twist # 36" have the deviation from ideal B conformation imparted to the 3' base pair alone. The generated structures were analyzed using the revised version of C ~ r v e sus~~-~~ ing a single linear axis with the helical parameters other than propellor of the 5' base pair of the step constrained to be those of ideal B-DNA. The parameters thereby obtained were very close to those used in this paper. The parameters used here are not invariant upon inversion of reading direction, and differ from parameters defined according to the mean-plane global or local Curves convention^.^^-^^ However, the definitions of helical parameters employed here are similar to those used by Olson and co-workers6 and by Sarai and co-workers." All steps were generated by altering the orientation of the central base pair only in an ideal BDNA heptamer. If the step has nonzero propellor, all of the base pairs in the heptamer have the same propellor as well. A minimum energy sugar-phosphate backbone corresponding to the fixed base coordinates was generated to provide an adequate description of the dielectric boundary. The minimization was performed using Jumna with a phosphate charge of -5, a bulk dielectric constant of 78, and a screening parameter slope of .16?3 These parameters were found by Lavery to give the most realistic structures for B-DNA.57The main role of the backbone in our base-stacking calculations is to help define the low dielectric region in the FDPB calculation of the electrostatic energy. Hence we expect the electrostatic energy to be relatively insensitive to the details of backbone conformation. In addition, the backbone energy minimization affects the torMETHODS sion angle of the thymine methyl group. The minimum energy orientation of the methyl group reflects Structures its torsion potential, and its van der Waals interThe structures of the interbase-pair steps analyzed action with its own base, surrounding bases, and the in this paper were generated using J ~ m n a ,except ~ ~ . ~ ~ sugar-phosphate backbone. The orientation of the as noted below. All helical parameters are defined methyl group is the only aspect of the interbase-pair with respect to the conventions used in the Jumna configuration that is not defined by the helical paprogram, except that the sign of the propellor is rerameters. Dielectric constant calculations were perversed to conform to the Cambridge convention^.^^ formed on 30-mers of ideal B-DNA based on the We define ideal B-DNA to have the same helical 1972 fiber structure of Arnott and c o - w ~ r k e r s ~ ~ parameters as the fiber B-DNA structure of Arnott generated using B~ilddna.~'The structure of

ELECTROSTATIC CONTRIBUTION T O BASE STACKING

CGCGAATTCGCG at 20°C obtained by Dickerson and coworkers60-62 was used. The structure was protonated using I n ~ i g h t I I . ~ ~ FDPB Calculations

The FDPB calculations of the electrostatic potential were performed using the Delphi program64as described p r e v i o ~ s l y . ~The ~,~ interior ~ , ~ ~ dielectric constant was taken to be 2, and the exterior dielectric constant was taken to be 78.5. A probe radius of 1.4 A was used. Energies were found not to vary significantly upon increase of the probe radius to 1.8 A. When the solution contained salt, the ion exclusion radius was 2 A.43Amber atomic radii were used, except that the hydrogen66radius was taken to be zero. It was found that results were insensitive to increasing this radius to 1 A. Calculations were done at several different grid units/A resolutions to ensure that the resolution was sufficiently high that artifacts due to placement on the grid were eliminated. In calculations involving phosphate-phosphate interactions, a resolution of 1.8grid units/A was found to be adequate, except for interactions between adjacent phosphates on the same chain, for which it was found necessary to use a resolution of 2.7 grid units/A. For phosphates farther apart than 22 A a resolution of .6 g/A was used. For base-stacking interactions the resolution was at least 3 grid units/ A. Focusing was performed to achieve this resolution when necessary. Debye-Huckel boundary conditions were used.36The presence of end effects was tested for by varying the interacting units on the oligomer. Interactions between units one or two base pairs from the ends have the same energy as equivalent units near the center of the oligomer. Therefore end effects are negligible. All calculations are done with the linear FDPB equation. The change in results on using the nonlinear equation in sample cases was found to be negligible. Although it is necessary to use the nonlinear equation for highly charged DNA,43the approximation introduced in the calculations in this paper, in which the maximum charge is -1, is small. Effective Dielectric Constants

The energy of interaction of two charged atoms, i and j , embedded in a constant dielectric medium in the vicinity of a dielectric discontinuity cannot be described by Coulomb's law using the dielectric constant of the medium they are in.26This is because the lines of force are perturbed by the presence of

147

the dielectric discontinuity. However, it is possible to characterize this interaction in terms of an effective dielectric constant ( or screening parameter ) E , , ~ .26.43 The quantity E ~ is , ~ the dielectric constant that, when used in Coulomb's law, gives the actual energy of interaction between the two charges. We calculate ee,/ by solving the linear FDPB equation with only atom i charged (with a unit charge). From this calculation we obtain 4LJ, the potential at atom j produced by the charge at atom i. Then eeY is given by

Note that for two interacting charges in a uniform dielectric medium of dielectric constant E , ceY = E. Since the linear Poisson-Boltzmann equation is a linear differential equation, $ Jis~ proportional ~ to qi. Therefore t," will be independent of qi . In the case of the phosphate-phosphate interactions, the effective dielectric constants were calculated based on the midpoints of the anionic oxygens. In the case of the interbase-pair interactions, the effective dielectric constants were calculated based on the atomic positions. The potential map was interpolated to give the potential at the interacting sites. The effective dielectric constants are also calculated according to the following simple dielectric models, all of which have been used in practice: 1. Vacuum modelz1: ceq = 1. (An effective dielectric constant of 4 has also been used.6) 2. Distance E,,) = r, (in A ) . 3. Debye-Huckel theory34,35,67: E,,, = 78.5/e-"'&~ 4. Hingerty model32:

(The Hingerty model is primarily a dielectric saturation model that incorporates dielectric discontinuity effects at short distances.) 5. Lavery te,J =

78 - 39 [ (sr,,)

+ 2sriJ + 21 e-srg) ( 3 )

(The Lavery model is based on the Hingerty model. Lavery found s = .16 to yield accurate structural parameters for B-DNA using This is also the value of s used by Jayaram et a1.68and Fenley et al.6') 6. Model of Sarai et a1.l': E ~ , =~ 5.01 {e-l 2 ( r , / - 3 . 0 ) 1 (The model of Sarai et al. is an exponentially screened Coulomb potential

148

FRIEDMAN AND HONIG

with the distance of closest approach of stacked bases subtracted from the interatomic distance. The authors tried various values of the exponential screening parameters and dielectric constant (numerator of the right-hand side of the above equation) in a statistical mechanical calculation of the structure of CGCGAATTCGCG. They found the above parameters to yield the closest match to the crystal structure.) Because the dielectric functions of Lavery51 and Sarai et al.'-'* are parameterized to reproduce experimental data, they may be viewed either as dielectric screening functions, or more generally as functions, which when used in conjunction with a charge set, constitute an effective potential of mean force that incorporates the effect of the solvent, including both dielectric and nondielectric effects. This consideration is especially relevant to the function of Sarai et al., which was parameterized and used in calculations in which the backbone is not includedg-12so that it may be thought to encode backbone as well as solvent effects. Hence the use of these functions in the modeling of DNA structure and in the calculation of electrostatic energies per se should be viewed as to some extent separate.

I

w

100

-

60

-

0

Stacking Interactions in DNA

FDPB stacking energies were calculated by computing the potential due to the charges of one base pair embedded in a DNA oligomer. The potential was obtained at the location of each of the atoms of the other base pair by interpolation from the potentials at the surrounding gridpoints. The electrostatic energy is then given by

where 4, is the potential at each of the atoms of base pair j due to the charges on all of the atoms at base pair i and the summation is over all of the atoms on base pair j . Energies of interaction according to the simple dielectric models were computed b y (5)

Unless otherwise stated, calculations were done with the all-atom AMBER charge set.66The charge set of Poltev and co-workers were also used in specific cases as a check.27In these two charge sets, the

.'

o.

1 0

@ .

2 0

30

Figure 1. Phosphate-phosphate effective dielectric constants vs distance in 1972 BDNA at 0 ionic strength. 0 , FDPB equation. Curve A, - - * - - * , ce,, = 78.5 (which is what Debye-Huckel theory reduces to at 0 salt). Curve B, -, cer, = r,,. Curve C, - - - -, Hingerty function, Eq. ( 2 ) . Curve D, - - - - - - -, Lavery function, Eq. ( 3 ) . The model of Sarai et al. is

- -

omitted because it is inappropriate to apply it to phosphate-phosphate interactions because it was parameterized for base-base interactions.

ELECTROSTATIC CONTRIBUTION T O BASE STACKING

charges are confined to the plane of the base. Possible effects of out-of-plane charges will be treated in a subsequent p~blication.~'

7.8(4.1)

H,

4.9(3.4)

149

,H 15.1(4.3) N I7.1(3.9)

c

4.6(3.4) N

18.(5.9)

5.3(3.7) 400 6.8(4.2)

6.7(4.1)

Figure 3. Effective dielectric constants between an adenine N1 atom and the atoms of an adenine 3' to it in 1972 fiber B-DNA. Numbers without parentheses, effective dielectric constants. Numbers in parentheses, interatomic distances (A).

a *

The LJ energy is given by

a a

where the -1 / r term represents the attractive (dispersive) interaction and the l / r term represents the repulsive term. Aij and Bij are constants.

c g

P

*

40 30 20

0

5

1 0

1 5

10

2 0

-

n 0

r.. (A)

1

2

3

4

5

6

7

8

9

'I

Figure 2. Phosphate-phosphate effective dielectric constants vs distance in 1972 fiber B-DNA at .15 ionic strength. 0 , FDPB equation. Curve A, - * - - - * , DebyeHuckel theory. Curve B, ---, tev= rij. Curve C, - - - -, Hingerty function, Eq. ( 2 ) . Curve D, - - - - - - -Lavery , function, Eq. ( 3 ) . The model of Sarai et al. is omitted because it is inappropriate to apply it to phosphate-phosphate interactions because it was parameterized for basebase interactions.

- -

Figure 4. Effective atom-atom dielectric constants describing the interaction of two adjacent adenines vs interatomic distance in 1972 fiber B-DNA at 0 ionic strength. 0 , FDPB equation. Curve A, - * - - - * , ce,, = 78.5 (which is what Debye-Huckel theory reduces to at 0 salt). Curve B, -, te,,= r , . Curve C, - - - -, Hingerty function, Eq. (2). Curve D, - - - - - - -,Lavery function, Eq. ( 3 ) . Curve E, . . . . . . . - ,function of Sarai et al.

-

. .

150

FRIEDMAN AND HONIG

Table I Base-Stacking Energy of Ideal B-DNA as a Function of Sequence and Dielectric Model (kcal/mol)

Step

FDPB

kJ= 1

AA AT TA GG GC CG TC CT AC CA Minimum Maximum Range

-0.18 0.33 -0.68 -0.06 -0.42 -0.80 -0.38 0.03 -0.13 -0.80 -0.80 0.33 1.14

0.65 2.54 -0.09 2.25 0.12 -2.38 0.54 1.28 1.05 -1.11 -2.38 2.54 4.92

teu

=

~ I J

-0.22 1.29 -0.67 0.48 -0.49 -1.84 -0.36 0.30 0.22 -1.34 -1.84 1.29 3.13

DebyeHuckel

Hingerty

0.01 0.03 0.00 0.03 0.00 -0.03 0.01 0.02 0.01 -0.01 -0.03 0.03 0.06

-0.16 0.53 -0.34 0.07 -0.27 -0.77 -0.23 0.08 0.06 -0.61 -0.77 0.53 1.30

The version of the Amber LJ parameters without the 10-12 interbase hydrogen-bond parameters was ~ s e d . ~Stacking ~ r ~ ~energies , ~ ~ were insensitive to inclusion or elimination of the 10-12 potential, including conformations having nonzero propellor. In specific cases, the LJ potential functions of Poltev and c o - w ~ r k e r were s ~ ~used ~ ~ ~as a check. The attractive component of the LJ energy is not

AA

Lavery

Sarai et al.

LennardJones

-0.74 1.84 -1.46 0.49 -1.25 -3.30 -0.96 0.24 -0.07 -2.62 -3.30 1.84 5.14

-0.71 1.90 -1.27 -0.32 -1.09 -2.56 -1.08 0.21 0.22 -2.20 -2.56 1.90 4.47

-15.99 -16.61 -15.47 -17.45 -17.96 -16.84 -16.74 -16.51 -17.16 -16.04 -17.96 -15.47 2.49

Interface Area

(A')

LJ Surface (cal/A2)

292.3 280.8 265.6 269.8 273.9 266.7 269.7 265.7 266.0 266.7 265.6 292.3 20.6

-54.7 -59.2 -58.2 -64.7 -65.6 -63.1 -62.0 -64.6 -67.8 -58.0 -65.6 -54.7 10.9

corrected for dielectric boundary effects because it results from very fast electronic motions.74 These motions are so fast that only the high-frequency electronic component of the electric susceptibility of the surrounding medium can respond to them.75 Since the high-frequency component is about the same for all molecular liquids and solids,76and is already present in the LJ parameters, no further

AT

Figure 5. Overlap of TA, AA, and A T steps in ideal B-DNA.

TA

ELECTROSTATIC CONTRIBUTION TO BASE STACKING

correction is necessary. The repulsive portion of the LJ interaction results from the local overlap of electron clouds and is therefore not subject to solvent effects.74The FDPB calculations take into account the dielectric response of the low dielectric DNA as well as the high dielectric medium. Hence, induced charge or dipole effects that appear as separate terms in quantum mechanical calculations 22 are automatically included in FDPB calculations. Changes in the LJ stacking energy accompanying a conformational change may be compensated for by interactions with the solvent. These interactions are estimated using a surface tension m ~ d e l . ~ The ~-~’ nonelectrostatic solute-solvent interaction is characterized by the microscopic alkane-water surface tension of 59 cal/ (mol A2).77 The curvature-corrected DNA solvent-accessible surface was calculated using the program GRASP7’ with a probe radius of 1.4 A.

-

RESULTS AND DISCUSSION Effective Dielectric Constants

0

the 1972 fiber structure5awere computed according to the FDPB equation and various simple dielectric models. These effective dielectric constants as a function of distance at 0 ionic strength are given in Figure 1and a t .15 ionic strength in Figure 2. It can be seen that none of the simple dielectric models, including ceG = 1 (not shown), accurately approximates the FDPB results, which are exact a t the level of continuum electrostatics. It may seem strange that effective dielectric constants larger than 78.5 were calculated by the FDPB method. This apparent anomaly may be understood as followsao: A fixed charge induces surface charge of opposite sign a t the dielectric boundary. This induced surface charge screens the potential due to the fixed charge felt by other fixed charges. In some cases, the screening is more effective than it would be if the fixed charges were embedded in a uniform medium of dielectric 78.5. Then the effective dielectric constant will be greater that 78.5.8’ The lack of correspondence between the Lavery screening function and the FDPB results is interesting in light of the Lavery function’s utility in Monte Carlo and counterion condensation

calculation^.^^^^^ All of the interatomic effective dielectric constants of two adjacent adenines in ideal B-DNA

The effective dielectric constants for phosphatephosphate interactions in ideal B-DNA based on I

151

I

I

I

I

I

I

I

I

I

1

I

I

I

I

I

I

I

1

2

3

4

5

6

7

8

9

I

I

I

I

1 0 1 1 1 2

Step Figure 6. Sequence stabilization energies of CGCGAATTCGCG. Step 1 is Cl-G2, step * ,LJ. 2 is G2-C3, etc. 0, -, electrostatic (FDPB). 0, *

------ -

152

FRIEDMAN AND HONIG

based on the 1972 fiber structure58at 0 ionic strength were calculated. The t,,, s characterizing the interaction of an adenine N1 atom and the atom of an adenine adjacent to it in the 3' direction are given in Figure 3. The eeslsare not a simple function of distance. The t,,]s of all of the atom-atom interactions between the two bases as a function of interatomic distance are given in Figure 4. The range of the te,]s is between 4 and 60. The Hingerty function fits the FDPB data better than the other models, although the fit is not especially good. Absolute Stacking Energies

The absolute stacking energies of the ten different kinds of ideal B-DNA steps a t 0 ionic strength according to the various dielectric models are given in Table I. Most of the electrostatic energies, when computed using the FDPB equation, are small, and their range is only 1kcallmol. This result indicates that charge-charge electrostatic stacking energies do not play a significant role in determining the sequence dependence of the stability of DNA. Of the simple dielectric models, the Hingerty function mimics the FDPB function most closely. The Hingerty function mimics the FDPB results more closely than is suggested by the comparison of dielectric constants in Figure 4. Apparently there is some cancellation of the factors associated with the individual charge-charge interactions. The Hingerty function is primarily derived from dielectric saturation, which, as discussed in the methods section we do not expect to occur significantly for DNA, while the FDPB equation calculates the effects of the dielectric boundary without including saturation effects. However, dielectric boundary effects and dielectric saturation both lead to an increase in the effective dielectric constant with increasing distance. Therefore, while the physical origin of the two ap-

Table I1 Electrostatic Contribution to the Sequence Stabilization Energy of CGCGAATTCGCG According to Various Dielectric Models Contribution

FDPB &,J =

1

-~ I J deb ye-Huckel Hingerty Lavery Sarai et al. tebJ

Energy (kcal/mol) .06 0.00

-.64 0.00 -.48

-.15 -3.14

proaches is different, their qualitative behavior is similar. It turns out that their quantitive consequences for base stacking is similar as well. The LJ contribution is strongly stabilizing, and is therefore dominated by the dispersion term. The range of the LJ interaction is 2.5 kcallmol. However, since the calculations omit LJ solvation energies, they cannot, by themselves, be used to predict the relative stabilities of DNA double helices of identical composition and different sequence. The base-pair-base-pair overlap areas were obtained by calculating the difference between the solvent-accessible surface of the isolated, backboneless base pairs and that of the backboneless step. These areas and the energy/A' of surface of the steps are also given in Table I. The average surface energy is -61 call (mol A'), in excellent agreement with the negative of the alkane work of cohesion," which is -60 call (mol A') .77 The areas of pyrimidine-purine steps are very close to those of purine-purine steps and purine-pyrimidine steps of identical composition. This result suggests that the experimentally observed undertwisting of AT steps and overtwisting of T A steps in alternating AT tracts" may not be caused by the steps having different overlaps in ideal B-DNAs3 because the associated energetic difference is only -1 kcallmol. Figure 5 shows that the notion that AT and TA steps have significantly different overlaps in ideal B-DNA is only apparent from molecular graphics if molecular volumes are not included. Indeed, a structure with an undertwisted T A step and an overtwisted AT step has been observed.'* The similarity of the LJ stacking energies of the AT and T A steps in ideal B-DNA conformation is also apparent from the results of several workers, 21,66*85 and has been noted explicitly in connection with the sequence dependence of the twist by Haran et al.I3 The step with the largest electrostatic stacking energy (in absolute value) is the CG step. When the ionic strength is raised from 0.0 to 1.0, the stacking energy remains at -230 kcallmol. The stacking energy is insensitive to ionic strength because most of the base-pair partial charges are far away from the solvent, unlike the phosphate charges, which are near the solvent. All calculations of stacking energy mentioned subsequently are done at 0 ionic strength. It is interesting to see if any single atom or group of atoms dominates the electrostatic interaction energy of two stacked bases. For the stacking of two adenines in ideal B-DNA based on the 1972 fiber structure, the total FDPB electrostatic energy is -.19 kcallmol. The total energy of interaction of each atom with all of the atoms on the adjacent base

-

ELECTROSTATIC CONTRIBUTION TO BASE STACKING

ranges between -1.13 and 1.79 kcal/mol. The atomatom interaction energies range between -11.52 and 12.52 kcal/mol. Hence the electrostatic energy is a small sum of positive and negative numbers that are often larger than the sum. Therefore, the interactions of no single atom or atom pair or small group of atoms or atom pairs dominate the electrostatic base-stacking energy.

The Sequence Stabilization Energy of CCCC AATTCGCC

Diffraction of fibers of random-sequence B-DNA yields the average structure of B-DNA. A crystal structure, on the other hand, represents the actual sequence-specific structure of DNA. The difference in energy between the crystal structure of a given

?

1

-E

n

0

\

Q

0

.Y

0

Y

% 0 L

!

I

Q)

C

W

-1

0

1

2

3

4

153

5

6

7

8

9 101112

Step Figure 7 . Sequence stabilization energies of CGCGAATTCGCG according to various of FDPB. +, * * * * * - . E ~ , , = 1. 0, . . * * * * * G" = rt,, +,. dielectric models. 0 , -, - - - - - , Debye-Huckel. m, - * - * - - * , Function of Sarai et al. 0, -, Hingerty function. A, -, Lavery function.

--- -

- - -

9

154

FRIEDMAN AND HONIG

sequence of DNA and the same sequence in the conformation of the fiber structure may be regarded as the energy that determines the sequence-specific conformation of DNA.13 We call this quantity the sequence stabilization energy. The size of the electrostatic and LJ contributions to the sequence stabilization energy is a measure of the extent to which each of these interactions determines the sequencedependent structure.13 The sequence stabilization stacking energy of CGCGAATTCGCG based on the crystal structure of Dickerson and co-workers and the 1980 fiber structure was calculated. (Note that the actual fiber structure and not an ideal structure based on it is used here.) The total (FDPB) electrostatic sequence stabilization energy is -.06 kcal/ mol. This result indicates that electrostatic stacking interactions contribute negligibly to the sequencespecific conformation of DNA. The total LJ sequence stabilization energy is -7.32 kcal/mol, a significant contribution. The contribution of the repulsive interactions to the sequence stabilization energy is -4.70 kcallmol. The contribution of attractive interactions is -2.62 kcal/mol. The change in the DNA-solvent surface area between the two structures was found to be 0 A'. Therefore, there are no destabilizing LJ interactions with the solvent

that compensate for the stabilizing LJ stacking interactions. The sequence stabilization energy on a stepwise basis is given in Figure 6. The magnitude of the electrostatic sequence stabilization energy is less than or equal to kT for all steps. Therefore, electrostatic stacking interactions do not determine the conformation of individual steps. The total sequence stabilization energy according to different dielectric models is given in Table 11. The sequence stabilization energy of individual steps according to different dielectric models is given in Figure 7. As before, the Hingerty model mimics the FDPB results reasonably well. Systematic Distortions of Ideal B-DNA

The effect of varying various structural parameters on stacking energies was studied. A series of AA steps were generated that have the ideal B-DNA conformation except that either roll, tilt, twist, or propellor was varied systematically. The steps were embedded at the center of ideal heptamers (except for steps with nonzero propellor, in which case each base pair has the same propellor). Similar series of structures of AT, TA, GG, GC, and CG steps were

10

-

a

n

0

-E

6

4

* El) L

a C w

2

0

-L

-30

-20

-10

0

10

20

30

Roll (0) Figure 8. Contributionsto the dependence of the stacking energy of an AA step on roll. Energies are relative to ideal B-DNA. a,-, electrostatic (FDPB). 0, * . * . LJ. A, - - - - - - -,electrostatic LJ.

+

- -- --

a ,

ELECTROSTATIC CONTRIBUTION T O BASE STACKING

also generated. For steps consisting of the same nucleotide, e.g., the AA step, oligomers in which the 5'-3' strand consisted only of that nucleotide, e.g., AT, were generated. For steps consisting of mixed nucleotides, alternating copolymers were generated, e.g., ATATATA. The change in the LJ energy and the electrostatic energy according to different dielectric models were calculated. The dependence of the LJ and electrostatic (FDPB) contributions to the stacking energy of an

AA step, and the sum of these contributions on roll, relative to the stacking energy of the ideal B-DNA conformation is given in Figure 8. It can be seen that the electrostatic contribution varies insignificantly with roll. The variation of the stacking energy with roll can, to an excellent approximation, be calculated using the LJ contribution alone. The LJ contribution is the sum of significant opposing attractive and repulsive components (not shown). The change in the solvent-DNA surface free energy is

L

0 -

-1

-

I I

# I

-

#

n

I

a

0

E

-

1

Ip

-2

-

-3

-

0

Y

Y

%

P Q)

c

w

! ! ! !

-4-

-5

! ! ! ! ! ! ! ! ! !

-

-

-6

-20

-30

-10

L I .

0

Roll

10

20

30

(O)

Figure 9. The dependence of electrostatic stacking energies on roll according to various FDPB. * * * * * * , dielectric models. Energies are relative to ideal B-DNA. 0 ,-, t," = 1.0, . . . . . . , c+, = rij. - - - - - - -Debye-Huckel. , m, - - - * - , Function of Sarai et al. 0,- - - - -, Hingerty function. A, -, Lavery function.

. . .

.

+,

155

*,

- - -

--

--

156

FRIEDMAN AND HONIG 4

0

I

I

I

20

4 0

60

Twist

f lI + * + A

80

100

(O)

Figure 10. The stacking energy of an AA step relative to ideal B-DNA as a function of twist. Symbols as in Figure 8.

insignificant compared to the change in the LJ interaction over the roll range studied. Similar results follow from repeating the calculation using the charge and LJ parameters of Poltev and co-workers. The electrostatic contribution will have a small ef-

fect on the statistical weights of the different roll contributions. The potential energy curves as functions of roll of TT, AT, TA, GG, GC, CG, AG, and GA steps were generated (not shown). All of these plots show the predominance of the LJ repulsive

L

-

20

-

10

-

n

z

\

(21

0 Y

v

) I P) L

0

c

W

’

O t -10

-30

I

I

1

I

I

1I

-20

-10

0

10

20

30

Propellor

(0)

Figure 11. The stacking energy of a n AA step relative to ideal B-DNA as a function of propellor. Symbols as in Figure 8.

ELECTROSTATIC CONTRIBUTION TO BASE STACKING

contribution. A study of the dependence of the stacking energy of AA, AT, TA, GG, GC, and CG steps on tilt (not shown) also indicates the predominance of LJ repulsive interactions that are not significantly compensated by interactions with the solvent. The electrostatic contribution to the roll energy according to various dielectric models is given in Figure 9. The Hingerty screening function mimics the FDPB results most closely. This tendency also holds for AT, TA, GG, GC, and CG (not shown), and so may be regarded to hold for other steps as well. The stacking energy of an AA step as a function of twist is given in Figure 10. Over the experimentally observed region of 25°-500,82,86-90the LJ component dominates. Several workers 10~18-20,24*91,92 have stated that the 36" twist of B-DNA is caused by a minimum (sometimes a shoulder) in the basestacking energy as a function of twist around 36" in which the electrostatic component plays a major role. Figure 10 shows that the electrostatic minimum is too shallow and broad in the vicinity of 36" to have a significant effect on the twist energy. Similar results pertain for AT, TA, GG, GC, and CG steps with the AMBER potential set (not shown), and for the AA step with the potential set of Poltev et al. DNA-solvent interactions reinforce LJ interactions over the experimentally observed twist range, so that the change in the LJ stacking interactions are not compensated by base-solvent interactions. The electrostatic energy of the same systems were studied as a function of dielectric model (not shown). It was found that the putative significant role of electrostatic energy in determining twist arises from the use of less accurate dielectric models, which exaggerate the significance of the electrostatic interaction. However, these models may have predictive value if regarded as effective potentials of mean force. As before, the Hingerty function proved to be an adequate way to calculate base-stacking energies. The stacking energy of an AA step as a function of propellor is given in Figure 11. The LJ component dominates the stacking energy. Both dispersive and repulsive contributions to the LJ stacking energy are significant. The change in stacking energy is destabilizing. Similar results pertain to GG, AT, TA, GC, and CG steps (not shown). These results contradict the contention of Dickerson et a1.61 and Calladine that propelloring occurs because it brings about enhanced attractive stacking interactions, and are in agreement with the conclusion of Shakked et al.13 and Spomer and Kyprg4that it does not. Inter-

157

actions with the solvent reinforces the repulsive interaction between the base pairs. The Hingerty function adequately models the electrostatic stacking energy as before.

SUMMARY AND CONCLUSIONS Electrostatic base-stacking interactions, at least for planar charge models, may be neglected to a good approximation. If greater accuracy is required, the Hingerty screening function should be used. A full FDPB treatment can be used for still greater accuracy, either for calculating energies or for calculating forces using the method of Sharp.93However, the magnitude of base-stacking energy usually will not justify the computational expense of this last method. LJ repulsive base-stacking interactions are a potentially significant determinant of DNA structure. The role of these interactions and base-solvent interactions will be the subjects of future papers.45 We thank our colleagues, Drs. B. Jayaram, Kim Sharp, Anthony Nicholls, Marilyn Gunner, and An-Suei Yang for software and helpful discussions. We thank Drs. Richard Lavery and Mihaly Mezei for furnishing software, and Gerald Manning for sending a preprint. This work was supported by the NIH (GM-41371) and NSF (DIR8720229).

REFERENCES 1. Dickerson, R. E. (1983) J. Mol. Biol. 166,419-441. 2. Marini, J. C., Levene, S. D., Crothers, D. M. & Englund, P. T. (1982) Proc. Natl. Acad. Sci. U S A 79, 7664-7668. 3. Patel, D. J. & Canuel, L. J. (1976) Proc. Natl. Acud. Sci. USA 73,3343-3347. 4. Simpson, R. T. & Kunzler, P. (1979) Nucleic Acids Res. 6,1387-1415. 5. Calladine, C. R. (1982) J. Mol. Biol. 161, 343-352. 6. Srinivasan, A. R., Torres, R., Clark, W. & Olson, W. K. (1987) J. Biomol. Struct. Dynam. 5,459-496. 7. Maroun, R. C. & Olson, W. K. (1988) Biopolymers 27, 561-584. 8. Maroun, R. C. & Olson, W. K. (1988) Biopolymers 27,585-603. 9. Sarai, A., Mazur, J., Nussinov, R. & Jernigan, R. L. (1988) in Structure & Expression Volume 3: D N A Bending and Curuature, Olson, W. K., Sarma, M. H., Sarma, R. H. & Sundaralingam, M., Eds., Adenine, Schenectady, NY, pp. 213-223. 10. Sarai, A., Mazur, J., Nussinov, R. & Jernigan, R. L. (1988) Biochemistry 27,8498-8502. 11. Sarai, A., Mazur, J., Nussinov, R. & Jernigan, R. L. (1989) Biochemistry 28,7842-7849.

158

FRIEDMAN AND HONIG

12. Mazur, J., Sarai, A. & Jernigan, R. L. (1989) Biopolymers 28, 1223-1233. 13. Haran, T. E., Berkovich-Yellin, Z. & Shakked, Z. (1984) J. Biomol. Struct. Dynam. 2,397-412. 14. Devoe, H. & Tinoco, I. (1962) J. Mol. Biol. 4, 500517. 15. Nash, H. A. & Bradley, D. F. (1965) Biopolymers 3, 261-273. 16. Pullman, B., Calverie, P. & Caillet, J. (1966) J . Theor. Biol. 12, 419-434. 17. Poltev, V. I. & Sukhorukov, B. I. (1967) Biofizika 12, 763-768. 18. Rein, R., Claverie, P. & Pollack, M. (1968) Znt. J. Quantum Chem. 2,129-144. 19. Polozov, R. V., Poltev, V. I. & Sukhorukov, B. I. (1975) J . Theor. Biol. 55, 491-503. 20. Ornstein, R. L., Rein, R., Breen, D. L. & MacElroy, R. D. (1978) Biopolymers 17, 2341-2360. 21. Kollman, P. A., Weiner, P. K. & Dearing, A. (1981) Biopolymers 20, 2583-2621. 22. Langlet, J., Claverie, P., Caron, F. & Boeuve, J. C. (1981) Znt. J . Quantum Chem. 19, 299-338. 23. Aida, M. & Nagata, C. (1986) Znt. J. Quantum Chem. 39, 1253-1261. 24. Bertran, J. (1972) J . Theor. Biol. 34, 353-361. 25. Poltev, V. I. & Sukhorukov, B. I. (1968)Biofizika 13, 941-949. 26. Gilson, M., Rashin, A., Fine, R. & Honig, B. (1985) J . Mol. Biol. 183, 503-516. 27. Zhurkin, V. B., Poltev, V. I. & Florent’ev, V. L. ( 1981) Molek. Biol. (Russ.) 14, 1116-1130. 28. Lavery, R., Sklenar, H., Zakrzewska, K. & Pullman, B. (1986) J. Biomol. Struct. Dynam. 3,989-1014. 29. Singh, U. C. & Kollman, P. A. (1984) J . Comp. Chem. 5 , 129-145. 30. Pearlman, D. A. & Kim, S. H. (1989) J . Mol. Biol. 211,171-187. 31. McCammon, J., Wolynes, P. G. & Karplus, M. (1979) Biochemistry 18,927-942. 32. Hingerty, B. E., Ritchie, R. H., Ferrel, T. L. & Turner, J. E. (1985) Biopolymers 24, 427-439. 33. Lavery, R. ( 1987) in Structure & Expression Volume 3: D N A Bending and Curvature, Olson, W. K., Sarma, M. H., Sarma, R. H. & Sundaralingam, M., Eds., Adenine, Schenectady, NY, pp. 191-211. 34. Debye, P. & Huckel, E. (1923) Physik. Z. 24, 185206. 35. von Kitzing, E. & Diekmann, S. (1987) Eur. Biophys. J . 15,13-26. 36. Klapper, I., Hagstrom, R., Fine, R., Sharp, K. & Honig, B. (1986) Proteins 1,47-59. 37. Gilson, M. K. & Honig, B. H. (1988) Proteins 3, 3252. 38. Gunner, M. R. & Honig, B. (1990) in Perspectives in Photosynthesis, Jortner, J. & Pullman, B., Eds., Kluwer Academic, Dodrecht, pp. 53-60. 39. Gunner, M. R. & Honig, B. (1991) Proc. Natl. Acud. Sci. USA, 88,9151-9155.

40. Rashin, A. A. & Honig, B. (1985) J. Phys. Chem. 89, 5588-5593. 41. Gilson, M. K. & Honig, B. (1988) Proteins 4, 7-18. 42. Jean-Charles, A., Nicholls, A., Sharp, K., Honig, B., Tempczyk, A., Hendrickson, T. & Still, C. (1990) J. Am. Chem. SOC.113,1454-1455. 43. Jayaram, B., Sharp, K. A. & Honig, B. (1989) Biopolymers 28,975-993. 44. Friedman, R. A. & Honig, B. (1989) in Book of Ab-

stracts: Sixth Conversation in Biomolecular Stereodynamics, Sarma, R. H., Ed., Institute of Biomolecular Stereodynamics, Albany, NY, pp. 95-96. 45. Friedman, R. A. & Honig, B., manuscript in preparation. 46. Drew, H. R. & Dickerson, R. E. (1981) J . Mol. Biol. 151, 535-556. 47. Westhof, E. (1988) Ann. Rev. Biophys. Chem. 17, 125-44. 48. Halle, B. & Piculell, L. (1982) J . Chem. Soc. Faraday Trans. 1 78,255-271. 49. Leyte, J. C. (1990) Makromol. Chem. Mucromol. Symp. 34,81-111. 50. Jayaram, B., Fine, R., Sharp, K. A. & Honig, B. (1989) J . Phys. Chem. 93,4320-4327. 51. Lavery, R. (1990) Jumna IZZe, Laboratorie de Bio-

chimie Theorique CNRS, Institut de Biologie PhysicoChimique, Paris. 52. Dickerson, R. E. (1989) J . Biomol. Struct. Dynam. 6, 627-634. 53. Arnott, S., Chandrasekaran, R., Birdsall, D. L., Leslie, A. G. W. & Ratliff, R. L. (1980) Nature 283, 743745. 54. Lavery, R. & Sklenar, H. (1989) J . Biomol. Struct. Dynam. 6,655-667. 55. Lavery, R. & Sklenar, H. (1988) J. Biomol. Struct. Dynam. 6,63-91. 56. Lavery, R. & Sklenar, H. (1989) Curves 2.1, Labor-

57.

58. 59. 60.

61.

atorie de Biochimie Theorique CNRS, Institut de Biologie Physico-Chimique, Paris. Lavery, R. ( 1990) Manual to Jumnu IZZe, Laboratorie de Biochimie Theorique CNRS, Institut de Biologie Physico-Chimique, Paris. Arnott, S. & Hukins, D. W. (1972) Biochem. Biophys. Res. Commun. 47,1504-1510. Mezei, M. (1986) Builddna 1.1, Department of Chemistry, Hunter College, New York. Drew, H. R., Wing, R. M., Takano, T., Broka, C., Tanaka, S., Itakura, K. & Dickerson, R. E. (1981) Proc. Natl. Acad. Sci. USA 78, 2179-2183. Dickerson, R. E. & Drew, H. (1981) J. Mol. Biol. 149,

761-786. 62. Drew, H. R. & Dickerson, R. E. (1981) Protein Data

Bank, 1BNA. 63. Biosym. ( 1990) ZnsightZZ 1.0, Biosym, San Diego, CA. 64. Sharp, K. A. & Nicholls, A. (1990) Delphi 3.0, Dept. of Biochemistry and Molecular Biophysics, Columbia

University, New York. 65. Gilson, M., Sharp, K. A. & Honig, B. (1988) J . Comp. Chem. 9,327-335.

ELECTROSTATIC CONTRIBUTION T O BASE STACKING

66. Weiner, S . J., Kollman, P. A., Nguyen, D. T. & Case, D. A. (1986) J. Comp. Chem. 7,230-252. 67. Bockris, J. 0. & Reddy, A. K. N. (1970) Modern Electrochemistry, Vol. 1, Plenum, New York, pp. 180-238. 68. Jayaram, B., Swaminathan, S., Beveridge, D. L., Sharp, K. & Honig, B. (1990) Macromolecules 23, 3156-3165. 69. Fenley, M. O., Manning, G. S. & Olson, W. K. ( 1990) Biopolymers 30, 1191-1203. 70. Bharadwaj, R., Friedman, R. A. & Honig, B. (1991) manuscript in preparation. 71. Cieplak, P. & Kollman, P. A. (1988) J. Am. Chem. SOC. 110, 3734-3739. 72. Dang, L. X. & Kollman, P. A. (1990) J. Am. Chem. SOC.112,503-507. 73. Poltev, V. I. & Shulyupina, N. V. (1986) J . Biomol. Struct. Dynam. 3, 739-763. 74. Maitland, G. C., Rigby, M., Smith, E. B. & Wakeham, W. A. ( 1987) Intermolecular Forces: Their Origin and Determination, Oxford University Press, Oxford, pp. 10-28. 75. Atkins, P. W. (1986) Physical Chemistry, Freeman, New York, pp. 576-584. 76. Gilson, M. & Honig, B. ( 1986)Biopolymers 25,20972119. 77. Nicholls, A., Sharp, K. A. & Honig, B. (1991) Proteins, 11,281-296. 78. Sharp, K. A., Nicholls, A., Fine, R. M. & Honig, B. (1991) Science 252,106-109. 79. Nicholls, A., & B. Honig (1991) GRASP 1.0, Department of Biochemistry and Molecular Biophysics, Columbia University, New York. 80. Sharp, K., Honig, B. & Harvey, S. (1989) Biochemistry 29,340-346. 81. Adamson, A. W. (1990) Physical Chemistry of Surfaces, Wiley, New York, p. 70.

159

82. Yoon, C., Prive, G. G., Goodsell, D. S. & Dickerson, R. E. (1988) Proc. Natl. Acad. Sci. U S A 85, 63326336. 83. Klug, A., Jack, A., Viswamitra, M. A., Kennard, O., Shakked, Z. & Seitz, T. A. (1979) J. Mol. Biol. 131, 669-680. 84. Quintana, J. R., Grzeskowiak, K., Yanagi, K. & Dickerson, R. E. (1991) J. Biomol. Struc. Dynam. 8, a171. 85. Claverie, P. (1968) J. Chim. Phys. 65,57-65. 86. Dickerson, R. E., Kopka, M. L. & Pjura, P. (1985) in Biological Macromolecules and Assemblies: Vol. 2Nucleic Acids and Interactive Proteins, Jurnak, F. & McPherson, A., Eds., Wiley, New York, pp. 38-126. 87. Frederick, C. A., Quigley, G. J., van der Marel, G., van Boom, J. H., Wang, A. H.-J. & Rich, A. (1988) J. Biol. Chem. 263,17872-17879. 88. Frederick, C. A., Williams, L. D., Ughetto, G., van der Marel, G. A., van Boom, J. H., Rich, A. & Wang, A. H.-J. (1990) Biochemistry 29, 2538-2549. 89. Heinemann, U. & Alings, C. ( 1989) J . Mol. Biol. 210, 369-381. 90. Nelson, H. C. M., Finch, J. T., Luisi, B. F. & R u g , A. (1987) Nature 330,221-226. 91. Rein, R. (1978) in Diutomics to Biopolymers, Pullman, B., Eds., Wiley, New York, pp. 308-362. 92. Fujita, H., Imammura, A. & Nagata, C. (1974) J. Theor. Biol. 45, 411-453. 93. Sharp, K. A., Nicholls, A., Friedman, R. & Honig, B. (1991) Biochemistry, 30,9686-9697. 94. Spomer, J. & Kypr, J. ( 1990) in Theoretical Chemistry & Molecular Biophysics, Volume I : DNA, Beveridge, D. L. & Lavery, R., Eds., Adenine, Schenectady, NY, pp. 271-284.

Received August 26, 1992 Accepted October 2, 1991