INTRODUCTION

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION CONCLUSIONS REFERENCES

Computer-aided modeling of DNA, RNA, and protein structures has become a partner with experiment in the elucidation of molecular structure-function relationships. Molecular mechanics and dynamics techniques are routinely used to generate feasible structures in solution, both in the refinement of high-resolution nuclear magnetic resonance (NMR) data and for de novo predictions of conformation and interactions. Advances in computer architecture and in force-field development support the increasing utility of these methods, with elaborate simulations of biomolecular dynamics (with full-scale solvation and explicit treatment of electrostatic interactions) now performed beyond the nanosecond range and free energies derived from these simulations in reasonable accord with experiment (Jayaram et al., 1999; Beveridge and McConnell, 2000; Kollman et al., 2000; Cheatham and Young, 2001; Yan et al., 2001). Nonetheless, the "multiple minimum" or "sampling" problem remains unsolved, precluding identification of both the lowest energy state and all energetically accessible conformations. Furthermore, simulations of molecular dynamics incorporating full-scale solvation are still too costly to sample the full conformational landscape.

Continuum electrostatic models address the issue of computational cost involved in explicit molecular calculations by greatly reducing the number of degrees of freedom of the solvent (reviewed in Cramer and Truhlar, 1999). There are a number of possible strategies for treating the electrostatic component of the solvation free energy, including numerical solution of the Poisson-Boltzmann equation (Sharp and Honig, 1990; Honig and Nicholls, 1995), implementation of the generalized Born approximation (Still et al., 1990; Cramer and Truhlar, 1992; Hawkins et al., 1995, 1996; Tsui and Case, 2001), and use of distance-dependent dielectric functions (reviewed in Mehler, 1996). These methods can be used in conjunction with other molecular data, such as the solvent-accessible surface area (Lee and Richards, 1971), to evaluate nonelectrostatic energy contributions (Hermann, 1972; Still et al., 1990; Cramer and Truhlar, 1992). Solution of the nonlinear Poisson-Boltzmann equation, although very accurate, is computationally expensive and not yet well suited for routine conformational searches (Baker et al., 2001). The generalized Born method, by contrast, is computationally efficient and benchmarked to reproduce the solvation free energies of small organic molecules (Best et al., 1997; Jayaram et al., 1998). Moreover, recent parameterizations and modifications of the generalized Born approximation appear to have corrected inaccuracies in the treatment of macromolecules (Srinivasan et al., 1999; Tsui and Case, 2000, 2001), although further tests will be of interest.

Distance-dependent dielectric functions, based on Lorentz-Debye-Sack theory (Sack, 1926, 1927; Debye, 1929; Lorentz, 1952), have a long history in the treatment of molecular interactions and have been widely used, because of their efficiency, to account implicitly for solvent effects in molecular mechanics calculations. According to these schemes, the dielectric constant between charged atoms increases with their distance of separation, tapering off at large interatomic displacements to the value of bulk solvent and thereby mimicking the expected diminution of electrostatic interactions as solvent molecules interpenetrate between the atom pairs. There is uncertainty, however, in the choice of a suitable functional form for the dielectric constant, particularly for the different environments of a macromolecule, such as the major and minor grooves of B-DNA. For example, a preferential build-up of counterions or co-ions around the double helix can change the dielectric medium and the distance-dependent interactions between charged atoms on the edges of the grooves (Young et al., 1998). The selection of a suitable dielectric function is thus key to the performance of a nucleic acid force field. The functional form introduced by Hingerty et al. (1985) is reasonably robust when compared against structural predictions based on more rigorous and expensive Poisson-Boltzmann calculations (Daggett and Kollman, 1990; Friedman and Honig, 1992, 1995).

A number of effective computational strategies for addressing the multiple minimum problem in the study of DNA and carcinogen-damaged DNA have been developed with the molecular mechanics program DUPLEX (Hingerty et al., 1989; Broyde and Hingerty, 1999). These include large-scale conformational searches for energetically favorable orientations of carcinogen with respect to DNA, progressive build-up of DNA structure from smaller mono- and oligonucleotide subunits, and restriction of nucleic acid conformational flexibility to the six exocyclic chain torsion angles and the puckering of the sugar rings. Bond lengths and bond angles are held fixed at standard equilibrium values and aromatic moieties are assigned planar conformations unless found otherwise by experiment. A sigmoidal, distance-dependent dielectric function (Hingerty et al., 1985) is used to treat solvation, and reduced partial charges are introduced on the pendant phosphate oxygens (Srinivasan and Olson, 1980) to mimic the effects of counterion condensation (Manning, 1978). Partial charges are taken from Ornstein and Rein (1979) with modifications noted in Hingerty et al. (1989). DUPLEX uses the linked atom algorithm (Scott and Scheraga, 1965) for coordinate generation from the bond lengths, bond angles, and dihedral angles, and the Powell algorithm (Powell, 1964) for energy minimization. The approach has proved useful both in prediction (Singh et al., 1991) and in refinement of DNA solution structures (Geacintov et al., 1997; Patel et al., 1998).

The NMR solution structures of DNA modified covalently at the C8 atom of guanine by the aromatic amine, 2-aminofluorene (AF) (Fig. 1 A), (Mao et al., 1998a,b) provide a useful starting point for reconsideration of the treatment of solvent and other aspects of the DUPLEX force field. Specifically, NMR studies of the AF-bound dodecamer duplex, d(CTCGG*CGCCATC)·d(GATGGCGCCGAG) (Sequence I, Fig. 1 B) with AF modification at G*, reveal a predominant, well-characterized conformer that places the AF in the B-DNA major groove while keeping all Watson-Crick base pairs (bp) intact. This conformer, present to the extent of ~90% of the population, is in equilibrium with another conformer that cannot be fully characterized with NMR distance restraints. Qualitative NMR data, however, indicate that the DNA in the second form adopts a structure with the AF inserted into the helix and the modified guanine displaced into the major groove. The latter conformer has been fully characterized in a different 11-bp duplex, d(CCATCG*CTACC)·d(GGTAGCGATGG) (Sequence II, Fig. 1 C) with AF modification again at G*. In this sequence context, the base-displaced, AF-intercalated conformer comprises ~70% of the conformational mixture. Moreover, the NMR data indicate that the base-displaced AF-intercalated states are essentially the same in the two sequence contexts. In both DNA adducts, there are two rapidly interconverting rotamers of comparable flexibility, with the C9-containing edge of the external, groove-bound AF directed toward either the 5'- or the 3'-side of the modified base and the C9-containing edge of the intercalated, base-displaced AF directed toward either the major groove or the minor groove of the double helix. The ~9:1 conformer ratio in Sequence I is indicative of a free energy difference of roughly 1.3 kcal/mol between the two states at 298 K, assuming equal conformational flexibilities. It should be noted that NMR structures within the bounds of the data represent ensemble averages rather than single points. This system serves as an excellent test case for the implicit treatment of solvent around DNA because the solvent exposure of carcinogen is so different in the two conformational states. We have thus evaluated and refined different distance-dependent dielectric functions for representative conformers of AF-bound DNA and further tested the capabilities of the expressions to reproduce various sequence-dependent structural features of unmodified B-DNA.

View larger version (24K): [in this window] [in a new window]   FIGURE 1   (A) Chemical structure of the AF-dG adduct (G*) showing the torsion angles (, ', ') that govern AF-DNA conformation. Angles are defined as follows: , O4'-C1'-N9-C4; ', N9-C8-N-C2; ': C8-N-C2-C1, where the C1, C2, and N are from 2-aminofluorene and the other atoms are from deoxyguanosine. (B, C) Sequences of high-resolution NMR structures, I (Mao et al.; 1998a) and II (Mao et al., 1998b), used to benchmark AF orientation in energy calculations involving the boxed 9-base pair sequences.

The right-handed, self-complementary B-DNA sequence, d(CGCGAATTCGCG)₂, known as the Dickerson-Drew dodecamer (Wing et al., 1980; Drew et al., 1981), has been crystallized under various conditions, and the many high-resolution crystal structures (Berger et al., 1998; Egli et al., 1998; Shui et al., 1998a,b; Tereshko et al., 1999a,b; Sines et al., 2000; Woods et al., 2000) now at hand in the Nucleic Acid Database (NDB; Berman et al., 1992) can be compared with DUPLEX-generated B-DNA structures based on different force fields. Among the most interesting aspects of the dodecamer structures are the subtle irregularities of the double-helix associated with the sequential positions of purine and pyrimidine bases along the complementary strands. This fine-scale information is not found in the idealized, regularly repeating B-DNA helical structure that can be derived from X-ray fiber diffraction measurements (Arnott et al., 1976). The sequence-dependent variation of base-pair geometry, backbone torsion angles, and groove dimensions computed on the basis of the DUPLEX force field is compared below against the corresponding conformational features observed in high-resolution Dickerson-Drew dodecamer structures. To improve the computational predictions, we reconsider the treatment of deoxyribose pseudorotation and the long-accepted strategy of minimizing the Lennard-Jones potential energy at a distance 0.2 Å greater than the van der Waals' separation (Brant et al., 1967; Olson and Flory, 1972; Srinivasan and Olson, 1980), in addition to modifying the sigmoidal, distance-dependent dielectric function.

Thus, we use the relative energies of AF-modified DNA duplexes in equilibrium between AF-external and AF-intercalated base-displaced conformers, and the sequence-dependent structural information of the Dickerson-Drew dodecamer as criteria for assessing the viability of DNA computational predictions based on the efficient and inexpensive treatment of solvation and counterions within the DUPLEX force field. A dielectric functional form that tapers more steeply than earlier functions to the dielectric constant of bulk water is found to reproduce fairly well the expected energy difference between the two conformers of AF-bound DNA. This function, together with a modified deoxyribose pseudorotation profile with reduced energy barriers between C2'- and C3'-endo minima, and a Lennard-Jones potential energy minimum at a distance 0.38 Å greater than the normal van der Waals' separation, accounts satisfactorily for the sequence-dependent structural features of Dickerson-Drew dodecamer crystals. The improved force field shows even better correspondence with the dodecamer structure of the same sequence deduced from NMR measurements in an aqueous dilute liquid crystalline phase (Tjandra et al., 2000). The robust performance of the revised force field, in conjunction with its efficiency through implicit treatment of solvent and counterions, provides a valuable tool for elucidating structure-function relationships of DNA, including those of chains modified by carcinogens and other ligands.

	METHODS

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION CONCLUSIONS REFERENCES

Computer simulation of the AF-modified dodecamer, Sequence I (Fig. 1 B), was performed using the NMR solution structure of the major AF-external conformer as a starting model (Mao et al., 1998a) (see Table 1). Only one of the two rotamers, the form with the C9-containing edge of AF directed toward the 5' side of the modified base (' = 160°, ' = 38°), was considered. The 12-mer was truncated to a 9-mer so that the modified G occupied the central position of the double helix (boxed segment in Fig. 1 B). Past experience with energy minimization of carcinogen-modified DNA (Shapiro et al., 1998) indicates that 9-mers are long enough to include all DNA-carcinogen interactions and concomitantly avoid, as much as possible, DNA end effects unrelated to adduct conformation. The secondary base-displaced, AF-intercalated conformer of the same sequence was modeled using the principal conformation detected in high-resolution NMR studies of the AF-modified Sequence II 11-mer (Mao et al., 1998b) as a starting state. The (' = 148°, ' = 138°) rotamer with the C9-containing edge of AF directed toward the minor groove was utilized. The 11-mer was then shortened by one residue at either end to create a 9-mer with a centrally modified G, and the sequence adjusted to match that of Sequence I (see Fig. 1, B and C). Subsequent minimization of the two 9-mers yielded structures very similar to the respective NMR parent structures with the key , ', ' torsion angles (Fig. 1 A) which govern the AF-DNA conformations departing no more than 7° from the starting states (Table 1).

We chose to model the minor (10%) AF-intercalated base-displaced structure of Sequence I rather than the secondary (30%) AF-external conformation of Sequence II because the latter form of the Sequence II complex is not as well defined by NMR distance restraints as the former arrangement of the Sequence I adduct. In addition, Sequence I is biologically more interesting in that it contains the NarI mutagenic hotspot for aromatic amine carcinogens (reviewed in Hoffmann and Fuchs, 1997). The (, ', ') rotamer choices noted above correspond to the lower energy (and better characterized) states found in the refinement of the NMR structures (Mao et al., 1998a,b). In the first tests of various dielectric functions, the energies of AF-external and AF-intercalated, base-displaced 9-mer duplexes were computed without further energy minimization, to avoid additional structural variations associated with the multiple minimum problem. Specifically, the comparison of identical structures for each tested dielectric function provides direct information on how the steepness of the dielectric function affects the computed conformational energy difference. This strategy also closely preserves the structures of the NMR-characterized states.

Crystal and DUPLEX-generated Dickerson-Drew dodecamer structures were analyzed at the local base-pair level using the RNA (Run Nucleic Acid) software package (Babcock et al., 1993) and at the level of torsion angles and groove geometry with the routines of Curves 5.2 (Lavery and Sklenar, 1988, 1989). Thirteen crystal structures (Table 2) from the Nucleic Acid Database (Berman et al., 1992) with resolution ranging from 0.9 to 1.9 Å and without any mismatches, flipped-out bases, bound drug molecules, or bulky modifiers were chosen to benchmark the DUPLEX-generated structures. The selected crystals fall into three groups. The first group includes three independent refinements of a structure determined in the 1980s (Drew et al., 1981; Holbrook et al., 1985; Westhof, 1987) with a resolution of ~1.9 Å. The second group contains four structures from the Williams laboratory solved in 1998-2000 (Shui et al., 1998a,b; Sines et al., 2000; Woods et al., 2000), with resolution in the range of 1.2-1.8 Å. The third group is composed of six structures from the Egli laboratory solved in 1998-1999 (Berger et al., 1998; Egli et al., 1998; Tereshko et al., 1999a,b), with resolution in the range of 0.9-1.6 Å, and with one or two 2'-deoxy-2'-fluoroarabino-thymidines in place of thymine in each strand. The chemically modified residues in the latter molecules do not affect the overall duplex structure, but do enhance the thermodynamic stability of B-form DNA (Rosenberg et al., 1993; Ikeda et al., 1998). The unprecedented high resolution of this set of crystal structures makes them particularly attractive as benchmarks for the DUPLEX force field.

In attempting to reproduce the sequence-dependent structural features observed in dodecamer crystals, we also considered four different modifications of the deoxyribose pseudorotation profile: 1) the current list of variablespuckering amplitude _m, internal bond angles, and energies as a function of the pseudorotation phase angle Pused in DUPLEX (Sasisekharan, 1973; Hingerty and Broyde, 1982); 2) the current input from DUPLEX with energies reduced by half; 3) the computational treatment originally offered by Sato (Sasisekharan, 1973; Sato, 1983); and 4) the Sato treatment with energies reduced by half. Profiles 1 and 4 are shown in Fig. 2 (gray and red lines, respectively). It should be noted that the bases are replaced by a hydrogen in the Sato treatment, so that the conventional high barrier in the western region, i.e., P  270°, (Saenger, 1984) of the pseudorotation cycle is not seen.

View larger version (32K): [in this window] [in a new window]   FIGURE 2   Deoxyribose pseudorotation profiles of puckering amplitude m, internal bond angles, and conformational energy of the sugar ring as a function of the pseudorotation phase angle, P. Note that the base is replaced by a hydrogen in the profile (Sasisekharan, 1973), and a gauche energy term (Olson, 1982) is included to account for the known preference of C2'-endo (P  180°) over C3'-endo (P  0°) puckered forms. Color coding: gray, original DUPLEX pseudorotation treatment, profile 1 (Hingerty and Broyde, 1982); red, current improved treatment, profile 4. Shaded areas correspond to the observed averages and standard deviations of geometric parameters in high-resolution crystal structures of mononucleosides and mononucleotides (Gelbin et al., 1996).

In the final stage of force-field improvements, we varied the location of the Lennard-Jones potential energy minimum over a range of distances 0.2-0.4 Å greater than the sum of the van der Waals' radii, first adding 0.2, 0.3, or 0.4 Å to the normal van der Waals' separation and then investigating finer (0.02-Å) increments of added distance between 0.3-0.4 Å.

	RESULTS

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION CONCLUSIONS REFERENCES

Search for a dielectric functional form that reasonably reproduces the energetics of AF-modified DNA conformers

Our approach to the enhancement of the DUPLEX force field entailed an iterative procedure. We first tested the functional form of the distance-dependent dielectric constant with a pair of AF-modified DNA 9-mers derived from NMR solution studies (see Methods).

The original distance-dependent dielectric function of Hingerty et al. (1985) was recast to the analytical form

(1)

Here (r) is the dielectric constant when the two atoms of interest are separated by a distance r, whereas c and e are variable quantities that are used to adjust the expression to virtually any sigmoidal form. Values of 2.5 for c and 2.0 for e produce the function used previously (Hingerty et al., 1985). The representative functional forms in Fig. 3 illustrate how (r) becomes steeper as c decreases and/or e increases.

View larger version (27K): [in this window] [in a new window]   FIGURE 3   Representative functional forms of the distance-dependent dielectric constant, (r), produced with selected values of the variables c and e in Eq. 1. The steepness of the function increases as c decreases and/or e increases. The original dielectric functional form (Hingerty et al., 1985) used by DUPLEX is shown in gray (Curve C); the functional form fitted to molecular dynamics simulations of Young et al. (1998) is shown in green (Curve B). The steep function used in current work is shown in red (Curve A). Curves D-F illustrate the change in (r) as c and e change.

The observed ~9:1 ratio between AF-external and base-displaced, AF-intercalated conformers of Sequence I corresponds to an expected energy difference of roughly 1.3 kcal/mol at 298 K, assuming comparable flexibility in the two conformational forms. The energies of representative NMR-derived states were computed with DUPLEX without minimization over the ranges of c and e shown in Table 3. The computed energy differences between the two conformers, in kcal/mol, are reported in the table for selected values of the two constants. We note that, as the steepness of the function increases (low c, high e), the energy difference between the two conformers becomes smaller. The original Debye treatment (Hingerty et al., 1985) with c = 2.5 and e = 2.0 (Fig. 3, curve C) gives an energy difference of 5.3 kcal/mol, over-favoring the AF-external conformation compared to the AF-intercalated, base-displaced state. We chose c = 0.5, e = 6.0 for further study because that combination most closely approximates the expected energy difference of ~1.3 kcal/mol (Table 3), and further steepening in the functional form (c = 0.5, e = 8.0) does not improve the energetics. As shown below, this function performs quite well in reproducing the sequence-dependent, base-pair and backbone conformational parameters observed in crystals of the Dickerson-Drew dodecamer.

View this table: [in this window] [in a new window]   TABLE 3   Energy differences, E (kcal/mol), between AF-external and AF-intercalated base-displaced DNA conformers for selected values of c and e in the distance-dependent dielectric function Eq. 1

Optimization of the deoxyribose energy and bond angle profiles to account better for the sequence-dependent structural parameters in crystal structures of the Dickerson-Drew dodecamer

In the next stage, we tested the revised force field with the new dielectric function by evaluating its performance in reproducing sequence-dependent base and backbone parameters observed in crystals of the Dickerson-Drew dodecamer. We found that the sequence-dependent pattern of most base-pair step parameters was reproduced reasonably well (data not shown). The DUPLEX-generated structure, however, was over-twisted, with an average dimer twist of 38.7 ± 1.6°, compared to an average value of 35.5 ± 4.4° in the crystal structures. In addition, the base-pair steps showed excessive negative roll into the minor groove, so that the minor groove width was much narrower on average (2.6 ± 1.5 Å) than that in the crystal structures (4.2 ± 1.2 Å), although not as narrow as with the original force field (Table 4). Further adjustments of the dielectric function did not improve these structural anomalies. Knowing that high twist is linked to a high phase angle of pseudorotation and negative roll (Olson and Zhurkin, 2000), we revisited the treatment of the pseudorotation phase angle in DUPLEX. The program previously used a modification of the deoxyribose pseudorotation potential (Hingerty and Broyde, 1982) devised by Sato (1983) and Sasisekharan (1973), in combination with a gauche energy term that is needed to reproduce the known preference of C2'-endo over C3'-endo puckered forms in deoxyribose (Olson, 1982). A test of the four deoxyribose pseudorotation profiles described in the Methods section showed that neither Sato's original deoxyribose treatment (Case 3 in the preceding section) nor the profile produced by halving the energy in the DUPLEX pseudorotation profile (Case 2 above) improved the overtwisting of the helix and the narrowness of the minor groove. A combination of the two (Case 4 above), however, satisfactorily produces a B-like minor groove and a reasonable sequence-dependent twisting pattern, without adverse effects on other measures of base-pair geometry, or the torsional parameters (data not shown). The change in barrier height follows from arguments developed to account for the undertwisting of DNA and the intermediate sugar puckering phase angles found in many computer simulations of B-DNA structures (Olson and Zhurkin, 2000). In the present work, the barrier is reduced to lower the twist, whereas the barrier must be raised in simulations based on other force fields, in which the computed structures are undertwisted.

Relocating the Lennard-Jones potential energy minimum to a distance 0.38 Å greater than the van der Waals' separation reasonably mimics the base-pair rise observed in crystal structures

At this stage, it was noted that the magnitude of base-pair rise and helical rise, nonsequence-dependent but important parameters for characterizing double helical structures (Sponer and Kypr, 1993; Hunter and Lu, 1997), were consistently smaller in the computer-generated dodecamer than the average values in the crystal structures. Consequently, the simulated dodecamer was visibly shorter than its crystal counterparts, with an average helical rise of 3.0 ± 0.0 Å compared to a value of 3.3 ± 0.1 Å in the crystal structures. In fact, the helical rise had diminished, following improvement of the dielectric function, from a value of 3.2 ± 0.2 Å obtained with the original force field (Table 4). The reduction in electrostatic interactions with the new dielectric function collapsed the DNA locally, allowing neighboring base pairs to approach more closely than in the original force field and in crystals. Changes in either the dielectric function or the deoxyribose pseudorotation profile did not alter either the local base pair rise or the helical rise significantly. Knowing that a large number of van der Waals' attractions can bring interacting moieties into closer contact than isolated atom pairs (Brant et al., 1967; Olson and Srinivasan, 1990), we investigated improving the rise values by relocating the Lennard-Jones potential energy minimum. Specifically, we reevaluated whether the Lennard-Jones energy minimum should remain at the distance 0.2 Å greater than the van der Waals' separation used in the first empirical energy studies of polypeptide (Brant et al., 1967) and polynucleotide (Olson and Flory, 1972) random coils, or whether it should be increased to a value more appropriate to double helical DNA. Notably, the effective van der Waals' radii of nucleotide atoms in the Poltev force field (Zhurkin et al., 1981), which was parameterized on the basis of the interactions of aromatic molecules that mimic the stacking of nucleic acid bases, are larger by as much as 0.2 Å than those in force fields originally developed to treat alkanes and peptides (Olson and Srinivasan, 1990). As detailed in the Methods section, we performed a series of tests in which the minimum in the Lennard-Jones terms was relocated to distances in the range of 0.2-0.4 Å greater than the sum of the van der Waals' radii. Preliminary results showed that changes in the location of the energy minimum over this range did not have any adverse effects on the sequence-dependent structural features of the double helix, but did, as expected, systematically increase the local rise between neighboring base-pair planes. An adjusted Lennard-Jones potential energy minimum at a distance 0.38 Å greater than the van der Waals' separation produced average local base-pair and helical rise values of 3.2 ± 0.1 Å, in remarkably good agreement with crystallographic observations (3.3 ± 0.1 Å). Additional displacement of the Lennard-Jones potential energy minimum to a distance 0.4 Å greater than the sum of van der Waals' radii, however, bent the computed DNA structure severely, indicating that further adjustment along these lines was not feasible.

Figures 4-6 and A1 (Appendix) compare selected sequence-dependent conformational parameters of the bases and chain backbones in the DNA dodecamer structures predicted on the basis of the original and fully improved DUPLEX force fields with the corresponding values in high-resolution crystal (Table 2) and NMR (Tjandra et al., 2000) structures. The sequence-dependent features of the revised model more closely mimic the experimentally observed trends than those based on the original potential. Table 4 shows similar improvements in the mean values and standard deviations of complementary base-pair orientation angles, base-pair step parameters, local helical angles and distances, chain torsion angles, deoxyribose pseudorotation phase angle, and major/minor groove widths/depths in the computer-generated dodecamers compared to the crystal data. The stereo images in Fig. 7 show the approximate global similarity of the structure based on the improved force field, the original Dickerson-Drew dodecamer (NDB_ID: BDL001; Drew et al., 1981), and the dodecamer structure in solution (PDB_ID: 1DUF, model 1, Tjandra et al., 2000). The original crystal structure is virtually indistinguishable at this level from the higher-resolution Dickerson-Drew dodecamer structures listed in Table 2 (with heavy-atom root-mean-square differences less than 0.5 Å from BDL001).

View larger version (28K): [in this window] [in a new window]   FIGURE 4   Variation of selected sequence-dependent base-pair parameters of high resolution Dickerson-Drew dodecamer crystals (Table 2), the dodecamer structure in an aqueous dilute liquid crystalline phase (PDB_ID: 1DUF, model 1) (Tjandra et al., 2000, and the DUPLEX-generated B-DNA structures of the same sequence, using original and improved force fields. Values calculated using RNA (Babcock et al., 1993). Color coding: blue circles, crystals; green inverted triangles, 1DUF; red squares, improved force field (FF); gray triangles, original force field. Crystal data are averages at corresponding sites in different structures, with the brackets denoting standard deviations.

View larger version (24K): [in this window] [in a new window]   FIGURE 5   Variation of groove widths and depths of high-resolution Dickerson-Drew dodecamer crystals (Table 2), the dodecamer structure in an aqueous dilute liquid crystalline phase (PDB_ID: 1DUF, model 1) (Tjandra et al., 2000), and the DUPLEX-generated B-DNA structures of the same sequence, using the original and improved force fields. Values calculated using Curves 5.2 (Lavery and Sklenar, 1988, 1989). Color coding and symbols are the same as in Fig. 4.

View larger version (43K): [in this window] [in a new window]   FIGURE 6   Variation of torsion angles and deoxyribose pseudorotation phase angle P of high-resolution Dickerson-Drew dodecamer crystals (Table 2), and the DUPLEX-generated B-DNA structures of the same sequence, using the original and improved force fields. Values calculated using Curves 5.2 (Lavery and Sklenar, 1988, 1989). Color coding and symbols are the same as in Fig. 4. Residue numbers reported in Curves 5.2 were shifted +1 for  and , in accordance with the IUPAC convention (Saenger, 1984).

View larger version (33K): [in this window] [in a new window]   FIGURE 7   Stereo views of the original Dickerson-Drew dodecamer crystal (NDB_ID: BDL001; top) (Drew et al., 1981), the high-resolution NMR solution structure of the same dodecamer in an aqueous dilute liquid crystalline phase (PDB_ID: 1DUF, model 1; middle) (Tjandra et al., 2000), and the DUPLEX generated B-DNA of the identical sequence obtained with the fully revised force field (bottom).

We returned to the AF-modified DNA conformer mixture as an additional check of the revised force field. Reminimization of the representative structures yielded an energy difference of 1.8 kcal/mol in favor of the AF-intercalated, base-displaced form, compared to the experimentally expected ~1.3 kcal/mol preference for the AF-external form.

We then recreated Table 3, using the preceding two newly energy-minimized states, to complete the optimization cycle. The results in Table 5 show that energies within ±3 kcal/mol from the theoretical target of E = 1.3 kcal/mol favoring the AF-external conformer are reproduced by functions within the shaded area (between Curves A and D of Fig. 3). This observation prompted us to compare the new steep dielectric function (Fig. 3, Curve A) with the dielectric expression proposed for DNA by Young et al. (1998), (Fig. 3, Curve B), and the function originally used in DUPLEX (Fig. 3, Curve C), represented, respectively, in Table 5 with E of 1.5 kcal/mol (c = 0.5, e = 1.0) and 2.3 kcal/mol (c = 2.5, e = 2.0), in reproducing sequence-dependent features of the Dickerson-Drew dodecamer. In these comparisons, the improved pseudorotation profile and relocated Lennard-Jones potential energy minimum were used. Figs. A2 and A3 (Appendix) show that the new steep dielectric function performs as well as or slightly better than the other two functions, in accounting for both local base-pair geometry and groove dimensions of the crystal structures.

View this table: [in this window] [in a new window]   TABLE 5   Energy differences, E (kcal/mol), between AF-external and AF-intercalated base-displaced DNA conformers for selected values of c and e in the distance-dependent dielectric function Eq. 1 with improved deoxyribose pseudorotation profile and relocated Lennard-Jones potential energy minimum

Finally, the conformational features of the dodecamer structure generated with the fully revised force field are in excellent agreement with the recent NMR solution structure of the Dickerson-Drew dodecamer in an aqueous dilute liquid crystalline phase (Tjandra et al., 2000), using the lowest energy and best representative conformer in the ensemble of five reported structures (PDB_ID: 1DUF; Berman et al., 2000). The results in Figs. 4-5 and A1 (Appendix) show that the conformational predictions of the revised force field match the details of the solution structure even better than those of the crystals.

	DISCUSSION

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION CONCLUSIONS REFERENCES

Dielectric function

Electrostatic interactions are among the most important forces in determining structural and functional features of proteins and nucleic acids. The calculation of the electrostatic contribution to the potential energy is dependent on many factors, including the assigned magnitude of partial charges, the behavior of solvent and counterions, the accessibility of solvent to solute, the surface hydrophobicity, and the polarizability of solute. In all-atom models of biopolymers surrounded by explicit waters, the dielectric constant is assigned a value of unity and the electrostatic interactions of all possible atom pairs are considered. This exact and computationally very expensive approach requires elaborate treatments, such as molecular dynamics or Monte Carlo simulations, to equilibrate the solvent continuously. Early solvent-implicit treatments of nucleic acids performed with limited computer resources, by contrast, assumed a dielectric constant of 2-4 to mimic the organic medium of the sugars and bases between successive negatively charged phosphates (Ramachandran and Srinivasan, 1970; Olson and Flory, 1972). Such treatments, although highly effective for homogeneous macromolecules such as randomly coiling chains composed of identical chemical repeating units (Flory, 1953), are very rough approximations for ordered biological macromolecules such as nucleic acid helices or globular proteins, which, respectively, contain grooves or clefts that can partially sequester counterions and solvent molecules (Jayaram et al., 1989, 1990). The distance-dependent treatment of the dielectric constant, where (r) = r (Brooks et al., 1983), is an improvement over the preceding model, but still does not account accurately for the screening of solvent between charged atoms (Jayaram et al., 1989).

As a further improvement of implicit modeling, Hingerty et al. (1985) offered a modification of Debye's original sigmoidal, distance-dependent dielectric function (Debye, 1929). This analytical expression extended the dependence of the dielectric constant on the distance r between charges to less than 3 Å, yielding a value of  = 1 at r = 0, and  = 78.3, the dielectric constant of water, as r approached 20 Å (Fig. 3, Curve C). Ramstein and Lavery (1988) reformulated this function by introducing adjustable parameters to specify the plateau value of the dielectric term at large distances and the location of its half value. Mazur and Jernigan (1991) subsequently presented an adaptable sigmoidal function that can simulate the functions of both Hingerty et al. (1985) and Ramstein and Lavery (1988), as well as other functional forms that mimic the salt-dependent B- to A-form transition of neighboring base pairs in a poly dG· poly dC duplex. Here we propose a further modification of the earlier function (Hingerty et al., 1985), with two empirical variables so that the function can adopt virtually any sigmoidal form. The form selected here on the basis of the equilibria in solution of AF-modified DNA duplexes derived from high-resolution NMR measurements, reaches solvent saturation rather fast, at an interatomic distance of ~3 Å. The electrostatic field is thus strongly damped and close to the dielectric screening of bulk water at a distance corresponding to just one solvation layer away from the charged DNA skeleton. This functional form is much steeper than our earlier dielectric expression, but is closer to functions suggested by others (Grahame, 1950; Takashima and Schwan, 1965; Pennock and Schwan, 1969; Harvey and Hoekstra, 1972; Young et al., 1998; Hassan et al., 2000), including the functional form for B-DNA proposed by Mazur and Jernigan (1991).

We find that steep dielectric functions systematically disfavor the AF-external conformer relative to the AF-intercalated, base-displaced one (Table 5, negative entries). This preference appears to stem from an overdampening of electrostatic interactions between phosphate groups across the minor groove in the intercalated conformer, the width of which is widened relative to the external conformer; P···P distances in the intercalated and external conformers are respectively: P5 ··· P17, 13.2 vs. 12.9 Å; P6 ··· P16, 15.4 vs. 13.3 Å; P7 ··· P15, 16.4 vs. 13.0 Å, where the numbers on the phosphorus atoms refer to the residues in Fig. 1 B. Specific distinctions of this type between different conformers of carcinogen-modified DNAs highlight the fact that the dielectric environment differs in various regions of DNA, both modified and unmodified, as previously noted by Young et al. (1998). Accordingly, benchmark structures can offer only approximate energetic criteria.

Our new steep dielectric function is remarkably similar to the dielectric behavior deduced by Young et al. (1998) from the computed positions of explicit waters and Na⁺ counterions around a B-form DNA oligonucleotide duplex during a 14-ns molecular dynamics trajectory, using an adaptation of Kirkwood dielectric theory (Kirkwood, 1939). The average dielectric constant ( = 59) at a distance of 2.8 Å from the DNA surfacereflecting the different dielectric behavior of interactions involving major groove ( = 53), minor groove ( = 51), and backbone ( = 66) atoms in the simulationsis thus somewhat smaller than the value of 78.2 in our revised function at the same distance. Currently available experimental information on the dielectric constant of solvent in the DNA groovesa value of  = 55 in the major groove according to fluorescence studies of dansylated-DNA-drug complexes (Barawkar and Ganesh, 1995) and a value of  = 20 in the minor groove based on the fluorescence properties of bound bisbenzimide molecules (Jin and Breslauer, 1988)suggests, however, that the average dielectric constant at 2.8 Å may be even lower than the values deduced by Young et al. The empirical Hingerty-Lavery sigmoidal function to which the dielectric behavior of the solvent was fitted by Young et al. (1998) is compared with our new steep dielectric function in Fig. 3 (Curves B and A, respectively). The set of dielectric functions in water, developed by Mehler and associates for peptides and proteins (Mehler, 1996; Mehler and Guarnieri, 1999; Hassan et al., 2000), reaches a bulk value of 78.3 at ~5 Å, suggesting further similarity in the distance-dependent dielectric activity of proteins and nucleic acids.

The fluctuations in the NMR solution structures of AF-modified duplexes (Mao et al., 1998a,b) coupled with the multiple minimum problem make energetic comparison of the two representative conformers used in this study only an approximate benchmark for force-field optimization. The dielectric function presented here, however, can be easily adjusted to fit different contexts as other reliable benchmarks arise. It appears that functional forms like the present expression or that of Young et al. (1998) will be applicable in other cases.

As pointed out by a reviewer, our proposed steep dielectric function resembles a two-step dielectric model implemented in the Delphi version (Nicholls et al., 1990) of the finite difference adaptation of Poisson-Boltzmann theory (Honig and Nicholls, 1995). To explore whether such a two-step dielectric might be useful instead of the more costly steep function, we recomputed the energy difference between the two conformers used to generate Table 5, using  = 4 if r < 2.71 Å and  = 78.3 if r > 2.71 Å. This produced an energy difference of 1.7 kcal/mol, very close to the expected value of 1.3 kcal/mol, suggesting that such a two-step dielectric might be worthwhile. Of course, considerable further testing of this approach would be needed. Systematic comparison with finite difference adaptations of Poisson-Boltzmann theory would be interesting future work.

Deoxyribose pseudorotation profile

The DUPLEX package previously used a modification of the Sato deoxyribose pseudorotation profile (Sasisekharan, 1973) where the puckering amplitude _m, and bond angles were shifted with respect to the phase angle P of pseudorotation (Fig. 2, gray line) (Hingerty and Broyde, 1982). This modification was made primarily to match experimental and theoretical data available (Olson, 1982) when the force field was first constructed. The current work, however, shows that this DUPLEX-generated Dickerson-Drew dodecamer was overtwisted with an extremely narrow minor groove, if this pseudorotation treatment is used. These deficiencies have now been remedied with the deoxyribose pseudorotation profile shown in Fig. 2 (red line), which corresponds to the series of structures originally proposed by Sato (1983) but with their energies halved. Local energy minima occur at P = 0° (C3'-endo) and P = 180° (C2'-endo) with a slight favoring of the C2'-endo region by ~0.6 kcal/mol imposed by the gauche term (Olson, 1982). The absence of the well-known high-energy barrier through the western region (P  270°) of the pseudorotation cycle (Saenger, 1984) reflects the omission of the base in the energetic treatment by Sato (see Methods). The steric effects, which give rise to this barrier, are accounted for explicitly by base-sugar interactions in the DUPLEX force field. Notably, the internal bond angles of the deoxyribose ring in the improved profile agree quite well with values (Fig. 2, shaded boxes) found in ultra-high-resolution crystal structures of mononucleosides and mononucleotides (Gelbin et al., 1996). The correct treatment of the deoxyribose is key not only to simulating B-DNA structures with proper twist (Fig. 4), as suggested by Olson and Zhurkin (2000), but also to reproducing proper groove widths (Fig. 5) when implicit solvent and counterions are used.

Lennard-Jones potential energy minimum location

It has been known for many years, in computing Lennard-Jones interactions, that the individual van der Waals' radii need to be adjusted to counter lattice compression effects (Brant et al., 1967; Olson and Flory, 1972; Srinivasan and Olson, 1980; Olson and Srinivasan, 1990). Without such adjustment, the interaction potential between neighboring base pairs is minimized at stacking distances up to 0.6 Å less than the normal B-DNA base-pair separation of 3.4 Å (Olson, 1978; Yoon et al., 1984). Minimizing the Lennard-Jones potential energy at a distance 0.2 Å greater than the normal van der Waals' separation, following the early studies of Brant et al. (1967) and Olson and Flory (1972), has become a standard practice. The low rise between neighboring base pairs in the Dickerson-Drew dodecamer structure generated with the new steep dielectric function and modified deoxyribose pseudorotation profile was largely corrected by relocating the Lennard-Jones potential energy minimum to a distance 0.38 Å greater than the sum of the van der Waals' radii.

The improved van der Waals' term closely resembles the expressions in other commonly used nucleic acid force fields, including one by Poltev (Zhurkin et al., 1981), AMBER (Cornell et al., 1995), and CHARMM (Foloppe and MacKerell, 2000; Mackerell and Banavali, 2000), with the exception of the depth of the energy well of the phosphorus-phosphorus Lennard-Jones interactions. Table A1 (Appendix) compares values of r_min, the distance between an atom pair at the potential energy minimum, and V_min, the well-depth in our improved force field, with the corresponding parameters of the AMBER, CHARMM, and Poltev force fields. The value of V_min for phosphorus-phosphorus (P ··· P) interactions in our improved force field is 1.06 kcal/mol, but only 0.58 to 0.20 in the other force fields. Modification of the depth of the P···P well to values comparable to those used in other force fields (0.30), however, has only a minimal effect on the conformational parameters of the dodecamer. Figures A4 and A5 (Appendix) show these results.

Solution versus crystal structures

A key issue in testing a force field involves selection of the appropriate benchmark. Although crystal structures represent the largest body of structural data with the least ambiguity, the question of just how these observations relate to solution structures is an actively debated question (Harvey et al., 1995; Dickerson et al., 1996; Shatsky-Schwartz et al., 1997; Beveridge and McConnell, 2000). Notably, the many crystal structures of the Dickerson-Drew dodecamer are very similar to one another, despite ongoing controversies concerning issues, such as the number of layers of water in the hydration spine (Shui et al., 1998b) and the types of ions that play crucial coordinating roles (Tereshko et al., 1999b; McConnell and Beveridge, 2000). The subtle interplay of different counterions, crystallization conditions, packing forces, and their effects on structural nuances remains to be unraveled. Gross end effects, such as the rupture of terminal base pairs, which then pack with neighboring helices (Liu et al., 1998), are thought to account for some of the irregularities in conformational parameters of the dodecamer crystals in Figs. 4 and 6. Notable sequence-dependence in certain base-pair step parameters, especially roll and twist (Fig. 4), has been shown to stem from packing forces (Subirana and Faria, 1997).

The present force-field treatment, with implicit solvent and counterions, is directed toward simulating reasonable solution structures. The collective crystal data are simply used as a benchmark of the sequence-dependent, conformational structure of DNA rather than as measures of precise atomic arrangements in the crystal lattice. This approach, with the goal of reproducing sequence-dependent trends within the range of observed data, transcends specific issues concerning counterions and the hydration spine. In this connection, it should be noted that, although the improved force field accounts reasonably for the sequence-dependence of widths and depths of the major and minor grooves, the magnitudes of the parameters deviate somewhat from crystal values (Fig. 5, red squares and blue circles). Moreover, computations using the less steep, distance-dependent dielectric function of Young et al. (1998) (Fig. 3, Curve B) in combination with the improved pseudorotation profile and relocated Lennard-Jones potential energy minimum perform similarly (Fig. A3, magenta squares; Appendix). Interestingly, a 15-ns molecular dynamics study of the Dickerson-Drew dodecamer with explicit solvent and counterions shows a comparable average minor groove width (McConnell and Beveridge, 2000), suggesting that the groove width may differ in solution compared to the crystalline state, possibly due to differing dielectric environments. In addition, both the DUPLEX-generated dodecamer structure and the dynamical simulations of McConnell and Beveridge (2000) reproduce the narrowness in the central AT-tracts of the DNA sequence, even though the average major groove width in the present study is smaller than the crystal average, and that of the molecular dynamics simulation is larger.

The recent NMR solution structure of the Dickerson-Drew dodecamer in an aqueous liquid crystalline phase (Tjandra et al., 2000) provides an opportunity for comparing the DUPLEX-generated dodecamer structure with a highly resolved DNA solution structure. As shown in Figs. 4-5, and A1 (Appendix), the computed structure is overall very close to the lowest energy and best representative member of the ensemble of five solution structures. In particular, the minor groove width is somewhat wider on average in the solution structure than in the crystals, and this is reproduced with the revised force field. Interestingly, the efficiently generated dodecamer produced by DUPLEX reproduces sequence-dependent conformational features in the NMR structure at least as closely as other force fields with state-of-the-art molecular dynamics simulations using explicit solvent and counterions (Figs. 4-6 and A1; Appendix) (Cheatham and Young, 2001).

	CONCLUSIONS

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION CONCLUSIONS REFERENCES

The improved force field in the present work satisfactorily reproduces the sequence-dependent structural features of the Dickerson-Drew dodecamer, namely the average geometry of complementary and neighboring base pairs, as well as the less clear sequence-dependent variation of torsion angles and pseudorotation parameters. The revised force field also reproduces the relative sequence-dependent trends of dimer tilt, roll, and twist angles found in the dynamical solution structure based on a 15-ns molecular dynamics study with AMBER (McConnell and Beveridge, 2000), although, as reported to date, the AMBER force field yields undertwisted DNA duplexes with neighboring base pairs unwound by ~3-4° compared to those in crystal structures (Cheatham and Young, 2001). We also note that Mazur (1998) has reported a successful reproduction of the original Dickerson-Drew dodecamer crystal structure (Wing et al., 1980; Drew et al., 1981) using a combination of partial explicit hydration in the grooves and implicit treatment of counterions according to Manning's theory (Manning, 1978), as in our approach, but with a linear distance-dependent dielectric function, (r) = r. The 5-ns trajectories based on this strategy converge to conformations close to that of the crystal structure in terms of root-mean-square atomic deviations, but the degree to which the sequence-dependent conformational parameters are reproduced remains to be reported.

Given the computational cost involved in molecular dynamics studies and the relatively robust performance of the revised DUPLEX force field, it is reassuring that the implicit treatment of solvent and counterions can still play an important role in molecular mechanics calculations. These computational shortcuts are particularly valuable for the extensive conformational searches with large numbers of energy minimization trials that must be carried out in studying carcinogen-DNA adducts; such detailed searches are beyond the present-day capabilities of molecular dynamics or Monte Carlo calculations with explicit solvent and counterions. In addition, NMR solution structures of carcinogen-damaged DNAs, which can be efficiently delineated in conjunction with molecular mechanics searches using distance restraints, benefit from a reliable force field with implicit treatment of solvent. Molecular dynamics simulations, using all-atom force fields such as AMBER (Case et al., 1999) with explicit solvent and counterions are, of course, essential for deciphering the structural codes that govern protein-DNA interactions, and for estimating conformational and binding free energies (Kollman et al., 2000). The force field treatments described here apply in the context of the DUPLEX molecular mechanics program, but similar approaches can be conveniently tested and modified, both in other force fields and with other benchmarks. We anticipate that these simple and efficient treatments will be useful in elucidating structure-function relationships for other applications, such as DNA-ligand interactions.

The Appendix contains five figures (Figs. A1-A5) presenting the variation of sequence-dependent base-pair parameters, groove widths, and depths in computed and experimental structures of the Dickerson-Drew dodecamer, along with one table (Table A1) presenting van der Waals' parameters, r_min and V_min, of nucleic acid atom pairs in different force fields.

View larger version (32K): [in this window] [in a new window]   FIGURE A1   Variation of additional sequence-dependent base-pair parameters of high resolution Dickerson-Drew dodecamer crystals (Table 2), the dodecamer structure in an aqueous dilute liquid crystalline phase (PDB_ID: 1DUF, model 1) (Tjandra et al., 2000), and the DUPLEX-generated B-DNA structures of the same sequence, using the original and improved force fields. Values calculated using RNA (Babcock et al., 1993). Color coding and symbols are the same as in Fig. 4.

View larger version (24K): [in this window] [in a new window]   FIGURE A2   Variation of selected sequence-dependent base-pair parameters of high-resolution Dickerson-Drew dodecamer crystals (Table 2), and the DUPLEX-generated B-DNA structures of the same sequence, using the original dielectric function (Hingerty et al., 1985) and the function of Young et al. (1998), both with the improved pseudorotation profile and relocated Lennard-Jones potential energy minimum. Values calculated using RNA (Babcock et al., 1993). Color coding: blue circles, crystals; gray triangles, original function; magenta squares, Young et al. (1998) function.

View larger version (21K): [in this window] [in a new window]   FIGURE A3   Variation of groove widths and depths of high-resolution Dickerson-Drew dodecamer crystals (Table 2), and the DUPLEX-generated B-DNA structures of the same sequence, using the original dielectric function (Hingerty et al., 1985) and the function of Young et al. (1998), both with the improved pseudorotation profile and relocated Lennard-Jones potential energy minimum. Values calculated using Curves 5.2 (Lavery and Sklenar, 1988, 1989). Color coding and symbols are the same as in Fig. A2.

View larger version (27K): [in this window] [in a new window]   FIGURE A4   Variation of selected sequence-dependent base-pair parameters of high-resolutoin Dickerson-Drew dodecamer crystals (Table 2), the dodecamer structure in an aqueous dilute liquid crystalline phase (PDB_ID; 1DUF, model 1) (Tjandra et al., 2000), and the DUPLEX-generated B-DNA structure of the same sequence, obtained with the improved force field using original and revised phosphorus-phosphorus (P ··· P) well depths. Values calculated uisng RNA (Babcock et al., 1993). Color coding: blue circles, crystals; green inverted triangles, 1DUF; gray triangles, original P ··· P interactions; red squares, revised P ··· P interactions.

View larger version (23K):