Spin-projected QM/MM Free Energy Simulations for Oxidation Reaction of Guanine in B-DNA by Singlet Oxygen.

Guanine is the most susceptible base to oxidation damage induced by reactive oxygen species including singlet oxygen (1 O2 , 1 Δg ). We clarify whether the first step of guanine oxidation in B-DNA proceeds via either a zwitterionic or a diradical intermediate. The free energy profiles are calculated by means of a combined quantum mechanical and molecular mechanical (QM/MM) method coupled with the adaptive biasing force (ABF) method. To describe the open-shell electronic structure of 1 O2 correctly, the broken-symmetry spin-unrestricted density functional theory (BS-UDFT) with an approximate spin projection (AP) correction is applied to the QM region. We find that the effect of spin contamination on the activation and reaction free energies is up to ∼8 kcal mol-1 , which is too large to be neglected. The QM(AP-ULC-BLYP)/MM-based free energy calculations also reveal that the reaction proceeds through a diradical transition state, followed by a conversion to a zwitterionic intermediate. Our computed activation energy of 5.2 kcal mol-1 matches experimentally observed range (0∼6 kcal mol-1 ).

Introduction

Reactive oxygen species such as n class="Chemical">singlet oxygen (1O2, 1Δg), hydroxyl radical (OH.) and superoxide anion (O2 −.) contribute to oxidation reactions in biological systems.[ , ] Singlet oxygen, commonly abbreviated as 1O2, is the first excited electronic state of oxygen molecule, lying 22.4 kcal mol−1 in energy above the triplet ground state (3Σg +). 1O2 can be generated by photosensitized reactions or be obtained from chemical sources such as hydrogen peroxide and endoperoxides. It exhibits a large reactivity toward electron‐rich olefins and aromatic compounds. The target biomolecules including nucleic acids, unsaturated lipids, and amino acids mostly undergo pericyclic reactions. While these oxidation reactions give rise to several diseases such as porphyria and skin cancer, they are useful in the natural product synthesis. In the present work, we focus on 1O2‐induced oxidative DNA damage implicated in various biological processes..[ , ] In DNA, guanine is the most susceptible base to oxidation damage, and is mainly oxidized to 8‐oxo‐7,8‐dihydroguanine (8‐oxoG). According to extensive experimental and computational studies, it is verified that the formation of 8‐oxoG is followed by a transient guanine 4,8‐endoperoxide intermediate (EP) generated by the [4+2] cycloaddition of 1O2 to the imidazole ring of guanine.[ , , , , , , , , , , , , ] It is still unclear whether the reaction proceeds through either a concerted mechanism or a stepwise mechanism through a zwitterionic (ZW) or a diradical (DR) intermediate (Scheme 1).

Scheme 1

Proposed two reaction mechanisms of an endoperoxide (EP) formation. Singlet oxygen (1O2) reacts with guanine to form either a concerted mechanism or a stepwise mechanism through a zwitterionic (ZW) or a diradical intermediate (DR).

Marchetti and Karsili studied the reaction mechanism of the formation of EP using an isolated system consisting of n class="Chemical">1O2 and guanine only (hereafter called QM‐only model). They performed multi‐configurational CASPT2/cc‐pVDZ calculations in the gas phase along an approximate reaction coordinate based on a linear interpolation in internal coordinates obtained at the spin‐restricted second‐order Møller‐Plesset (RMP2)/cc‐pVDZ level.[ , , ] Their results suggest that the [4+2] cycloaddition is virtually barrierless and proceeds via a concerted mechanism. Dumont et al. performed quantum mechanical and molecular mechanical molecular dynamics (QM/MM MD) simulations to investigate the reaction mechanism in aqueous solution and in B−DNA.[ , ] The spin‐restricted density functional theory (RDFT) method, namely the RLC−BLYP/6‐31+G* method, suggested a stepwise mechanism, in which the reaction starts with an attack of the proximal oxygen atom (Op) of 1O2 on the C8 atom of guanine to form ZW, and the following nucleophilic attack of the distal oxygen atom (Od) of the peroxo group on the C4 atom leads to EP. They concluded that the B−DNA environment stabilizes ZW and significantly enhances the formation of EP because the calculated barrier for the rate‐limiting process is small (6.2 kcal mol−1) and the overall reaction is strongly exergonic (∼−60 kcal mol−1).

However, it has been pointed out that energy profiles based on spin‐restricted calculations are in general not accurate for the open‐shell electronic structure of n class="Chemical">1O2, due to the inherent inability to describe static correlation effects arising from the double degeneracy of π* orbitals.[ , , , , , ] The broken‐symmetry spin‐unrestricted DFT (BS−UDFT) calculation can handle such static correlation, but its broken‐symmetry singlet state significantly suffers from the spin contamination from triplet ground state.[ , , ] Several methods have been proposed to remove spin contamination from electronic energies for broken‐symmetry singlet states.[ , , ] Thapa et al. investigated both the zwitterionic and diradical pathways using the QM‐only model. Yamaguchi's approximate spin projection (AP) method described by eq. (1) (see Computational Method) was applied to a range of DFT functionals including LC−BLYP. Solvent effects were implicitly included by means of SMD solvation model. The single‐point SMD/NEVPT2/6‐31+G** calculations on the SMD/BS‐UB3LYP‐optimized stationary points were also performed for comparison. They found that inclusion of implicit solvent effects can stabilize ZW, and the stepwise zwitterionic pathway was found to be more preferable. It was also indicated that the effect of spin contamination is present in transition states and intermediates. The activation enthalpy relative to isolated reactants predicted by AP−ULC−BLYP (8.7 kcal mol−1) was comparable to that by NEVPT2 (7.1 kcal mol−1), assuring the applicability of the AP method to the investigation of guanine oxidation as well as other singlet oxygen reactions.

Motivated by these computational results, we wish to provide a comprehensive description of guanine oxidation in B−DNA by applying the AP method to QM/n class="Disease">MM MD simulation. The adaptive biasing force (ABF) method was employed to determine free energy landscape. In light of the previous studies mentioned above and computational efficiency, the LC−BLYP/6‐31G* model was chosen to treat the QM region (Figure 1a).[ , , , , , ]

Figure 1

a) QM/MM model with bulk water molecules (left) and QM region (right) used in this study. b) Final snapshots (DR, ZW, EP) from the QM/MM MD equilibration exploited as the initial configurations for ABF simulations, with isovalue surfaces of positive (orange) and negative (yellow) spin density (isovalue=0.01).The RLC−BLYP/6‐31G* method was used for ZW and EP, and the BS−ULC−BLYP/6‐31G* method was used for DR.

a) QM/MM model with bulk water molecules (left) and QM region (right) used in this study. b) Final snapshots (DR, ZW, n class="Gene">EP) from the QM/MM MD equilibration exploited as the initial configurations for ABF simulations, with isovalue surfaces of positive (orange) and negative (yellow) spin density (isovalue=0.01).The RLC−BLYP/6‐31G* method was used for ZW and EP, and the BS−ULC−BLYP/6‐31G* method was used for DR.

The present approach allows to make a more accurate and unbiased analysis, retaining the benefits of accounting for both the B−DNA environment and static correlation. With regard to singlet oxygen reactions, the use of AP−UDFT is advantageous over more rigorous multi‐configurational calculations and much faster semi‐empirical calculations.[ , ] Multi‐configurational calculations are prohibitively expensive for QM/MM free energy calculations, since several tens of thousands of electronic structure calculations are needed to persist several tens of picoseconds. We notice that semi‐empirical methods cannot describe the n class="Chemical">1O2 reduction process appropriately, even if a parameter set tuned to open‐shell systems is utilized. To illustrate the applicability of the spin projection method to QM/MM MD simulations, the performance of AP−ULC−BLYP is compared with that of RLC−BLYP and BS−ULC−BLYP. We demonstrate how inclusion of static correlation and removal of spin contamination improve the accuracy of free energies. Emphasis is also placed on whether the reaction proceeds via either ZW or DR. Although several QM/MM MD studies have been reported for open‐shell systems such as metalloproteins,[ , , ] to the best of our knowledge, this is the first study that deals with a reaction on the open‐shell singlet surface using spin projection.

Results and Discussion

Singlet‐Triplet Energy Gap for O2

The adiabatic singlet‐triplet energy gap (ΔE ST) for 1O2 has often been examined to check the accuracy of theoretical methods for describing of its electronic structure.[ , , , , ] We show how the electronic energies of the singlet state vary between n class="Gene">RLC−BLYP, BS−ULC−BLYP, and AP−ULC−BLYP. The 6‐31G* basis set was used. The energy of the reference triplet ground state was calculated at the BS−ULC−BLYP/6‐31G* level, with an assumption that the state is free from spin contamination. The ΔE ST values obtained using RLC−BLYP, BS−ULC−BLYP, and AP−ULC−BLYP are 39.5, 11.0, and 22.0 kcal mol−1, respectively. It is seen that the RLC−BLYP calculation is unsatisfactory because of static correlation error. The resulting ΔE ST value of 39.5 kcal mol−1 significantly overestimates the experimental value of 22.4 kcal mol−1. The singlet state calculated by BS−ULC−BLYP has an BS value of 1.0029, which is typical for diradical species containing an equal mixture of singlet and triplet. The severe spin contamination by the triplet ground state causes an underestimation of the ΔE ST value (11.0 vs. 22.4 kcal mol−1). The AP−ULC−BLYP calculation successfully corrects the energy of the singlet state, leading to a good agreement with the experimental data (22.0 vs. 22.4 kcal mol−1). On the basis of the results, we wish to propose a hypothesis that the three different calculations will give rise to different free energy profiles for singlet oxygen reactions. The RLC−BLYP calculation is likely to underestimate activation barriers and overestimate exergonicity of a target reaction, whereas the opposite is true for the BS−ULC−BLYP calculation. The AP−ULC−BLYP is expected to lie between the two calculations and to provide a much better description.

Diradical Pathway

The QM/MM ABF simulations using BS−ULC−BLYP and AP−ULC−BLYP for the QM region were started with the DR intermediate shown in Figure 1. We will hereafter refer to, for example, QM(BS−ULC−BLYP)/MM ABF as BS−ULC−BLYP for convenience. The electronic structure of DR can be viewed as a perfect diradical with an BS value of 1.0215 calculated at the BS−ULC−BLYP/6‐31G* level. The reaction coordinate was described by the C8‐Op distance. The one‐dimensional potentials of mean force (1D‐PMFs) for the diradical pathway calculated at the BS−ULC−BLYP and AP−ULC−BLYP levels are shown in Figure 2. Also displayed in Figure 2 are representative snapshots of the reactant state (RS), transition state (TS1), and DR intermediate state (DR) and corresponding n class="Gene">spin density populations obtained with AP−ULC−BLYP. Analogous results are found in their BS−ULC−BLYP counterparts (RS, TS1, DR). The obtained RS has a C8⋅⋅⋅Op separation of ∼3.0 Å and the Op‐Od bond length of 1.19±0.03 Å. In TS1, the C8‐Op bond‐forming distance is shortened to 1.95±0.03 Å, while the Op‐Od bond is slightly elongated to 1.23±0.03 Å. Compared with these values in TS1 (1.90±0.02 and 1.26±0.03 Å), TS1 can be indicated as a more reactant‐like transition state. Upon the syn‐addition of 1O2, the DR intermediate is generated, with C8‐Op and Op‐Od bond distances of 1.43±0.03 Å and 1.30±0.03 Å, respectively. These critical structural parameters are comparable to those in DR used as the initial structure. All computed BS values are ∼1.0 throughout the simulation, confirming that the reaction follows the diradical pathway. Spin densities on the 1O2 and guanine moieties also characterize the diradical nature of the three states.

Figure 2

Top: 1D‐PMFs for the diradical pathway obtained with the BS−ULC−BLYP/MM ABF and AP−ULC−BLYP/MM ABF simulations. Bottom: Representative snapshots of RS, TS1, and DR calculated by the AP−ULC−BLYP/MM ABF simulations, with isovalue surfaces of positive (orange) and negative (yellow) spin density are plotted (isovalue=0.01).

Top: 1D‐PMFs for the diradical pathway obtained with the BS−ULC−BLYP/MM ABF and AP−ULC−BLYP/MM ABF simulations. Bottom: Representative snapshots of RS, TS1, and DR calculated by the AP−ULC−BLYP/MM ABF simulations, with isovalue surfaces of positive (orange) and negative (yellow) n class="Gene">spin density are plotted (isovalue=0.01).

The superior performance of AP−ULC−BLYP over BS−ULC−BLYP is found in the estimation of relative free energies. The activation and reaction free energies computed with BS−ULC−BLYP are 12.3 and 2.5 kcal mol−1, whereas those with AP−ULC−BLYP are decreased to 7.1 and −5.6 kcal mol−1. Clearly, the C8‐Op bond formation is found to be more facile and the exergonicity of the reaction is reproduced by eliminating the large effects of n class="Gene">spin contamination.

Zwitterionic Pathway

The QM/MM‐based ABF simulations for the zwitterionic pathway were started with the ZW intermediate shown in Figure 1, which has a closed‐shell electronic structure ( BS=0). The C8‐Op distance was used for the reaction coordinate in line with the diradical pathway. Figure 3 depicts the 1D‐PMFs corresponding to the zwitterionic pathway obtained with the QM/MM ABF simulations (QM=n class="Gene">RLC−BLYP, BS−ULC−BLYP, and AP−ULC−BLYP).

Figure 3

Top: 1D‐PMFs for the zwitterionic pathway obtained with the RLC−BLYP/MM ABF, BS−ULC−BLYP/MM ABF, and AP−ULC−BLYP/MM ABF simulations. Bottom: Representative snapshots of RS, TS2, R/U, and ZW calculated by the AP−ULC−BLYP/MM ABF simulations, with isovalue surfaces of positive (orange) and negative (yellow) spin density (isovalue=0.01).

Top: 1D‐PMFs for the zwitterionic pathway obtained with the RLC−BLYP/MM ABF, BS−ULC−BLYP/MM ABF, and AP−ULC−BLYP/MM ABF simulations. Bottom: Rn class="Gene">epresentative snapshots of RS, TS2, R/U, and ZW calculated by the AP−ULC−BLYP/MM ABF simulations, with isovalue surfaces of positive (orange) and negative (yellow) spin density (isovalue=0.01).

The RLC−BLYP method provided the reaction composed of RS, TS2, and ZW. In RS, the C8⋅⋅⋅Op sn class="Gene">eparation is 2.58±0.08 Å, which is comparable to the length of 2.66±0.03 Å reported in a previous QM/MM MD study. These values are shorter than those computed with BS−ULC−BLYP and AP−ULC−BLYP for the diradical pathway (∼3.0 Å), but the Op‐Od bond length of 1.20±0.02 Å seems to support that the interaction between 1O2 and guanine is negligibly small, The C8⋅⋅⋅Op distance of 2.06±0.08 Å in TS2 is slightly longer than that for the anti‐addition of 1O2 in the B−DNA environment (∼1.95 Å) reported by Dumont and co‐workers. The C8‐Op and Op‐Od bond distances of 1.31±0.02 Å and 1.46±0.03 Å in ZW are compatible with the initial ZW structure (Figure 1). The activation and reaction free energies are calculated to be 2.7 and −26.7 kcal mol−1, while those for the anti‐addition are 6.2 and −25.0 kcal mol−1. It indicates that syn‐addition is more favorable over anti‐addition in the B−DNA environment. The result agrees well with a previous QM‐only model calculation that employed a range of theoretical methods and showed that the syn‐addition is favored by 4∼6 kcal mol−1 over the anti‐addition judging from the activation barrier.

The BS−ULC−BLYP method using the Guess=Mix keyword (see Computational Method) gave only closed‐shell solutions before the onset of the restricted/unrestricted instability, due to the closed‐shell nature of ZW served as the starting point. As the sampling progressed, however, open‐shell electronic structures having non‐zero BS values gradually appeared. Such symmetry breaking led to a lowering of the potential energy due to spin contamination from the triplet ground state. With regard to the simulations using AP−ULC−BLYP, calculations for the triplet state were also conducted to compute the n class="Gene">spin‐projected energies in regions beyond the onset of the restricted/unrestricted instability. Then, every BS−ULC−BLYP/MM calculation yielded a broken‐symmetry solution with complete diradical character ( BS∼1.0) in the vicinity of a region with a C8⋅⋅⋅Op distance of ∼1.6 Å (denoted as R/U and R/U, Figure 3). It may be because the MM region also underwent gradual geometric changes so that it can fit into the open‐shell electronic structure of the QM region. This behavior caused the resampling of the already sampled points, but no minimum corresponding to DR started to appear. It means that the ZW intermediate is much more stable than the DR intermediate, being in good agreement with previous studies using QM calculations.

Consequently, 1D‐PMFs obtained using BS−ULC−BLYP and AP−ULC−BLYP methods resulted in having only one minimum corresponding to ZW and ZW. On the basis of their QM‐only model calculations, Thapa et al. indicated that a zwitterionic intermediate (virtually identical to ZW) can be formed directly via a transition state having a weak diradical character with BS=0.34 and the C8⋅⋅⋅Op distance of 1.87 Å.[ , ]

Such behavior is not observed in our computed results. We find that the nature of electronic structure changed from open‐shell to closed‐shell in R/U and R/U after the diradical transition states (TS2, TS2). The C8⋅⋅⋅Op distances in TS2 and TS2 (1.98±0.05 and 1.98±0.04 Å) are longer than that for the QM‐only model. From this comparison, it can be argued that the explicit B−DNA environment is likely to stabilize the open‐shell electronic structure and to lead to an earlier transition state. The simulations using BS−ULC−BLYP can reproduce the free energy profile qualitatively. Both the relative energies in TS2 (9.3 kcal mol−1) and ZW (−18.4 kcal mol−1) are lower than those observed in TS1 (12.3 kcal mol−1) and DR (2.5 kcal mol−1) for the diradical pathway. However, one must ken class="Gene">ep in mind that these values contain spin contamination. The computed activation barrier of 9.3 kcal mol−1 exceeds the experimentally observed range of 0∼6 kcal mol−1 for reactions of singlet oxygen.

The AP−ULC−BLYP‐based simulations can remove the spin contamination effect, thereby decreasing the activation and reaction free energies to 5.2 and −22.0 kcal mol−1. Unlike BS−ULC−BLYP, the activation free energy of 5.2 kcal mol−1 rn class="Gene">eproduces well the experimentally observed range (0∼6 kcal mol−1). Compared with the results for the diradical pathway, the zwitterionic pathway is found to be more exergonic, with a smaller activation barrier. As such, the zwitterionic pathway is expected to occur.

Evaluation of Spin Contamination and Static Correlation

Comparison is made between the free energy profiles computed with RLC−BLYP, BS−ULC−BLYP, and AP−ULC−BLYP methods, to quantitatively evaluate the static correlation and n class="Gene">spin contamination errors. The differences between the three QM methods in the ΔE ST value for 1O2 are reflected in the calculated activation and reaction free energies, supporting the aforementioned hypothesis. For the diradical pathway, spin contamination errors (BS−ULC−BLYP vs. AP−ULC−BLYP) on the activation and reaction free energies are calculated to be 5.2 and 8.1 kcal mol−1, respectively. To characterize the spin contamination effect along the reaction coordinate, vertical ΔE ST values in the three states were additionally evaluated on the AP−ULC−BLYP geometries taken from the snapshots depicted in Figure 2. The values calculated by BS−ULC−BLYP (AP−ULC−BLYP) are 10.9 (21.9), 5.8 (11.7), and 0.4 (0.8) kcal mol−1 for RS, TS1, and DR, respectively (see also Table S1 in the Supporting Information). It underlines that the magnitude of spin contamination decreases in the order RS>TS1>DR. The AP−ULC−BLYP method is needed to be used for sampling configurations close to RS and TS1, where the potential energies differ significantly with or without AP correction. Otherwise, the reaction is predicted to be endergonic, which is qualitatively inconsistent with previous experimental and computational findings.

Concerning the zwitterionic pathway, the spin contamination errors (BS−ULC−BLYP vs. AP−ULC−BLYP) on the activation and reaction free energies turn out to be 4.1 and 3.6 kcal mol−1. These values are large, but less prominent than those for the diradical pathway. This is presumably due to the closed‐shell nature and remarkable stability of ZW. As most of the sampled configurations were in the vicinity of ZW, the n class="Gene">spin contamination effects in RS and TS2 may not significantly affect the free energy landscape. A similar trend is observed in the static correlation errors (RLC−BLYP vs. AP−ULC−BLYP) on the activation and reaction free energies (i. e., 2.5 and 4.7 kcal mol−1). The error on the activation barrier is smaller as compared to that arising from spin contamination. The present analysis apparently provides quantitative evidence that the 1D‐PMF computed with RLC−BLYP can show better agreement with the experimental data, as suggested in previous computational studies.[ , ] Although one might argue that RLC−BLYP is computationally efficient for studying guanine oxidation, the small activation barrier stems from the fact that RS does not correspond to the correct energetic reference point (details in Supporting Information). As such, AP−ULC−BLYP is more reliable in that it is capable of characterizing the energetically correct reference state (RS) and diradical TS2, and of providing the smooth conversion from the open‐shell surface to the closed‐shell one.

Ring‐closure Reaction

Then, we explored the subsequent intramolecular ring‐closure reaction indn class="Gene">ependently starting from EP presented in Figure 1. Since no open‐shell species is considered to be involved in this process, the QM subsystem was treated by the RLC−BLYP/6‐31G* method. The RLC−BLYP/MM‐based ABF simulations were carried out using the reaction coordinate described by the C4‐Od distance. Figure 4 shows the obtained 1D‐PMF. The obtained reactant state is essentially identical to ZW. We call the state ZW hereafter. Representative snapshots of the reactant state, transition state, and intermediate state (denoted as ZW, TS3, and EP) are also depicted.

Figure 4

Top: 1D‐PMFs for the ring closure reaction obtained with the RLC−BLYP/MM ABF simulations. Bottom: Representative snapshots of ZW, TS3 and EP.

Top: 1D‐PMFs for the ring closure reaction obtained with the RLC−BLYP/MM ABF simulations. Bottom: Rn class="Gene">epresentative snapshots of ZW, TS3 and EP.

The C4‐Od bond formation requires an activation barrier of 4.2 kcal mol−1 in TS3, where the C4⋅⋅⋅Od is shortened to 3.03±0.07 Å from 3.53±0.04 Å in ZW. Similarly, the Op‐Od distance is decreased to 1.43±0.03 Å from 1.45±0.04 Å, whereas the O8‐Od is stretched to 1.36±0.03 Å from 1.31±0.02 Å. A previous QM/MM MD study that used two C8‐Op and C4‐Od bonds successively as reaction coordinates predicted that the C4‐Od bond formation stn class="Gene">ep is almost barrier‐free with a large reaction free energy of −37.0 kcal mol−1. Meanwhile, our results suggest that formation of ER from ZW is slightly facile than the C8‐Op bond formation in terms of the smaller activation free energy (4.2 vs. 5.2 kcal mol−1) and the resulting EP lies −15.4 kcal mol−1 below ZW. Nevertheless, our findings are consistent with the notion that the reaction proceeds via a stepwise mechanism, the first C8‐Op bond formation is the rate‐limiting step, and the overall reaction is strongly exergonic (∼−37.4 kcal mol−1).

Conclusions

In this study, cycloaddition reaction of 1O2 with n class="Chemical">guanine in B−DNA is explored by means of QM/MM ABF calculations, to assess whether the reaction proceeds via either a zwitterionic or a diradical intermediate. We demonstrate how inclusion of static correlation via the broken‐symmetry spin‐unrestricted approach and removal of spin contamination by Yamaguchi's AP method improve the accuracy of free energies. Our simulations highlight that the free energy profile vary depending on correlation treatment in the QM region. The AP−ULC−BLYP calculation gives a realistic free energy profile by accounting for the effects of static correlation and spin contamination. It turns out that the static correlation errors on the activation and reaction free energies are 2.5 and 4.7 kcal mol−1 for the zwitterionic pathway, while spin contamination errors on these quantities range from 3.6 to 8.1 kcal mol−1. These errors are sufficiently large not to be neglected. We emphasize that the AP method is, in particular, essential for investigating the diradical pathway, where the spin contamination effects are pronounced. According to the free energy profiles obtained with the accurate and robust QM(AP−ULC−BLYP)/MM ABF simulations, the zwitterionic pathway is expected to occur. Our computed results match the experimental findings, with the notion that the reaction proceeds via a stepwise mechanism, the first C8‐Op bond formation is the rate‐limiting step, and the overall reaction is strongly exergonic. The computed activation free energy of 5.2 kcal mol−1 shows agreement with the experimentally observed range (0∼6 kcal mol−1).

Computational Method

The X‐ray crystallographic structure of the B−DNA dodecamer d(CGCGAATTCGCG)2 was taken from the Protein Data Bank (PDB code: 1BNA, resolution 1.9 Å). We constructed three QM/MM models (DR, ZW, and EP), each of which has the same QM region consisting of the second n class="Chemical">guanine from the 5’‐end and 1O2 (17 atoms). Hydrogen atoms were added by using the CHARMM‐GUI input generator. The N and C atoms in the N‐glycosidic bond were set to the boundary between the QM and MM subsystems for both models as presented in Figure 1. The initial structures of the QM regions for ZW and EP were optimized at the RLC−BLYP/6‐31G* level of theory, while the BS−ULC−BLYP/6‐31G* model was used to optimize that for DR. Geometry optimizations for the QM subsystem were carried out by using the Gaussian 09 program package. The systems were solvated within a 70 Å cubic box of TIP3P water molecules and neutralized by 22 Na+ ions using the Visual Molecular Dynamics (VMD) program. Then, classical 1 ns MD simulations were performed for each model in the NPT ensemble at 300 K with a time step of 2.0 fs by using NAMD. A cut‐off of 14 Å was employed with a switching for nonbonded interactions, and electrostatic interactions were treated by the particle mesh Ewald method. During the MD simulations all the atoms in the QM region represented by the CHARMM General Force Field were kept fixed. The CHARMM36 force field and TIP3P water models were exploited to describe the rest of the system.

All QM/MM MD simulations were conducted with NAMD. The QM sn class="Chemical">ubsystem was described by the Gaussian 09 program package, whereas the MM subsystem was treated by the same parameters as in the classical MD simulations mentioned above. The energies of the open‐shell singlet state were corrected using Yamaguchi's approximate spin projection (AP) method given by:

with

where and represent the total energies for the broken‐symmetry singlet and triplet states obtained with BS−ULC−BLYP. The and are the corresponding expectation values of total n class="Gene">spin angular momentum. We made a wrapper script that supports the interface to the Gaussian 09 program and adapts to the standardized NAMD QM/MM interface. The forces on the MM point charges (F MM) derived from the interaction with the QM atoms were calculated by means of the electric field (E MM) at point charges (q MM) as:

The Charge, Density, and Prop=(Read, Field) keywords were applied to calculate as with AMBER‐Gaussian and DL_POLY‐Gaussian interfaces.[ , ] The present NAMD‐Gaussian interface gives the same results as the native NAMD‐ORCA interface, provided that QM calculations are performed at the same level of theory (see Supporting Information including Table S2 for details). All QM/MM MD simulations were conducted in the NVT ensemble with a time stn class="Gene">ep of 0.5 fs. The link atom method was used with the Z3 scheme for manipulating MM point charges.

In this study, we explored three reactions; (i) formation of DR, (ii) formation of ZW, and (iii) formation of EP via ZW. First, the three systems (DR, ZW, and n class="Gene">EP) were equilibrated without any constraints for 2.0 ps, respectively. The RLC−BLYP/6‐31G* method was employed for the closed‐shell ZW and EP. The BS−ULC−BLYP/6‐31G* method was used to compute DR with an open‐shell character. Starting with the final snapshots obtained from the equilibration step (Figure 1), 1D‐PMFs for the reactions (i)–(iii) were computed, coupled with the ABF method implemented as a part of the collective variables module of NAMD. The RLC−BLYP, BS−ULC−BLYP, and AP−ULC−BLYP methods were exploited as the QM regions for the reaction (ii). The QM region for the reaction (i) was treated by BS−ULC‐LYP and AP−ULC−BLYP because the reaction takes place on the open‐shell surface, whereas RLC−BLYP was sufficient for the reaction (iii) occurring on the closed‐shell surface. The C8‐Op bond distance was chosen as the reaction coordinate for the reactions (i) and (ii), in a range from 1.36 to 3.10 Å, and from 1.28 to 3.10 Å, respectively. The reaction coordinate for the reactions (iii) was set to the C4‐Od bond distance ranging from 1.40 to 3.80 Å. For all cases, the bin width was set to 0.02 Å, and a threshold of 100 samples prior to the application of the bias was used. Concerning the reaction (ii), the restricted/unrestricted instability necessarily takes place, since the reactant 1O2 is a diradical and the ZW intermediate is closed‐shell. Thus, during the ABF simulation, the initial guess for every QM step was generated by using the Guess=Mix keyword (see also Supporting Information). ABF simulations were run for 20∼40 ps in total (i. e., 40,000∼80,000 QM calculations).

Conflict of interest

The authors declare no conflict of interest.

As a service to our authors and readers, this journal provides supporting information supplied by the authors. Such materials are peer reviewed and may be re‐organized for online delivery, but are not copy‐edited or typeset. Technical support issues arising from supporting information (other than missing files) should be addressed to the authors.

Supplementary

Click here for additional data file.