1108
Chem. Res. Toxicol. 2005, 18, 1108-1123
Free-Energy Perturbation Methods to Study Structure and Energetics of DNA Adducts: Results for the Major N2-dG Adduct of Benzo[a]pyrene in Two Conformations and Different Sequence Contexts Sushil Chandani, Chiu Hong Lee, and Edward L. Loechler* Biology Department, Boston University, Boston, Massachusetts 02215 Received December 20, 2004
The potent mutagen/carcinogen benzo[a]pyrene (B[a]P) is activated to (+)-anti-B[a]PDE, which induces a variety of mutations (e.g., G f T, G f A, etc.) via its major adduct [+ta]B[a]P-N2-dG. One hypothesis is that adducts (such as [+ta]-B[a]P-N2-dG) induce different mutations via different conformations, probably when replicated by different lesion-bypass DNA polymerases (DNAPs). We showed that Escherichia coli DNAP V was responsible for G f T mutations with [+ta]-B[a]P-N2-dG in a 5′-TGT sequence (Yin et al., (2004) DNA Repair 3, 323), so we wish to study conformations of this adduct/sequence context by molecular modeling. The development of a CHARMM-based molecular dynamics (MD) simulations protocol with free-energy calculations in the presence of solvent and counterions is described. A representative base-pairing and base-displaced conformation of [+ta]-B[a]P-N2-dG in the 5′TGT sequence are used: (1) BPmi5, which has the B[a]P moiety in the minor groove pointing toward the base on the 5′-side of the adduct, and (2) Gma5, which has the B[a]P moiety stacked with the surrounding base pairs and the dG moiety displaced into the major groove. The MD output structures are reasonable when compared to known NMR structures. Changes in DNA sequence context dramatically affect the biological consequences (e.g., mutagenesis) of [+ta]B[a]P-N2-dG. Consequently, we also developed a MD-based free-energy perturbation (FEP) protocol to study DNA sequence changes. FEP involves the gradual “fading-out” of atoms in a starting structure (A) and “fading-in” of atoms in a final structure (B), which allows a realistic assessment of the energetic and structural changes when two structures A and B are closely related. Two DNA sequence changes are described: (1) 5′-TGT f 5′-TGG, which involves two steps [T:A f T:C f G:C], and (2) 5′-TGT f 5′-TGC, which involves three steps [T:A f T:2AP f C:2AP f C:G], where 2AP (2-aminopurine) is included, because T:2AP and C:2AP retain more-or-less normal pairing orientations between complementary bases. FEP is also used to evaluate the impact that a 5′-TGT to 5′-UGT sequence change might have on mutagenesis with [+ta]-B[a]P-N2-dG. In summary, we developed (1) a CHARMM-based molecular dynamics (MD) simulations protocol with free-energy calculations in the presence of solvent and counterions to study B[a]P-N2-dG adducts in DNA duplexes, and (2) a MD-based free-energy perturbation (FEP) protocol to study DNA sequence context changes around B[a]P-N2-dG adducts.
Introduction Polycyclic aromatic hydrocarbons (PAHs) are a class of substances produced by incomplete combustion and are found ubiquitously in the environment, with the beststudied example being benzo[a]pyrene (B[a]P)1 (1-7). B[a]P is a potent mutagen and carcinogen that is metabolized to its carcinogenic form, (+)-anti-B[a]PDE, which reacts with duplex DNA, giving one major adduct, [+ta]-B[a]P-N2-dG) (Figure 1) (8, 9). (Other adducts can * Corresponding author: Biology Department, Boston University, Boston, MA 02215. Phone, 617-353-9259; fax, 617-353-6340; e-mail,
[email protected]. 1 Abbreviations: B[a]P, benzo[a]pyrene; (+)-anti-B[a]PDE, 7R,8Sdihydroxy-9S,10R-epoxy-7,8,9,10-tetrahydrobenzo[a]pyrene; [+ta]B[a]P-N2-dG, the major adduct of (+)-anti-B[a]PDE, formed by trans addition of N2-dG to (+)-anti-B[a]PDE; [-ta]-, [+ca]-, and [-ca]-B[a]PN2-dG, other stereoisomers of B[a]P-N2-dG; PAH, polycyclic aromatic hydrocarbon; NMR, nuclear magnetic resonance; BPmi5, a B[a]P-N2dG conformation with the B[a]P moiety in the minor groove pointing toward the base on the 5′-side of the adduct; Gma5, a B[a]P-N2-dG conformation with the B[a]P moiety stacked with the surrounding base pairs and with the dG moiety displaced into the major groove; MD, molecular dynamics.
Figure 1. Structure of [+ta]-B[a]P-N2-dG. Torsion angles are defined as follows: R′, N3-C2-N2-C10(BP); β′, C2-N2-C10(BP)-C9(BP); and χ, C2′-C1′-N9-C4.
form via other modes of metabolism (10-13).) Most carcinogens are active by causing mutations, and it is generally believed that mutations induced by PAHs in general, and B[a]P in particular, can cause mutations relevant to cancer causation (representative references
10.1021/tx049646l CCC: $30.25 © 2005 American Chemical Society Published on Web 06/23/2005
Free-Energy Perturbation of B[a]P adducts
include refs 14-20), and are important in human cancer (e.g., consider refs 21-24). Mutagenesis studies have been done with racemic (()anti-B[a]PDE (25-29), as well as with pure (+)-anti-B[a]PDE in Escherichia coli (30-32) and in mammalian Chinese hamster ovary (CHO) cells (33-36). In E. coli, base substitutions were most prevalent with (+)-anti-B[a]PDE (30-32) and principally occurred at G:C base pairs, where the qualitative and quantitative patterns of mutagenesis were influenced by DNA sequence context (see below). Furthermore, studies from our laboratory (3741) and other laboratories (42-46) have shown that DNA sequence changes can affect both the qualitative and quantitative patterns of mutagenesis by the major adduct [+ta]-B[a]P-N2-dG and related adducts. One hypothesis for this mutational complexity is that an adduct can adopt multiple conformations in DNA, where each conformation can cause a different kind of mutation, with adduct conformation being controlled by factors such as DNA sequence context (31, 47-49). Furthermore, it is clear that the DNA polymerase involved in adduct bypass affects the bypass outcome (50-63), including in the case of PAH adducts (64-77). A variety of physical studies, notably using NMR (see below) and fluorescence spectroscopy (78-86) have shown that [+ta]-B[a]P-N2-dG (and related adducts) can adopt multiple DNA conformations. Molecular modeling/computational chemistry also has proven to be an effective method to probe the conformational complexity of B[a]PN2-dG adducts (87-104; discussed at greater length in refs 100, 105). Broyde and colleagues have studied this most extensively (79, 80, 84-86, 100-111), notably in their pioneering studies of these adducts when modeled in solution (106, 107) and their most recent studies on adduct conformations in the active sites of DNA polymerases (108-111). On the basis of molecular modeling studies, we suggested that [+ta]-B[a]P-N2-dG could adopt at least 16 classes of conformations in ds-DNA (112-117). To study these conformations, we developed a four-step molecular modeling procedure (“24P”), and we showed it to be remarkably accurate in predicting the lowest energy conformations in five cases where the lowest energy conformation is known from a NMR study (116). However, 24P is an in vacuuo method with obvious limitations. Herein, we describe a CHARMM-based protocol to study [+ta]B[a]P-N2-dG in water using molecular dynamics (MD) with free-energy calculations, which follows closely the established methods of Broyde and colleagues (106-111). We chose a representative base-pairing (BPmi5) and base-displaced (Gma5) conformation for the study herein. “BPmi5” indicates that the B[a]P moiety of the adduct is in the minor groove and pointing toward the base on the 5′-side of the adduct, while the dG moiety forms a moreor-less normal G:C base pair. On the basis of NMR studies, [+ta]-B[a]P-N2-dG preferentially adopts the BPmi5 conformation in two DNA sequence contexts (78, 79). Gma5 is one of four possible “base-displaced” conformations, where “Gma5” indicates that the dG moiety is in the major groove and the “5” designates that the pyrene-face of the B[a]P adduct shown in Figure 1 is pointing toward the base on the 5′-side. Gma5 was preferentially adopted by [-ca]-B[a]P-N2-dG (86). In addition, when [+ta]-B[a]P-N2-dG was studied by NMR in a 5′-TG sequence, the dominant conformation was BPmi5 (78); however, a minor conformation was also
Chem. Res. Toxicol., Vol. 18, No. 7, 2005 1109
present, and it was noted probably to be base displaced, although the structure could not be determined (78). Gma5 is a reasonable candidate for this minor conformation, since it is adopted by [-ca]-B[a]P-N2-dG (86), which has the identical adduct bond stereochemistry as [+ta]B[a]P-N2-dG, and since Gma5 was computed to be the lowest energy, base-displaced conformation in our previous modeling studies of [+ta]-B[a]P-N2-dG in several 5′TG sequences (113, 116, 117), although none of these arguments are definitive. A 5′-TGT sequence was used herein, because we recently showed that [+ta]-B[a]PN2-dG induces an even mixture of G f T and G f A mutations, where the former but not the latter depend on E. coli DNAP V (118), and we wished the work herein to be relevant to our mutational studies. DNA sequence context can affect the biological and conformational properties of B[a]P-N2-dG adducts. The mutagenic pattern of [+ta]-B[a]P-N2-dG can change dramatically (37-46); for example, in different sequence contexts, its mutations can be >95% G f T (37), ∼95% G f A (40), or ∼80% -1 frameshifts (119-121). Subtle changes in sequence context affect conformation; for example, B[a]P-N2-dG adducts changed conformation when a 5′-CG sequence was changed to 5′-meCG sequence (122, 123). Sequence context can also affect subtler parameters, such as the roll angle in the vicinity of the adduct (124). It is not clear whether sequence context affects adduct properties directly or indirectly; for example, DNA sequence could affect water structure, which in turn could affect adduct conformation. The structure and energetics of adduct conformations in different DNA sequence contexts can be studied by doing independent MD simulations in each sequence context. However, two MD trajectories may evolve differently, such that differences in structures (and energies) may reflect the individualities of the trajectories more than inherent differences. When two structures differ only slightly (e.g., with related sequence contexts) free-energy perturbation (FEP) has often proven valuable and offers an alternative approach to circumvent the vagaries of independent MD trajectories. FEP (also called “thermodynamic simulation” or “thermodynamic integration”) can be thought of as the gradual fading-out of an ensemble of atoms in an initial state (A) concomitant with the fading-in of replacement atoms in the final state (B), using a series of closely spaced, hybrid intermediates (125, Materials and Methods). The total free-energy change in going from A to B is the sum of the free-energy changes between the intermediates, since free energy is a state function. FEP methodology has been used to investigate the properties of proteins and protein folding, as well as protein affinity for small molecules, ions, hormones, and peptide inhibitors (125-137). FEP has also been used to study DNA mismatches in duplex DNA (138) and the role of sequence context on the B- to Z-DNA transition (139, 140). FEP has helped rationalize the affinity of transcription factors for different DNA sequences (141, 142), as well as contributed to our understanding of T7 DNA polymerase mechanism (143, 144), and the role of sequence context on RNA hairpin stability has also been studied with FEP methods (145, 146). In summary, the goals of the study herein were twofold. First, we developed a CHARMM-based protocol to study [+ta]-B[a]P-N2-dG in a base-paired conformation (BPmi5) and a base-displaced conformation (Gma5) using MD in water with free-energy calculations. Second,
1110
Chem. Res. Toxicol., Vol. 18, No. 7, 2005
we developed FEP methods to study the relative free energies of B[a]P-N2-dG adducts in closely related DNA sequence contexts, in particular to study [5′-TGT f 5′TGC] and [5′-TGT f 5′-TGG] for [+ta]-B[a]P-N2-dG in the BPmi5 and Gma5 conformations. We also used FEP to estimate that [+ta]-B[a]P-N2-dG in the Gma5 conformation is destabilized slightly less in going from a 5′TG to a 5′-UG sequence compared to the BPmi5 conformation, a finding that is consistent with the observation that the dG moiety of Gma5 stacks favorably on the methyl group of the 5′-T in the 5′-TG sequence. This effect is modest, however, and the finding is considered in light of our recent mutagenesis findings that G f T mutations declined ∼2-fold in going from 5′-TGT to 5′-UGT.
Materials and Methods Viewing and manipulations of structures were done with Insight II 97.0 from Accelrys, Inc., a subsidiary of Pharmacopeia, Inc., as we have done in the past (147). Structures were viewed on O2 workstations (Silicon Graphic, Inc.), while computations were conducted on either an Origin 200 (Silicon Graphic, Inc.) or one of Boston University’s IBM Power4 systems (p690 or p655). CHARMM was obtained from Dr. Martin Karplus (Harvard University), and the CHARMM 27 force field was used (148-150). We followed the method of Yan et al. (106, 107) in our MD simulations with free-energy calculations in the presence of solvent and counterions, except where noted. A nine-base pair duplex (5′-CGCTGTCAT/5′-ATGACAGCG without 5′-terminal phosphates) was used, as preliminary work showed that 9-mers gave DNA that appeared reasonable in water, and 9-mers are a sensible compromise between simulating long stretches of ds-DNA (minimizing end effects) versus permitting calculations to be completed in reasonable time frames. Following our previously established protocol (112-117), a starting structure for [+ta]-B[a]P-N2-dG in either a base-paired conformation (BPmi5) or a base-displaced conformation (Gma5) was incorporated into a DNA molecule derived from canonical fiber diffraction B-DNA coordinates (151). In all cases, starting structures were constructed to be consistent with NMR structures. For example, the hydrogen on the hydroxyl groups in the B[a]P moiety was placed in the same orientation as observed by NMR for BPmi5 (79), and the orientation of the dC moiety opposite [-ca]-B[a]PN2-dG was positioned in the same orientation as observed by NMR for Gma5 (86). Parameters for [+ta]-B[a]P-N2-dG were based on our previous work (112-117) and were consistent with the CHARMM 27 force field. Identical parameters were used for both the BPmi5 and Gma5 conformations. Point charges were at the center of each atom, and each nucleotide had a total charge of -1, as prescribed by CHARMM. A sodium ion was placed 4 Å from each phosphorus (16 total) in the oxygenphosphorus-oxygen plane. To optimize the position of the sodium ion, conjugate gradient minimization was performed first with (50 steps) and then without (200 steps) fixing the positions of the atoms in the DNA, where the latter also resolves, for example, van der Waals contacts without moving the structure significantly toward an in vacuuo structure. TIP3P water (152) was then added, and subsequent calculations were performed with periodic boundary conditions at atmospheric pressure with a 0.2 ps coupling parameter at constant temperature (cpt). Each 9-mer DNA duplex was centered in a water box, whose dimensions were 48 Å × 40 Å × 40 Å, and which extended at least 8 Å from the DNA in each dimension. The structures for BPmi5/TGT, Gma5/TGT, BPmi5/UGT, and Gma5/ UGT contained 2279, 2272, 2277, and 2266 water molecules, respectively. The water was allowed to relax into position by fixing the coordinates of the DNA and the sodium ions and by performing conjugate gradient minimization (tolerance gradient
