Article pubs.acs.org/JCTC
Dominant Pathways of Adenosyl Radical-Induced DNA Damage Revealed by QM/MM Metadynamics Paweł Wityk,†,§ Miłosz Wieczór,‡,§ Samanta Makurat,† Lidia Chomicz-Mańka,† Jacek Czub,*,‡ and Janusz Rak*,† †
Faculty of Chemistry, University of Gdańsk, Wita Stwosza 63, 80-952 Gdańsk, Poland Department of Physical Chemistry, Faculty of Chemistry, Gdańsk University of Technology, Narutowicza 11/12, 80-233 Gdańsk, Poland
‡
S Supporting Information *
ABSTRACT: Brominated nucleobases sensitize double stranded DNA to hydrated electrons, one of the dominant genotoxic species produced in hypoxic cancer cells during radiotherapy. Such radiosensitizers can therefore be administered locally to enhance treatment efficiency within the solid tumor while protecting the neighboring tissue. When a solvated electron attaches to 8-bromoadenosine, a potential sensitizer, the dissociation of bromide leads to a reactive C8 adenosyl radical known to generate a range of DNA lesions. In the current work, we propose a multiscale computational approach to elucidate the mechanism by which this unstable radical causes further damage in genomic DNA. We employed a combination of classical molecular dynamics conformational sampling and QM/MM metadynamics to study the thermodynamics and kinetics of plausible reaction pathways in a realistic model, bridging between different time scales of the key processes and accounting for the spatial constraints in DNA. The obtained data allowed us to build a kinetic model that correctly predicts the products predominantly observed in experimental settingscyclopurine and β-elimination (single strand break) lesionswith their ratio and yield dependent on the effective lifetime of the radical species. To date, our study provides the most complete description of purine radical reactivity in double stranded DNA, explaining the radiosensitizing action of electrophilic purines in molecular detail as well as providing a conceptual framework for the computational modeling of competing reaction pathways in biomolecules.
1. INTRODUCTION Thanks to its efficacy, radiotherapy is one of the most common modalities used in cancer treatment.1 Unlike UV radiation, ionizing radiation can penetrate inside the patient’s body, killing tumor cells located in deeper tissues. The hardness of ionizing radiation leads inevitably to increased mutation rate and even secondary cancers in the surrounding healthy cells that absorb part of the therapeutic dose.2,3 Moreover, since the cells adjacent to the tumor have a normal oxygen level, they are about 3 times more sensitive to the effects of irradiation than the hypoxic solid cancer cells.4 This apparent stalemate can be resolved by the directed application of radiosensitizers compounds that locally increase the damaging effects of ionizing radiation, allowing for lower therapeutic doses and hence for limiting the devastating side effects of radiotherapy.5 One important class of such radiosensitizing agents that can be incorporated into DNA in vivo are the widely known halogen derivatives of pyrimidine nucleosides, 5-bromo(BrdU) and 5-iodouridine (IdU).6−12 However, not only halouridines should be able to sensitize DNA to radiationinduced hydrated electrons (ehyd−) in vivo. In analogy to BrdU and IdU, other nucleotide haloderivativessuch as 8brominated purines (BrdA, BrdG) or 5-brominated cytosine © 2017 American Chemical Society
(BrdC)have been shown to undergo dissociative electron attachment (DEA) that results in the formation of the bromide anion (Br−) and a reactive radical located at the C5 position of pyrimidines or C8 of purines.13 One trait common to these potential sensitizers is that this reactive moiety subsequently initiates a cascade of reactions leading to the formation of DNA lesionschemical modifications that compromise the stability and readability of genetic information stored in DNA.14,15 In particular, primary reactions of the radical involve the adjacent hydrogen atoms of the deoxyribose moiety, creating many possible stabilization pathways that can compete with one another. These individual primary reactions then lead to the formation of different products, preselecting for lesions that might elicit different cellular responses. To date, several studies attempted to compare the thermodynamic and kinetic feasibility of individual pathways by employing minimal QM model systems either with implicit solvent or in vacuo.16,17 A growing body of research, however, highlights the crucial influence of proper modeling of reaction environment on the reliability of the calculated kinetic and thermodynamic Received: September 19, 2017 Published: November 9, 2017 6415
DOI: 10.1021/acs.jctc.7b00978 J. Chem. Theory Comput. 2017, 13, 6415−6423
Article
Journal of Chemical Theory and Computation
Figure 1. (A) View of the full simulation system consisting of a 21-base-pair dsDNA helix solvated with explicit water and Na+/Cl− ions. (B) The QM subsystem comprising the DNA trinucleotide. Hydrogen atoms that undergo abstraction in primary reactions are marked in red, orange, blue, and green, consistently with color-coding throughout the manuscript. An isosurface of positive spin density localized on the C8 carbon is shown in purple, and sodium ions are shown as pink spheres. (C) Dissociative electron attachment to 8-bromoadenine that leads to the formation of a reactive radical in sensitized DNA. (D) Schematic representation of all primary (left) and secondary (right) reactions explored in this work. For detailed representation of the reactions, see Figures S1−S7 in the SI.
properties.18−20 A unified approach, incorporating both spatial restraints imposed by the geometry of DNA and the electrostatic effects due to the inhomogeneous solvent distribution into a single multipathway model, would therefore be needed to properly predict the fate of the reactive species and thus reveal the molecular mechanism as well as the resulting products of DNA sensitization by halogenated nucleotides. From the sensitization point of view, the most relevant products of DEA and subsequent reactions are the DNA lesions that induce a strong cellular response. Chief among them are DNA single strand breaks (SSBs) since they may be readily converted into potentially lethal double strand breaks (DSBs) either (i) during repair of base damage by the base excision repair (BER) mechanism, when a transient single-strand nick is combined with a pre-existing SSB on the opposite strand,21 or (ii) when a replication fork collapses upon encountering an SSB in the S phase.21 Other types of damage initiated by nucleobase radicals are inter-22,23 and intrastrand24−29 cross-links that impair replication by stalling the replication fork or introducing
mutations, as well as destabilize the DNA double helix.30−35 Finally, cyclopurine lesions, unique to purine nucleobases, are formed in the stabilization pathway of the 8-yl radical of adenosine or guanosine.36−38 If not repaired, these products can impair gene transcription and protein−DNA binding.39 Here, we describe the outcome of our QM/MM metadynamics study on a double stranded DNA fragment comprising the reactive C8 adenosyl radical. Using combined multiscale state-of-the-art molecular dynamics (MD) simulations, we studied possible hydrogen atom transfer reactions between the radical and the adjacent deoxyribose moieties, as well as the following secondary reactions that lead to DNA lesions. The calculated free energy profiles, along with kinetic data derived from classical MD, allowed us to extract relative rate constants and construct a system of coupled kinetic equations that predicts the yields of individual products. Consistent with experimental data,36,7 our results suggest that the cyclization of purine to 5′,8-cyclo-2′-deoxyadenosine and the formation of single strand breaks in the β-elimination pathway are the major stabilization paths of the adenosine-8-yl 6416
DOI: 10.1021/acs.jctc.7b00978 J. Chem. Theory Comput. 2017, 13, 6415−6423
Article
Journal of Chemical Theory and Computation radical. Our model also shows how the yield and ratio of products might depend on the (environment-dependent) lifetime of the radical species as well as aligns with the growing body of research40,41 highlighing the relevance of extensive conformational sampling in comprehensive modeling of chemical reactions within biomacromolecules.
2. METHODS 2.1. Simulation System. The DNA 21-mer was built in an ideal B-DNA conformation using the X3DNA utility f iber.42 Its sequence (5′-CGGCGTGAGCACGCTGAGACT-3′) corresponds to a fragment of the BRCA1 oncogene coding region. The hydrogen atom at the C8 position of the central adenine (bolded in the sequence) was removed to yield the adenosyl radical produced, e.g., during debromination of 8-BrdA.16 For the purpose of classical MD simulations, the adenosyl radical moiety was reparametrized to account for changes in charge distribution by calculating the ESP fitted charges at the HF/631G(d,p) level of theory (i.e., consistent with the Amber force field) using Gaussian g09-D. Amber99sb force field parameters43 with the parmbsc0 correction for nucleic acids44 implemented in the Gromacs suite45 were used for all classical MD simulations. The DNA molecule was embedded in a 60 × 60 × 110 Å perpendicular box and solvated with 8118 TIP3P water molecules. A total of 64 sodium and 24 chloride ions were added to ensure charge neutrality at the physiological concentration of 0.154 M (for the full simulation system, see Figure 1A). The Plumed plugin46 was used to apply a harmonic potential restraining the XY distance between the top and bottom base pairs so that the DNA molecule remained in an upright position (aligned along the Z-axis). After initial energy minimization, a 1-μs equilibrium simulation was run to obtain a proper ensemble of DNA conformational states. 2.2. QM/MM Starting Structures. The long classical MD simulation served to generate initial frames for the subsequent QM/MM metadynamics simulations aimed at examining the reactions initiated by the adenosyl radical in DNA. To prevent metadynamics from introducing large distortions to the DNA structure due to the deposited bias, we selected the starting geometries that already favored the considered reactions as indicated by low C8−Hx (x = 2′, 3′, 5′, 2′(n)) distances. In practice, the initial frames were selected based on histograms of C8−Hx distances so as to correspond to the first percentile of the distance distribution; i.e., geometries were chosen in which the C8−Hx distance was smaller than in 99% of trajectory frames (see Figure 2). 2.3. QM/MM Procedure. All QM/MM simulations were performed using the cp2k suite.47 Notably, cp2k employs an auxiliary plane-wave basis set to achieve a significant speedup of the calculations, which in practice precludes the use of modern hybrid DFT functionals for larger systems due to the slow calculation of HF exchange. Since the PBE functional typically employed in cp2k ab initio MD simulations is known to perform poorly in calculations of hydrogen transfer barrier heights,48 we used the more advanced Minnesota local functional M06-L, shown to yield more than 2-fold smaller mean unsigned errors on the standard test set (4.15 vs 9.31 kcal/mol). To this end, cp2k was compiled with the LibXC library providing access to an extended set of functionals. To obtain a correct representation of the core density, the Goedecker−Teter−Hutter (GTH) pseudopotentials for each atomic species were reparametrized for M06-L using the ATOM utility provided with cp2k. GTH pseudopotentials were
Figure 2. (A) Distributions of distances between the C8 radical and individual hydrogen atoms considered in the primary reactions. Inset: zoom on the peak in the distribution corresponding to the reactive (H5′R) state. (B) Positions of H5′ atoms with respect to the C8 radical, observed in the equilibrium simulation, with red spheres denoting the reactive (H5′R) and blue spheres, the nonreactive state.
used in conjunction with the DZVP atomic basis set, and plane waves were cut off at 300 Ry. The D3 correction by Grimme49 was employed to ensure proper description of dispersion forces, in particular those related to π−π stacking of DNA bases. 2.4. Free Energy Calculations. In QM/MM metadynamics simulations, the QM subsystem contained three DNA bases (the central adenine/adenosyl radical and two flanking cytosine residues) connected by three phosphodiester bonds (see Figure 1B), and three backbone-bound sodium ions were included to neutralize the net charge of the QM box, yielding a total of 96 atoms. Link atoms were introduced to cap dangling C4′−C5′ bonds at the QM/MM boundary. The QM system was embedded in a 20 × 20 × 20 Å cubic box, with reflective boundaries located at 1.5 Å from the box edges. The Quickstep module of cp2k was used to calculate energies and forces in the QM subsystem, and interactions in the MM system were calculated using the FIST module with the Amber99SB/ parambsc0 force field parameters. Simulations were performed in the NVT ensemble with a time step of 0.5 fs. The CSVR thermostat with a time constant of 5 fs was used to maintain the temperature at 300 K, with separate treatment of the QM and MM subsystems. The GEEP module in cp2k was used to generate the Gaussian expansion of the MM potential, with 12 Gaussians used in the expansion. All primary (hydrogen transfer) metadynamics simulations were preceded by geometry optimization at the QM/MM level and 7.5 ps equilibration in cp2k with no restraints. The actual well-tempered metadynamics runs were carried out with a bias factor of 35, with individual simulation lengths of 150 ps corresponding to 300 000 simulation steps. The historydependent bias was deposited every 25 steps as Gaussians with a height of 0.5 kcal/mol and scaling factor of 0.5. Quartic walls were introduced to prevent systems from sampling beyond the chemically relevant region of the configurational space. Initial frames for the simulations of secondary reactions (DNA lesions) were selected after the first hydrogen transfer was observed. Since even well-tempered metadynamics tends to produce free energy profiles that fluctuate around the actual profile in finite-length runs, the obtained free energies were probability density-averaged over the last 50 ps to produce converged profiles using a modification of the approach proposed by Laio et al.: τ1 1 G (ξ ) = log exp[−βG(ξ ; τ )] dτ β(τ1 − τ0) τ0
∫
6417
DOI: 10.1021/acs.jctc.7b00978 J. Chem. Theory Comput. 2017, 13, 6415−6423
Article
Journal of Chemical Theory and Computation where G(ξ;τ) is the free energy profile estimated from metadynamics at time τ, β is equal to (kBT)−1, and τ0 and τ1 denote the interval over which averaging was performed.50 In primary reactions, the reaction coordinate was defined as the difference between the Cx−Hx and C8−Hx distances (x = 2′, 3′, 5′, 2′(n)), with Cx being the deoxyribose carbon atom from which hydrogen abstraction took place and Hx, the hydrogen atom being transferred. In secondary bond-breaking or bond-forming reactions, interatomic distances between the atoms involved in the bond of interest were used as reaction coordinates (see Figure 3 for a schematic description). The free
distribution of profiles evaluated at the free energy barrier and in the product state. 2.5. Kinetic Modeling. For all chemical reactions, preexponential factors in the transition state theory equation kon = ν exp[(−ΔG⧧)/(kBT)], which can be thought of as vibrational frequencies in the reactant state, were calculated from the MD trajectories as an inverse of the autocorrelation times calculated for the reacting partners’ distance.51,52 From the obtained free energy profiles and autocorrelation MD data, rate constants were calculated and used to solve a set of coupled ordinary differential equations describing the kinetics of the studied reactions53 (see Figure S9 and eq S1). Rate constants for the reverse reactions were calculated on the basis of equilibrium constants, K = exp[(−ΔG°)/(kBT)], as koff = kon/K. To account for the eventual termination of the radical reactions, we introduced an additional set of closed-shell (stable) species that can be formed from any radical (unstable) species in a termination reaction, with a rate constant kter of 1, 103, or 106 s−1 chosen based on the work of Chatgilialoglu et al.36 Populations of all thus defined species (see Table S1 for an exhaustive list of chemical moieties and parameters used in the modeling) were propagated numerically starting from the C8 radical as the initial species (100% at t = 0) for at least 10/kter (10, 10−2, or 10−5 s, in respective order). Among the studied hydrogen transfer reactions, only the C8−H5′ distance histogram indicated the presence of two distinct populations, the reactive H5′R (C−H distance centered about 3 Å) and nonreactive H5′NR (C−H distance centered above 5 Å). The hydrogen abstraction from the C5′ position was therefore modeled as a two-step process, involving the initial transition from the nonreactive to the reactive conformation (with rate constants kn→r for the forward and kr→n for the backward transition) and the subsequent chemical step. The rate constants kn→r and kr→n were calculated from mean transition times, with states (reactive vs nonreactive) inferred from the interatomic distances as hidden variables using a two-state hidden Markov model (see Figure S8). In practice, though, these rate constants were typically 2−3 orders of magnitude higher than for other reactions, so the hydrogen abstraction from C5′ was modeled using an effective preexponential factor νeff = (kn→r/kr→n)ν, i.e., by assuming instantaneous equilibration of the two states. For error analysis, the population modeling was repeated 100 times with ΔG‡ values sampled from the normal distribution according to the standard deviation taken from the bootstrap analysis.
Figure 3. Free energy profiles for the primary (hydrogen abstraction) and secondary (DNA lesions) reactions. Solid lines correspond to the free energies calculated using QM/MM metadynamics and the dotted lines, to Boltzmann-inverted histograms from classical MD simulations. “S” marks the substrate state in which the radical is localized at C8. The asterisk (*) marks the intermediate state in which the radical is localized on the sugar moiety. “P” marks one of the considered product states (lesions): abasic site (ABS), β-elimination product (BEL), ketone (KET), aldehyde (ALD), or cyclopurine (CYC).
energy profiles for possible products of 2′/2′(n) hydrogen abstractionβ-elimination and abasic site formationwere assumed to be independent of nucleotide sequence and hence were calculated only once, for the case of H2′ abstraction within the central nucleotide of the QM subsystem. In effect, the corresponding segments of the profiles in Figure 3 are identical. To extend the free energy profiles beyond the (reactive) interval sampled by metadynamics, all strand break products were parametrized at the MM level consistently with the initial simulation of the adenosyl radical, and additional 200 ns equilibrium simulations were performed for the DNA 21mer bearing each lesion. The extended free energy profiles for both primary and secondary reactions (see Figure 3) were obtained by Boltzmann inversion from histograms of interatomic distances in the classical simulations and served to provide a better estimate of the reaction free energy, in particular to account for the entropic contribution associated with the conformational freedom of the newly formed strand ends. To obtain the standard errors for the calculated free energies, a bootstrap analysis was performed by drawing with replacement a set of (total time/free energy autocorrelation time) free energy profiles from the metadynamics run and averaging the profiles aligned in the substrate state. The procedure was repeated 100 times, and the reported standard errors were calculated as the standard deviation of the obtained
3. RESULTS AND DISCUSSION Strand breaks and cyclopurine lesions are the products most frequently reported in experimental studies on DNA radiation damage14,21,54 Hence, much effort has been devoted to the elucidation of mechanisms leading to formation of these lesions.6 It is argued that this demanding and costly research could be greatly facilitated by a robust computational technique that would allow for testing of the underlying assumptions in a biologically relevant environment. Here, we show that by performing a comprehensive set of QM/MM and classical MD simulations, it is possible to reliably test the hypothesized mechanisms of DNA damage. Upon treatment with ionizing radiation in a hypoxic aqueous solution, 8-bromo-2′-deoxyadenosine (BrdA) incorporated into DNA undergoes efficient dissociative electron attachment (DEA) leading to bromide elimination and the formation of 6418
DOI: 10.1021/acs.jctc.7b00978 J. Chem. Theory Comput. 2017, 13, 6415−6423
Article
Journal of Chemical Theory and Computation Table 1. Computed Thermodynamic and Kinetic Properties of the Considered Reactionsa ΔG° [kcal/mol] ΔG‡ [kcal/mol] ν [109 × 1/s] k [1/s]
3′
2′(n)
2′
5′
ALD
BEL
−10 ± 0.9
−13 ± 0.7
−12 ± 0.5
−17 ± 0.7
−8 ± 1.6
cyclo
5 ± 1.1
ABS
−1 ± 1.3
KET
7 ± 0.8
1 ± 0.7
24 ± 1.1
15 ± 0.7
22 ± 0.5
8 ± 0.7
14 ± 1.1
23 ± 0.7
15 ± 0.8
17 ± 0.5
13 ± 0.3
4.18 1.11 × 10−8
12.31 1.27 × 10−1
18.35 1.42 × 10−6
0.02 2.77 × 101
4.08 5.58 × 10−2
104 1.43 × 10−8
104 1.03 × 10−2
104 3.55 × 10−4
104 3.01 × 10−1
The table reports free energy difference (ΔG°) and kinetic barrier heights (ΔG‡) taken from the free energy profiles, as well as the pre-exponential factors (ν) inferred from the corresponding autocorrelation times and rate constants (k) calculated from the transition state theory equation. ΔG’s are reported ± standard error as calculated by bootstrap analysis of metadynamics free energy profiles estimated at different simulation times. a
states, as shown in Figure 2A,B and Figure S8. The state found predominantly in the simulation is characterized by a large mean C−H distance of ∼5 Å and was hence dubbed “nonreactive” (H5′NR). Meanwhile, the other stateeven though significantly less populatedcorresponds to a mean C−H distance of only ca. 3 Å and should favor the reaction due to spatial proximity of the atoms involved (the “reactive” state, H5′R). We note that the observed time scale of the H5′NR → H5′R transition, on the order of 80 ns, highlights the importance of initial conformational sampling both for the choice of initial frames for QM/MM MD simulations and for the subsequent modeling of reaction kinetics. To investigate the thermodynamics and kinetics of the actual chemical reactions following the formation of the C8 adenosyl radical, we employed well-tempered metadynamics within a QM/MM MD framework. In case of primary reactions, the difference between the Cx−Hx and C8−Hx distances was chosen as the reaction coordinate, so that negative values of the coordinate correspond to the substrate and positive, to the product. After the first observed hydrogen atom abstraction by the C8 radical, the structure of the product (with the radical now on the sugar moiety) was used as a starting frame for the simulations of secondary reactions in separate metadynamics runs, according to the scheme presented in Figure 1D. Here, as secondary reactions involve the formation of DNA lesions, the lengths of the breaking/forming bonds were chosen as the reaction coordinate. In addition, equilibrium classical simulations of the final products were performed to probe the distribution of interatomic distances in broken DNA strands. The free energy profiles obtained in this way, shown in Figure 3, combine all considered reactions leading to specific DNA lesions, with primary and secondary reactions merged together after alignment of the respective free energy minima. As a result, the calculated thermodynamic driving forces and kinetic barriers allow direct comparison of the feasibility of individual pathways. As seen in Figure 3, while kinetic barriers (ΔG‡) for the primary reactions are different in all cases, ranging from 8 to 24 kcal/mol, the hydrogen abstraction itself is thermodynamically favorable with the free energy change (ΔG°) below −10 kcal/ mol (see also Table 1). This strongly suggests that under biologically relevant conditions, i.e., when individual pathways compete with one another at physiological temperature, the proportion of the formed DNA lesions is governed by relative rates of competing processes (kinetic control). Indeed, the high free energy barriers for the reverse reactions (>25 kcal/mol) would result in rate constants lower than 10−9/s, rendering the backward processes too slow to be observed. The profiles clearly indicate that both the most kinetically accessible and the most thermodynamically favored primary reaction is the H5′ hydrogen abstraction, characterized by a free energy barrier of 8
the C8 adenosyl radical (Figure 1C). The fate of this unstable aromatic radical is assumed to initially involve the abstraction of a hydrogen atom from the adjacent sugar moiety in what is hereafter referred to as primary reactions. This determines further reactivity of the system by activating multiple secondary reaction pathways that ultimately lead to DNA breaks and irreversible damage (see Figures S1−S7). As suggested in prior experimental and theoretical works,6 the following secondary reactions can follow the initial hydrogen abstraction to lead to DNA damage: (i) purine cyclization (CYC, see Figure S6), (ii) β-elimination (BEL see Figures S4 and S5), (iii) abasic site formation (ABS see Figures S1 and S2), and (iv/v) aldehyde (ALD see Figure S3) or ketone (KET see Figure S7) formation (see Figure 1D for schematic overview). It is worth noting that these pathways are markedly different from the ones typically observed in nonsensitized DNA, where the initial radical is most likely to reside on the solvent-exposed C4′ atom55 that is conformationally inaccessible to nucleotide-based sensitizers. Notably, the C8 adenosyl radical formed as the initial product of DEA is remarkably unstable, despite its location in a planar conjugated ring system. This can be explained by drawing a parallel between the reactivity of adenosyl radicals and the stability of the respective radicals in the imidazole system. Indeed, in the dehydroimidazole system, the Ncentered radical was shown to be the most stable form due to its π character, contrary to the σ-type C-centered radical,56 putting a rough estimate of 20 kcal/mol on the energy difference between the unstable C8 and the most stable N7 adenosyl radical. Such reactivity renders the C8 radical genotoxic, a desired trait in the design of radiosensitizers. To investigate how the adenosyl radical can initiate the damage of the DNA strand, we first performed a long, 1 μs unbiased classical MD simulation to search the conformational space for states that favor the primary reactions and would later serve as starting geometries for QM/MM calculations. In particular, we focused on the distribution of distances between the C8 adenosyl radical and the four deoxyribose hydrogen atoms, chosen based on previous studies13,17 (see Figure 2), as an important parameter in determining the feasibility of primary reactions. We found that the atoms with the lowest mean C8− Hx distance are the 2′ hydrogen atoms located on either the same nucleotide as the C8 radical orsurprisinglythe neighboring one in the 5′ direction (further referred to as H2′(n)). Due to an orders-of-magnitude difference between time scales accessible to electronic structure calculations and time scales involved in structural transitions, we selected initial frames based on histograms of C8−Hx distances so that they corresponded to the top first percentile of the respective distance distribution. Interestingly, the distance distribution of C8−H5′ reveals the existence of two distinct conformational 6419
DOI: 10.1021/acs.jctc.7b00978 J. Chem. Theory Comput. 2017, 13, 6415−6423
Article
Journal of Chemical Theory and Computation kcal/mol and a reaction free energy of −17 kcal/mol. The subsequent formation of a C5′ radical can result in formation of (i) cyclopurine or (ii) strand break (P−O bond rupture) with a terminal aldehyde (see Figure 3 and Figures S1−S7). Notably, the former type of damage also seems to be the most frequently observed type of damage induced by purine sensitizers,36,37 in agreement with our predictions. We speculate that the lowest free energy barrier as well as the most negative ΔG° found for H5′ abstraction might be attributed to the greater mobility of the C5′ atom, not constricted by the five-membered deoxyribose ring as in the case of C2′ and C3′. This not only allows H5′ to approach the C8 radical with more flexibility but also enables the newly formed alkyl radical to assume its preferred geometry and hence stabilize this intermediate state. Though the abstraction of H2′ from the adenosine nucleotide bearing the C8 radical is prevented by a high free energy barrier of 22 kcal/mol (Figure 3), the virtually identical reaction involving the neighboring sugar in the 5′ direction is much more kinetically favorable, with a free energy barrier of 15 kcal/mol. This primary H2′(n) atom abstraction might therefore compete with H5′ abstraction to produce a C2′ radical that can engage in two different secondary reactions: (i) C−N glycosidic bond rupture that leads to the formation of an abasic site or (ii) β-elimination that produces a strand break with a 5′ phosphate group and a double C2′−C3′ bond within the sugar moiety, as illustrated in Figures 3 and S5. It is worth noting that the almost identical ΔG° obtained for H2′(n) and H2′ can serve as a control: obtaining almost identical values (−12 vs −13 kcal/mol) for two independent simulations with virtually chemically identical products proves that the profiles are well converged after the 150 ps runs. On the other hand, the corresponding difference in barrier heights (15 vs 22 kcal/ mol) are large, and to explain them we resorted to equilibrium classical MD simulations. As revealed by a principal component analysis (PCA) performed on a subset of the entire system (adenosyl radical nucleotide and the neighboring sugar moiety), the dominant conformational motion of the DNA strandrevealed by the first principal componentinvolves close approach of the C8 radical to the H2′(n) atom, as shown by the orange line in Figure 4. Meanwhile, even though the H2′ atom often remains in closer proximity to C8, in the dominant conformational mode the two atoms do not move toward one another but rather pass by each other (blue line in Figure 4). This is also reflected in the previously discussed distance distributions (Figure 2), where the distribution in the case of H2′ is significantly narrower than for H2′(n), again implying a huge difference in the mobility along the respective reaction paths. This leads to an interesting example of how kinetics in large biological systems can be affected by the local intrinsic deformability of the biomolecule. Finally, the highest kinetic barrier (24 kcal/mol) corresponds to the abstraction of the H3′ atom. Here, the barrier height correlates directly with the equilibrium distribution of C8−H3′ distances, shifted significantly toward higher values (ca. 4.8 Å). Hence, in order for the C8−H3′ distance to be sufficiently small for the abstraction reaction to occur, the local conformation has to be significantly distorted, which comes at an energetic cost and contributes to the barrier height. Due to the kinetic limitation of this primary reaction, the subsequent thermodynamically favorable rupture of the O−P bond and release of the 3′ phosphate group or the formation of the ketone lesion (see Figure 3 and S7) are also extremely unlikely
Figure 4. PCA analysis performed on the minimal subset of the DNA strand (top) from the long equilibrium classical MD simulation. The interatomic distances (bottom) are calculated from an artificial trajectory showing progression along the dominant principal component. Note that in the H5′R reactive state, the dominant conformational motion of the DNA would also facilitate abstraction of the H5′ atom.
to occur since the competing primary reactions are markedly faster. As we found the hydrogen abstraction reactions to be controlled kinetically, the difference in barrier heights between the two most feasible reaction pathwaysH5′ and H2′(n) should be crucial for predicting the amount of products observed in experimental settings. Assuming identical preexponential factors for all reactions, the calculated barrier height difference of ca. 7 kcal/mol (see Table 1) would translate to a rate constant difference of 4−5 orders of magnitude, virtually ruling out the possibility of H2′(n) abstraction. However, to obtain a more reliable estimate, these rate constants should be corrected by including the long-time conformational dynamics of DNA extracted from equilibrium MD simulations. The difference in pre-exponential factors can be inferred from autocorrelation times to account for different frequencies with which the adenosyl radical approaches the hydrogen atoms. In addition, rate constants will be affected by the transition rate between the reactive and nonreactive states for H5′. Together, these corrections change this predicted difference to ca. 2 orders of magnitude, leaving the products of both pathways amenable to experimental detection. This conclusion agrees with the dominant view that the C5′ radical is one of the most abundant radicals present in the irradiated BrdA solutions (as opposed to the C4′ radical, formed in high yields upon irradiation of nonsensitized DNA). Indeed, in recent studies concerning the stability of the C5′ radical, cyclopurine lesions are predominantly observed, in line with our results.31,32 Due to the high stability of the formed C5′−C8 bond, cyclopurine lesions can accumulate in DNA over a long time,57 also in agreement with our findings (highly negative ΔG of formation of cyclopurine, see Table 1 and Figure 3). In addition, experimental data obtained for 6420
DOI: 10.1021/acs.jctc.7b00978 J. Chem. Theory Comput. 2017, 13, 6415−6423
Article
Journal of Chemical Theory and Computation
Figure 5. Yields of all products of radical stabilization reactions, obtained for different effective termination rate constants. Yields were obtained by numerical integration of the corresponding kinetic equations. In calculation of the model parameters, the heights of free energy barriers were drawn from the normal distribution (μ = ΔG‡, σ = δ(ΔG‡)) to account for standard errors in free energy calculations, as described in detail in the Methods and Table S1. After 100 independent runs, standard deviations of log10 of species’ final populations were calculated to obtain the error estimate (marked by vertical bars).
amount of lesions (
