Excited-State Properties and Relaxation Pathways of Selenium-Substituted Guanine Nucleobase in Aqueous Solution and DNA Duplex.

The excited-state properties and relaxation mechanisms after light irradiation of 6-selenoguanine (6SeG) in water and in DNA have been investigated using a quantum mechanics/molecular mechanics (QM/MM) approach with the multistate complete active space second-order perturbation theory (MS-CASPT2) method. In both environments, the S1 1(nSeπ5*) and S2 1(πSeπ5*) states are predicted to be the spectroscopically dark and bright states, respectively. Two triplet states, T1 3(πSeπ5*) and T2 3(nSeπ5*), are found energetically below the S2 state. Extending the QM region to include the 6SeG-Cyt base pair slightly stabilizes the S2 state and destabilizes the S1, due to hydrogen-bonding interactions, but it does not affect the order of the states. The optimized minima, conical intersections, and singlet-triplet crossings are very similar in water and in DNA, so that the same general mechanism is found. Additionally, for each excited state geometry optimization in DNA, three kind of structures ("up", "down", and "central") are optimized which differ from each other by the orientation of the C═Se group with respect to the surrounding guanine and thymine nucleobases. After irradiation to the S2 state, 6SeG evolves to the S2 minimum, near to a S2/S1 conical intersection that allows for internal conversion to the S1 state. Linear interpolation in internal coordinates indicate that the "central" orientation is less favorable since extra energy is needed to surmount the high barrier in order to reach the S2/S1 conical intersection. From the S1 state, 6SeG can further decay to the T1 3(πSeπ5*) state via intersystem crossing, where it will be trapped due to the existence of a sizable energy barrier between the T1 minimum and the T1/S0 crossing point. Although this general S2 → T1 mechanism takes place in both media, the presence of DNA induces a steeper S2 potential energy surface, that it is expected to accelerate the S2 → S1 internal conversion.

Introduction

In the prebiotic age, the protection of the ozone layer was not efficient and the flux of UV radiation on Earth was much higher than it is nowadays, resulting in an hostile atmosphere capable of causing n class="Chemical">severe damage to DNA.[1] In these conditions, canonical purine (adenine and guanine) and pyrimidine (cytosine, thymine, and uracil) nucleobases emerged as very photostable compounds, able to preserve the genetic code from the deleterious effects of UV radiation. Photostability originates from an efficient excited-state decay that releases the excess of energy via ultrafast and efficient radiationless deactivation processes.[2−5] Today, it is well established that low-energy conical intersections that allow the system to return to the electronic ground state in a short time scale are responsible for the excited-state decay efficiency.[6−11]

Small chemical modifications can drastically alter the inherent photostability of nucleobases.[12] Particular attention has been devoted to n class="Chemical">thiobases analogues—in which oxygen atoms are replaced by sulfur—because, unlike their canonical counterparts, intersystem crossing results in high quantum yields of triplet states.[13−20] For instance, 6-thioguanine (6tG) has an intersystem crossing quantum yield of ca. 60%.[17,21] The radical nature of the triplet states render thiobases potential applications as photosensitizers for photodynamic therapy against a variety of ailments.[22] Examples are 2,4-dithiothymine,[23] 6-thio-2-deoxyguanine,[24] or 2,6-dithiopurine.[25] A natural extension in this direction are seleno-nucleobases, where selenium replaces oxygen or sulfur. Current studies[26−29] point out that Se-nucleobases can form stable RNA,[30] DNA duplex,[31] and G-quadruplex structures.[32] They exhibit advantageous red-shift absorption spectra in comparison with their thio-counterparts,[33,34] and most importantly, due to their heavier atom, a faster intersystem crossing has been observed compared to their thio-analogs—albeit with a shorter triplet state lifetime.[34]

Theoretical theory is particularly well suited to unravel the photophysical mechanisms responsible for the efficient electronic population of excited states.[35] However, all previous studies in sulfur- and n class="Chemical">selenium-substituted nucleobases[36−43] have systematically excluded the effect of the DNA environment due to the complexity involved, despite it could strongly influence the photophysical properties of its chromophores.[44]

In this paper, we employ hybrid Quantum Mechanics/Molecular Mechanics (QM/MM) techniques[45,46] to achieve a realistic description of the biomolecular environment by including the effect of solvation and that of the double helix DNA background on 6-selenoguanine (6n class="Chemical">SeG). According to previous theoretical studies in gas phase,[47] in 6SeG the S21(ππ*) state transfers its population to the S11(nπ*) state via a conical intersection, from which the triplet states are formed. More recently, three selenium-substituted uracils (2SeU, 4SeU, and 2,4SeU) have also been investigated in gas phase and similar excited-state relaxation pathways to triplet states are reported.[48,49] The primary goal of this contribution is to examine whether the main photophysical events of 6SeG are affected by the presence of water and/or biological media.

Computational Details

System Setup

The 6SeG molecule was studied both in n class="Chemical">water solution and in a DNA environment. In both cases, first classical molecular dynamics (MD) simulations were performed with the Amber16 package[50] in order to obtain initial structures to use in subsequent QM/MM calculations. Intra- and intermolecular parameters were extracted from the generalized Amber force field (GAFF),[51] while the force field parameters of the selenium atom were taken from the parametrization of the 2-selenouridine molecule embedded in aqueous solution.[52] The RESP (Restrained Electrostatic Potential) procedure[53] was performed to get the set of atomic charges at the B3LYP/cc-pVDZ level.[54−58] In the water simulations, the TIP3P model[59] was used with a system comprising one solute molecule surrounded by 5406 water molecules inside of a truncated octahedron simulation box of 12 Å (Figure a). As no constraint was applied during the simulation, a reduced time step of 0.5 fs was employed. Initially, the system was minimized with the steepest descent algorithm, followed by heating in the NVT ensemble for a total time of 100 ps, scaling the temperature from 0 to 300 K. Afterward, an equilibration regime was applied for 1 ns in the NPT ensemble in normal conditions of temperature and pressure (P = 1 bar and T = 300 K). The final production was carried out for 10 ns, from which the last snapshot was employed to investigate the photophysics of 6SeG in water with a QM/MM approach.

Figure 1

QM/MM system setup. (a) 6SeG in water solution. (b) 6SeG embedded in DNA environment. Two different QM regions are indicated: only 6SeG (QM1) or 6SeG paired with cytosine (QM2). All remaining atoms are included in the MM subsystem.

QM/MM system setup. (a) 6n class="Chemical">SeG in water solution. (b) 6SeG embedded in DNA environment. Two different QM regions are indicated: only 6SeG (QM1) or 6SeG paired with cytosine (QM2). All remaining atoms are included in the MM subsystem.

In order to include the DNA environment, a duplex with the sequence 5′-ATGGTGCAC-3′ and 3′-TACCACGTG-5′ was employed (Figure b). This n class="Chemical">sequence was first constructed by using a nucleic acid builder tool in Amber16 package as in previous experimental studies.[60] The DNA duplex system was solvated in a water box of 56 × 59 × 69 Å. DNA, waters, and counterions were described with default ff99SB force fields[61] and TIP3P model.[59] Then, 1000 MM minimization cycles with frozen DNA were carried out, followed by 2000 cycles without any geometric constraints. The minimized system was heated to 300 K in a 40 ps equilibration MD simulation, followed by 10 ns MD simulations. As in water solution, for the QM/MM study the final MD snapshot was used, in which the oxygen atom of the fourth G base was replaced with a selenium atom (5′-ATG-SeG-TGCAC-3′:3′-TAC-C-ACGTG-5′).

QM/MM Calculations

The QM/MM calculations were performed using either MOLCAS8.0[62,63] or OpenMOLCAS[64,65] interfaced with TINKER6.3.[66] In aqueous solution, the QM region contained the 6SeG molecule while all the n class="Chemical">water molecules were described with MM. In DNA, two different QM regions were constructed to explore effects of interbase hydrogen-bonding interaction, see Figure b. The first (QM1) includes only 6SeG while the second contains the 6SeG paired with a cytosine base (6SeG-Cyt). In both models, all remaining atoms from the DNA duplex, sugar, phosphates, water molecules, and counterions are left in the MM region.

Regardless of the environment (only water or solvated DNA) the QM region was described with the state-averaged complete active space n class="Chemical">self-consistent field (SA-CASSCF) method[67] for geometry optimizations and the multistate second-order perturbation theory (MS-CASPT2) approach[68,69] on top to refine single-point energies. A total of 4 singlets and 3 triplets were included in the state-average procedure. During the SA-CASSCF calculations, C, N, O, and H atoms were described with cc-pVDZ and the Se atom with the aug-cc-pVDZ basis set.[58,70] The active space consisted of 12 electrons distributed over the 10 orbitals; see Figure for the active orbitals in water and Figures S1 and S2 for the orbitals within the DNA environment. The MS-CASPT2 energy calculations were carried out employing larger atomic basis sets, specifically cc-pVTZ for C, N, O and H, and aug-cc-pVTZ for Se. In this case, the active space was augmented to 14 electrons in 12 orbitals with σ and σ* orbitals located on the C=Se moiety; see Figure .

Figure 2

Orbitals included in the active spaces used for the QM/MM calculations in water. In green are shown the orbitals of the (12,10) active space used in the optimizations. The two additional σ and σ* orbitals located on the C=Se moiety are only included to refine single point energies.

Orbitals included in the active spaces used for the QM/MM calculations in n class="Chemical">water. In green are shown the orbitals of the (12,10) active space used in the optimizations. The two additional σ and σ* orbitals located on the C=Se moiety are only included to refine single point energies.

An imaginary level shift of 0.2 au was employed to avoid the issues of intruder states.[71] The IPEA shift value was n class="Chemical">set to zero.[72,73] The Cholesky decomposition technique for two-electron integrals was employed to accelerate CASSCF and MS-CASP2 calculations.[74] Spin–orbit couplings were obtained at the MS-CASPT2 level within the atomic mean-field (AMFI) approximation.[75−77] The effective spin–orbit couplings are expressed as follows:in which Ψ and Ψ are the perturbatively modified electronic wave functions of the corresponding singlet and triplet states; H, H, and H are the x, y, and z components of the spin–orbit operator.

We would like to emphasize that only adiabatic states are considered throughout this work. In this state representation, population transfers among states of the same multiplicity are induced by nonadiabatic couplings, which become extremely large clon class="Chemical">se to conical intersections. Population transfer between different multiplicities is mediated by spin–orbit couplings, whose magnitude is independent of the energy gaps. Hence, in addition to the location of relevant minimum-energy singlet–triplet crossings, we also computed the spin–orbit coupling matrix elements at these crossings.

For comparison, in addition to the QM/MM calculations of 6SeG in the pren class="Chemical">sence of explicit water molecules, we also carried out computations with the polarizable continuum model (PCM) that considers solvent effects only implicitly.[78,79]

To save computational effort, only the last saved snapshot from the MD simulations was used as a starting point for the optimization of excited state geometries with QM/MM. However, in order to verify whether this approach suffices to predict deactivation pathways, we also performed single point vertical excitation energy calculations considering ten randomly selected snapshots taken from the classical simulation in both environments (see Tables S1–S3). The superposition of the selected snapshots shows that the nearest water molecules are well distributed around the 6SeG (see Figure S3), which means that the selected snapshots are representative and give a good representation of the solvent configurations. Since the calculated root-mean-square deviation is small for the two lowest singlet states, with the largest deviation of 0.07 eV—within the accepted error of the method—we conclude that the use of one snapshot for the following discussion of the relaxation pathways is reasonable.

Results and Discussion

Equilibrium Structures

Figure shows the optimized geometries of 6SeG in the electronic ground state (S0) in n class="Chemical">water solution (a) and in DNA (b), as obtained with the QM(CASSCF)/MM approach. In aqueous solution, 6SeG has a structure nearly planar with the two hydrogen atoms of the amino group out of the six-membered ring plane by ca. 26°. The C=Se bond length is computed to be 1.790 Å, slightly shorter than that predicted in gas phase,[47] 1.821 Å.

Figure 3

QM(CASSCF)/MM optimized structures of 6SeG in the electronic ground state, in aqueous solution (a), in DNA (b), and in a pair with cytosine (c). Selected bond lengths are in Å, and atomic numbering is given.

QM(CASSCF)/MM optimized structures of 6SeG in the electronic ground state, in aqueous solution (a), in DNA (b), and in a pair with n class="Chemical">cytosine (c). Selected bond lengths are in Å, and atomic numbering is given.

The S0 structure of 6SeG embedded in DNA (Figure b) and that within the explicit 6n class="Chemical">SeG-Cyt base pair (Figure c) are also essentially planar. Noticeable is that the C=Se bond lengths are 1.827 and 1.848 Å, respectively, much longer than the corresponding C=S bond length observed in thio-counterpart (∼1.650 Å).[16,38,80] As expected, the explicit consideration of three hydrogen bonds in 6SeG-Cyt influence the geometric parameters of the 6SeG moiety. Accordingly, the hydrogen-bonding interaction related to the Se and N10 atoms is much weaker than that related to the N11 atom, as suggested by their corresponding bond lengths. The two hydrogen bonds, Se···H20 and N19···H21, are 2.635 and 2.357 Å, whereas the O12···H22 bond length is 1.933 Å. Since the former two hydrogen bonds are much weaker, there are small effects on the geometric parameters around the C2–Se1 moiety. Likewise, the stronger O12···H3 hydrogen-bonding interaction affects the nearby bond lengths; for example, the C9–N11 and N8–C9 bond lengths are 1.365 and 1.288 Å in 6SeG, but 1.312 and 1.336 Å in 6SeG-Cyt.

Overall, the S0 structures of 6SeG in both water and DNA environments are similar except for the C9–N11 and C=Se bond lengths–likely due to the pairing cytosine.

Spectroscopic Properties

Table collects the vertical excitation energies of 6SeG computed in n class="Chemical">water solution and within DNA, without and with the explicit consideration of the paired cytosine (QM1 vs QM2). We shall first discuss the results in water solution. The first singlet excited state, S11(nSeπ5*), corresponds to a dark state predicted at 3.13 eV in implicit solution and 2.62 eV with explicit solvent. Relative to the S0 state, this state is due to a single electronic excitation from the nSe nonbonding orbital localized on the Se atom, to the π5* antibonding orbital. As mentioned above, there is a difference of around 0.5 eV between the results obtained with the two solvation models. This difference is expected since the degrees of freedom of the solvent are neglected in the implicit solvation model. As known in the literature,[82] lone pairs states are more affected by the explicit solvation, especially in polar solvent, due to the formation of solute–solvent hydrogen bonds.

Table 1

Vertical Excitation Energies (ΔE, eV) and Related Oscillator Strengths (f) of 6SeG in Water Computed at the MS-CASPT2/PCM and MS-CASPT2(14,12)/MM Levels of Theory and in DNA with QM1(MS-CASPT2(14,12))/MM and QM2(MS-CASPT2(14,12))/MM (in Parentheses)

	water				DNA
	MS-CASPT2/PCMa		QM(MS-CASPT2)/MMb		QM(MS-CASPT2)/MMc
states	ΔE	f	ΔE	f	ΔE	f
S1 1(nSeπ5*)	3.13	0.000	2.62	0.000	2.56 (2.79)	0.000 (0.006)
S2 1(πSeπ5*)	3.29	0.577	3.51	0.570	3.40 (3.29)	0.525 (0.517)
S3 1πSeπ6*	4.73	0.006	4.43	0.016
T1 3(πSeπ5*)	2.60		2.55		2.35 (2.48)
T2 3(nSeπ5*)	2.88		2.59		2.51 (2.80)

Structure optimized with MS-CASPT2(12,10).

Structure in water optimized with QM(CASSCF(12,10))/MM.

Structures in DNA optimized with QM(CASSCF(12,10))/MM: 6SeG (6SeG-Cyt).

The S21(πSeπ5*) state is computed at 3.29 eV in implicit n class="Chemical">water and at 3.51 eV with explicit solvent. This is a singly excited configuration from the πSe orbital localized on the Se atom to the π5* orbital, and it is the bright state (f = 0.577). It is encouraging that our calculation with QM(MS-CASPT2)/MM level is in excellent agreement with experimental data that reported an absorption maximum of 6SeG in aqueous buffer solution at 357 nm (3.47 eV).[34] Thus, this electronic transition is responsible of the first peak of the electronic absorption spectrum of 6SeG in water (see Figure S4).

The third singlet excited state, S31(πSeπ6*), is around 0.9 eV higher in energy above the S21(πn class="Chemical">Seπ5*) state at the QM(MS-CASPT2)/MM level of theory. The corresponding oscillator strength is 0.016—about 1 order of magnitude lower than that observed for the electronic transition to the S21(πSeπ5*) state, meaning that its contribution to the first absorption band will be small. Accordingly, it will not be considered hereafter.

We have also calculated the two lowest-lying triplet excited states. They have the same electronic transition characters as the singlet counterparts S11(nSeπ5*) and S21(πn class="Chemical">Seπ5*), i.e., T13(πSeπ5*) and T23(nSeπ5*). Their vertical excitation energies at the QM(MS-CASPT2)/MM level of theory are 2.55 and 2.59 eV, respectively, and are overestimated with the PCM model, particularly the T2. The next triplet state, which is described by a T33(πSeπ6*) transition, is situated around 1 eV above the S21(πSeπ5*) state; thus, they are not expected to be accessed from the bright state and are not considered further.

We now proceed to discuss the excitation energies of 6SeG in DNA. As can be seen, both QM1 and QM2 regions (i.e., 6SeG and 6SeG-Cyt) predict the same energetic order of the excited singlet and triplet states and wave function character as in solution. The S1 state is spectroscopically dark with nπ* character, and the S2 state is bright due to the ππ* character. In the S11(nSeπ5*) state the π5* orbital is delocalized all over including the C=Se double bond. In the S21(πSeπ5*) state, the π5* orbital is the same as that in the S11(nSeπ5*) state, while the πSe orbital is mainly composed of the p_z orbital of Se with some minor contribution from the C3–C7 bond in either 6SeG or 6SeG-Cyt in DNA. It is worth mentioning that the lowest-energy states of 6SeG-Cyt are basically described by local excited states on 6SeGua and Cyt, respectively, due to the high energy difference between the excited states computed for each of them separately, which are below about 3 eV for 6SeG and 5 eV for Cyt. This is confirmed by QM(MS-CASPT2)/MM calculated vertical excitation energies on 6SeG-Cyt: the S3 (4.64 eV) and S4 (4.77 eV) states are local excitations on cytosine and charge transfer states, respectively, which are on average 1.4 eV higher than S2 (see Table S6). Therefore, charge transfer states are placed at higher energetic regions and the state character of the 6SeG-Cyt S1 and S2 lowest states are preserved in relation to those observed for the QM1 scheme.

Curiously, the hydrogen-bonding interaction that appears when the n class="Chemical">cytosine is treated quantum mechanically does not change the characters of the S1 and S2 states as aforementioned but it considerably affects their energies. The QM1(MS-CASPT2)/MM vertical excitation energies of S1 and S2 states are predicted to be 2.56 and 3.40 eV and are increased and decreased to 2.79 and 3.29 eV, respectively, with QM2(MS-CASPT2)/MM. This means that the hydrogen-bonding interaction stabilizes the S21(πSeπ5*) state by 0.11 eV but destabilizes the S11(nSeπ5*) state by 0.23 eV. A similar situation is found in T13(πSeπ5*) and T23(nSeπ5*) states. Their vertical excitation energies are predicted as 2.35 and 2.51 eV with QM1 but 2.48 and 2.80 eV with QM2, respectively.

It is apparent that the QM(MS-CASPT2)/MM approach predicts similar vertical excitation energies regardless of water or DNA environment, except for the fact that the n class="Chemical">hydrogen-bonding interactions destabilizes the S1 and T2 states by 0.23 and 0.29 eV at the QM(MS-CASPT2)/MM level.

Excited-State Minima

Figure shows the lowest-lying minimum-energy structures optimized in aqueous solution using the QM(CASSCF)/MM method. The S1 minimum (S1-MIN) is located 2.44 eV adiabatically above the S0 minimum at the QM1(MS-CASPT2)/MM level of theory. It is characterized by a pyramidalization of the Se atom (by 24°) and an elongated C=n class="Chemical">Se bond length (by 0.12 Å) with respect to the S0 minimum. The S2-MIN is predicted at 3.34 eV above the S0 minimum. It retains the Se atom in the molecular plane and also shows a stretched C=Se bond (by 0.118 Å) in comparison to the S0 minimum, while the amino group is rotated making its two hydrogen atoms almost perpendicular to the molecular plane.

Figure 4

Minimum-energy structures of 6SeG in water solution optimized at the QM(CASSCF)/MM level of theory. Selected bond lengths are in Å.

Minimum-energy structures of 6SeG in n class="Chemical">water solution optimized at the QM(CASSCF)/MM level of theory. Selected bond lengths are in Å.

The T1-MIN and T2-MIN also have the Se atom pyramidalized by about 30°. The C=n class="Chemical">Se bond length is longer in the T2 minimum than in the T1 (1.941 vs 1.936 Å). These minima are calculated 2.38 and 2.45 eV above the S0 minimum. Experimentally, the triplet state lifetime in 6SeGua is reported to be 835 times shorter than that observed for 6tG in aqueous solution.[34] Earlier theoretical investigations in 6tG pointed out the existence of two minimum structures in the T1 state,[81] but only a single minimum is predicted for 6SeGua. Hence, we suggest that the absence of the second minimum on T1 PES could be one of the reasons for the shorter triplet state population.

For completeness, these minima have also been optimized with the PCM model. In general, the results are very similar to thon class="Chemical">se obtained with QM/MM but the pyramidalization on the Se atom becomes larger (S1-MIN, 35°; S2-MIN, 44°; T1-MIN, 42°; T2-MIN, 41°). The observed larger pyramidalization angles in relation to those computed with the QM/MM method can be explained by considering how the solute–solvent hydrogen bonds are modeled in each approach. In addition, the adiabatic excitation energies are also influenced, with the PCM model predicting values of 2.99, 3.18, 2.40, and 2.71 eV for the S1-MIN, S2-MIN, T1-MIN, and T2-MIN structures, respectively. Overall, in the particular case of 6SeG, both solvation models well agree about the relative energy order of the minima.

The excited-state minimum-energy structures have also been optimized in DNA at the QM1 and QM2(CASSCF)/MM levels of theory; see Figure . The corresponding S1, S2, T1, and T2 minima can be classified into three groups, differing by the orientation of the C=n class="Chemical">Se group–above (“up”), below (“down”) or within the molecular plane (“central”), denoted as U, D, and C, respectively. The Se atom leaves the molecular plane pointing toward the guanine in the “up” structure, and pointing toward the thymine in the “down” structure (see Figure S9). As they are energetically quasi-degenerated (Tables S3 and S4), only the “up” conformations (referred as S1-MIN-U, S2-MIN-U, T1-MIN-U, and T2-MIN-U) will be discussed.

Figure 5

Minimum-energy structures of 6SeG in the S2, S1, T2, and T1 states in DNA optimized at QM1(CASSCF)/MM (top) and QM2(CASSCF)/MM (down) levels of theory. In these structures, the C=Se group is twisted out of the molecular plane, which is referred to as “up” configuration and indicated by “U” in the structure labels (see “central” and “down” in Figures S5–S8). Selected bond lengths are in Å.

Minimum-energy structures of 6SeG in the S2, S1, T2, and T1 states in DNA optimized at QM1(CASSCF)/MM (top) and QM2(CASSCF)/MM (down) levels of theory. In then class="Chemical">se structures, the C=Se group is twisted out of the molecular plane, which is referred to as “up” configuration and indicated by “U” in the structure labels (see “central” and “down” in Figures S5–S8). Selected bond lengths are in Å.

Conspicuously, the QM1(CASSCF)/MM calculations predict a C=Se bond length elongated in all the structures, by more than 0.30 Å in S2-MIN-U and by about 0.15 Å in the other minima, with respect to the S0 minimum. In comparison, the other bond lengths change less than 0.03 Å. The calculations with QM2(CASSCF)/MM are consistent, with changes in the C=n class="Chemical">Se bond ranging from 0.15 to 0.32 Å and the other bond differing by less than 0.03 Å. This implies that the presence of the pairing cytosine hardly influences the geometries of the 6SeG moiety. Interestingly, the electronic excitations obtained with both QM regions are also very similar, as it has been also observed in other nonnatural base pairs in DNA.[83] The adiabatic excitation energies of the S1, S2, T1, and T2 are to be, 2.46 (2.58), 3.08 (2.86), 2.27 (2.29), and 2.36 (2.53) eV at the QM1(MS-CASPT2)/MM level (QM2(MS-CASPT2)/MM in parentheses). Structures and energies of the other conformations can be found in Figures S5–S8 and Tables S4 and S5.

Crossing Points

Crossing points among different PES are key to disentangling excited-state decay pathways. Unlike excited-state minima, the optimization of crossing points is much more difficult. Here, we have first optimized two-state crossing points at the QM(CASSCF)/MM level of theory and after analyzing the nature of the involved electronic states, a single-point energy calculation was done with QM(MS-CASPT2)/MM, verifying whether the two states keep degenerated. In some can class="Chemical">ses, some two-state minimum-energy crossing points resulted in three-state degenerate points.

In water, three minimum-energy crossing points have been found; n class="Chemical">see Figure . The first one corresponds to a crossing between the S2 and S1 states (S2/S1). This conical intersection is very similar to the S2-MIN; as a matter of fact, both structures are quasi-degenerate in energy (S2-MIN lies at 3.34 eV while S2/S1 is at 3.34/3.29 eV). For this reason, we computed numerical frequencies at the SA-CASSCF level of theory (the same level of theory employed for optimizations) in all environments, to make sure that it is a true minimum and not simply the lowest point on the S1/S2 intersection seam. Similar to the three-state crossing point found for 6SeG in gas phase,[47] as well as for other canonical and modified nucleobases either in gas phase or in solution,[84−87] the search for a conical intersection between the two low-lying T1 and T2 states, returned a three-state degenerate crossing point on which the S1 state is also degenerated. This S1/T1/T2 three-state crossing structure has the Se atom out of the molecular plane by ca. 28° and a C=Se bond length of 1.949 Å.

Figure 6

Minimum crossing points of 6SeG in solution obtained with QM1(CASSCF)/MM. Selected bond lengths are in Å.

Minimum crossing points of 6SeG in solution obtained with QM1(CASSCF)/MM. n class="Chemical">Selected bond lengths are in Å.

The last crossing point found in aqueous solution is the singlet/triplet crossing T1/S0, between the S0 and T1 PES. This structure shows a distinctive large pyramidalization at the Se atom (ca. 61°)—larger than in any of the other optimized critical points.

Figure shows the obtained crossing points of 6SeG in DNA; as with the local minima, only the “up” structures are shown, unless otherwin class="Chemical">se stated (see Figures S10 and S11 for other conformations). The two S2/S1 conical intersections (S2/S1-U) are very similar, except for C=Se bond and its orientation relative to the molecular plane. As expected, the presence of the hydrogen bond shortens the C=Se bond (from 2.263 to 2.212 Å); however, in comparison, changes in the other bond lengths can be disregarded. The energetics of the S2/S1-U geometries are also affected by the consideration of the base pair (see Table S5). Accordingly, while the S2 and S1 states are at 3.28 and 3.20 eV for QM1(MS-CASPT2)/MM, they are stabilized to 2.90 and 2.86 eV for QM2(MS-CASPT2)/MM; see Table S2.

Figure 7

Minimum crossing points of 6SeG in DNA obtained with QM1(CASSCF)/MM (top) and QM2(CASSCF)/MM (down) levels of theory. Selected bond lengths are in Å.

Minimum crossing points of 6SeG in DNA obtained with QM1(CASSCF)/MM (top) and QM2(CASSCF)/MM (down) levels of theory. n class="Chemical">Selected bond lengths are in Å.

Singlet/triplet T1/S0 crossing points were also optimized for 6SeG in DNA. The T1/S0-U shows a C=n class="Chemical">Se bond length much longer than that in the S0 minimum (2.004 (QM1)/2.015 (QM2) Å in the T1/S0-U versus 1.827 (QM1)/1.848 Å in S0-MIN). The S0 and T1 energies at T1S0-U were calculated to be 2.58 and 2.62 eV (QM1) and 2.43 and 2.40 eV (QM2), i.e., in contrast to the S2/S1-U conical intersection. Here, the presence of cytosine has a negligible effect on the energies.

In addition to the two-state intersection structures, two different three-state intersection structures were also identified for 6SeG in DNA: the S2/S1/T2-C and the T2/T1/S0-U. The S2/S1/T2-C crossing structure is characterized by a much longer C=Se bond length (1.950 and 1.926 Å) and a Se–N11–N4–N6 dihedral angle (180° and 175°), considering QM1 and QM2 regions, respectively. The most striking structural feature of S2/T2/S1-C, with respect to the other intersection structures, is that the NH2 group becomes nearly perpendicular to the molecular plane. Moreover, the C=Se bond is essentially in the molecular plane, unlike the more pyramidalized feature of S2/S1-U and T1/S0-U. The energies of the S2/S1/T2-C structure are 4.11/4.07/4.10 eV at the best QM2(MS-CASPT2)/MM level of theory, indicating that this crossing point is not accessible.

The T2/T1/S0 crossing points differ on the orientation of the C=Se, and accordingly, they are labeled as n class="CellLine">T2T1S0-U and T2T1S0-D. From them, only T2/T1/S0-U is shown in Figure (see Figure S10 for the -D conformation). As can be seen, T2/T1/S0-U presents a longer C=Se bond length (2.625(QM1) or 2.680(QM2) Å) than in the other crossing points. The estimated energies of the T2/T1/S0 states are respectively 3.06/2.96/3.03 eV or 3.04, 3.03, and 2.97 eV with the smaller QM1 region, showing that the hydrogen bonding interaction has very small effects on the energetics.

Excited-State Relaxation Pathways

Finally, having obtained all the important critical points of 6SeG in n class="Chemical">water and in DNA, we proceed now to discuss the excited-state relaxation mechanisms, assisted by the linear interpolation internal coordinates (LIIC) technique. Note that even when the energy barriers computed with the LIIC technique are typically overestimated, these scans allows us to postulate the most plausible deactivation mechanisms, which can be further confirmed only by means of a detailed nonadiabatic molecular dynamics study. Additionally, linear interpolation in internal coordinates scans allow us to verify the consistency of the active space employed for the computation of all points, whether or not there are energy barriers due to intermediate geometries connecting the initial and final points, and to enable a graphical presentation which resembles the familiar potential energy surface concept. Figure displays the obtained pathways for 6SeG in water.

Figure 8

Plausible relaxation pathway of 6SeG in water calculated at the QM(MS-CASPT2)/MM level of theory. Relative energies of minima and crossing points are in eV. Spin–orbit couplings at the singlet/triplet crossings are given in cm–1.

According to our results, the primary UV absorption event populates the bright S2 (1πSeπ5*) state. After that, the system evolves barrierless along the S2 PES toward the S2/S1 conical intersection, by dissipating 0.17 eV excess of energy. Such an energetic profile is also confirmed using the minimum-energy path technique with QM(CASSCF)/MM (see Figure S12). The S2/S1 conical intersection coincides with the S2-MIN that simultaneously is also a three-state intersection, where the S1, S2, and T2 electronic states are quasi-degenerated. Therefore, in the vicinity of S2-MIN, the S2 electronic population can be transferred to both S1 and T2 states, through a conical intersection or a singlet/triplet crossing, respectively. The latter process benefits from a large spin–orbit coupling of ca. 440 cm–1 between S2 and T2. However, which of these pathways is preferred should be elucidated with the help of dynamical simulations.

Presuming a partial population transfer from the S2 to S1 state, the system can then reach the S1 minimum (S1-MIN), located 0.90 eV below the S2/S1 conical intersection region, as suggested by the downhill LIIC pathway. Starting from the S1-MIN region, we investigate whether the S1/T2/T1 region, located only 0.17 eV above the S1-MIN region, could be accessed. The computed LIIC profile does not show any barrier along this path and taking into account that the SOC is also substantial (S1/T1 = 380 cm–1 and S1/T2 = 155 cm–1), we hypothesize that population transfer to the T1 state should be favorable. In principle, internal conversion to the ground state and the ISC process could occur in parallel. However, the optimized minimum energy crossing point between the S1 and the ground state was computed to be located 1 eV above the S1 min region. Based on that, the direct relaxation pathway from the S1 to the S0 state is suppressed, in agreement with what is observed in other similar systems.[16,48] Therefore, direct relaxation pathways to the S0 were not considered by us.

The effects of sulfur-to-n class="Chemical">selenium substitution in the photophysics promotes a substantial decrease of the triplet lifetime, as observed by the experiment.[34] Our LIIC scans shows that the ISC process is favorable for 6SeG in water; however, from static calculations, we are not able to quantify how much more efficient this process is in relation to 6tG. Thus, only nonadiabatic dynamics simulations could be used in order to better understand why the 6SeG triplet lifetime is shorter. From here, there is a barrierless pathway toward the T1-MIN (0.26 eV below the three-state S1/T2/T1 crossing region), from which radiative and nonradiative processes can occur. The system can emit phosphorescence at 1.63 eV. Alternatively, as shown by the LIIC pathway, the singlet/triplet T1/S0 crossing located adiabatically 0.26 eV above the T1 minimum, allows an intersystem crossing to the S0 state. The calculated spin–orbit coupling for this crossing (283 cm–1), despite smaller than the previous values, could be sufficient for an efficient deactivation. Although, with MS-CASPT2 the energy gap at the T1/S0 crossing is still 0.33 eV at the CASSCF optimized structure, we assume that population to the ground state is viable. For comparison, in gas phase, the T1/S0 crossing point is predicted 0.11 eV above the T1-MIN at the MS-CASPT2 level of theory[47] – a little lower than that in aqueous solution.

The excited-state deactivation mechanisms of 6SeG (QM1) and 6n class="Chemical">SeG-Cyt (QM2) in DNA have been also investigated, see Figure . In DNA, three possible deactivation mechanisms can be followed from the Franck–Condon region. One of them involves a barrierless relaxation to the S2 minimum with the up configuration i.e. S2-MIN-U, by releasing 0.32 eV of energy, as computed at the QM1(MS-CASPT2)/MM level (Figure a). Once in the S2-MIN-U region, the system evolves toward an internal conversion region with the S1 state (S2/S1-U, with the up configuration too, which is located adiabatically 0.20 eV above the S2-MIN-U region. From here, population is transferred to the S1 state, finally reaching barrierlessly the S1 minimum S1-MIN-U, in which the up configuration is kept. Analogously, compared to the 6SeG in water, the S1 of 6SeG and 6SeG-Cyt in DNA need to overcome 1.03 and 0.67 eV energy barriers to reach the S1/S0 conical intersection region, whereas LIIC paths connecting the S1 and triplet manifolds are nearly barrierless; consequently, we believe that S1 directly deactivating to the ground state is inefficient. The S1 minimum region also represents a three-state crossing region among the S1, T2, and T1 electronic states, on which the population can be transferred from the S1 to T1 or T2 states, via intersystem crossings (S1/T1 or S1/T2). The S1 → T1 intersystem crossing process is more likely due to its larger S1/T1 spin–orbit coupling of 439.4 cm–1, in accord to the El-Sayed rule.[91] Once on the T1 state, the system evolves to T1-MIN-U, from which the intersystem crossing with the S0 state, T1/S0 (T1/S0-U, SOC = 273 cm–1), can be reached after overcoming an energetic barrier of ca. 0.30 eV. Due to this barrier, the system will be trapped a while in the T1 state, before hopping to the S0 state.

Figure 9

Plausible relaxation pathway of 6SeG in DNA calculated at the (a) QM1(MS-CASPT2)/MM and (b) QM2(MS-CASPT2)/MM levels of theory. Relative energies of minima and crossing points are in eV. Spin–orbit couplings at the singlet/triplet crossings are given in cm–1.

Alternatively, the S2 state can also relax toward another S2 minimum, the S2-MIN-D (see Figure S13). In such a can class="Chemical">se, and as discussed previously, after reaching the T1-MIN-D, it is necessary to overcome an energy barrier of 0.26 eV to access the T1/S0-D crossing point. Finally, if instead, the S2 state relaxes to the S2-MIN-C minimum, its further relaxation to the S1 state via the S2/S1/T2-C three-state intersection point encounters a comparable barrier in the S2 state because this process involves a large rotation of the N11H2 group (see Figure S14). Accordingly, this pathway does not compete with the other two. However, if S2-MIN-C can be first converted either into S2-MIN-U or S2-MIN-D in the S2 state, the system initially populated at S2-MIN-C still can decay to the T1 state efficiently.

The calculations of the excited-state relaxation pathways of 6SeG in DNA using the GC ban class="Chemical">se pair (Figure b) are similar to those described above, with few exceptions. The S2 potential energy is much steeper, so that the initial S2 relaxation pathway toward its minimum, S2-MIN-U, should be faster. The internal conversion process from the S2 to S1 states via S2/S1-U is essentially barrierless (0.02 eV here versus 0.20 eV without C). This is mainly attributed to the hydrogen bond interaction between the 6SeG and C. Intermolecular hydrogen bonds stabilize the 6SeG ππ* excited state level, which significantly reduces the excited state lifetime. In fact, several similar phenomena were reported in the literature.[88−90] For example, Temps et al. observed that the formation of guanine and cytosine Watson–Crick base pairs results in ultrafast nonradiative transition.[88] Around the S1-MIN-U minimum, the system will undergo intersystem crossing to both T2 and T1 states. But, as mentioned above, the S1 → T1 intersystem crossing is preferred because it has a larger S1/T1 spin–orbit coupling of 399 cm–1, justified by the El-Sayed rule.[91] Once in the T1 state, the barrier separating T1-MIN-U from T1/S0-U increases to 0.14 eV, i.e., 0.14 eV higher than that without C. For completeness, we mention that when the S2 system relaxes to its other minimum S2-MIN-D, similar excited-state relaxation pathways are identified (see Figure S15) and the corresponding pathways from S2-MIN-C are not studied here because of the high energies of the S2/T2/S1–C triplet crossing point, which makes it unlikely (see above).

Figure displays a schematic summary of our proposed relaxation mechanisms. In short, the environment surrounding the 6n class="Chemical">SeG greatly affects how fast the internal conversion from the bright state to the lowest singlet state takes place. A high energy barrier is computed in the gas phase,47 disappearing in water. Depending on which structure will be primarily accessed in DNA, the population transfer to the S1 PES may be faster or slower. Afterward, the population will be easily transferred to the triplet states via ISC, since the internal conversion to the S0 is less important.

Figure 10

Possible deactivation mechanisms of 6SeG in gas phase, water, and DNA.

Conclusions

A multiscale QM/MM strategy including the very accurate MS-CASPT2 method for the QM part has been employed to explore the excited-state properties and excited-state decay pathways of 6SeG in n class="Chemical">water and in a DNA environment, and compare these to those in gas phase. In the DNA environment, the effect of the pairing cytosine was also investigated by including it into the QM region. Spectroscopically, in all the three environments, 6SeG is characterized by the same two lowest excited singlet states, i.e. S11(nSeπ5*) and S21(πSeπ5*), corresponding to a dark and the spectroscopically bright excited states, respectively. Two triplet electronic states, i.e., T13(πSeπ5*) and T23(nSeπ5*), were found to be energetically lower than the bright S2 state. We find that the order of the singlet and triplet states and their characters does not depend on the environment; moreover, the vertical excitation energies in water and DNA are very similar to each other, as long as only the 6SeG is considered quantum mechanically. The inclusion of explicit hydrogen-bonding interactions stabilizes the S21(πSeπ5*) state and destabilizes the S11(nSeπ5*) state.

For comparison, the inclusion of solvent is also taken into account by implicit methods to find that PCM model overestimates both S11(nSeπ5*) and T23(n class="Gene">nSeπ5*) states, in relation to the values computed with the QM/MM method.

In water, we have found that the excited-state relaxation pathway corresponds to an efficient relaxation from the initially populated S21(πn class="Chemical">Seπ5*) state to the T13(πSeπ5*) state. After irradiation, the S2 state quickly relaxes from the Franck–Condon region to an S2 minimum, near which a S2/S1/T2 three-state intersection region induces both S2 →S1 internal conversion and S2 → T2 intersystem crossing processes. Subsequently, it is expected that both the S1 or T2 states decay further to the T1 state via either S1 → T1 intersystem crossing or T2 → T1 internal conversion. Both intersystem crossing processes benefit from substantial spin–orbit couplings (S1/T2 = 440 cm–1 and S1/T1 = 380 cm–1). Once in the T1 state, the system should be trapped for some time, due to the existence of a small energy barrier between the T1 minimum and the T1/S0 crossing point, estimated to be 0.26 eV.

The inclusion of DNA environment increases the complexity of the excited-state relaxation pathways of 6n class="Chemical">SeG due to the existence of different conformational structures, that depend on whether the C=Se group is oriented above (“up”), below (“down”) or within (“central”) the molecular plane. In general, it is possible to identify a preferred pathway, with resemblances to that in water. For example, starting from the up configuration, the excited S2 system first arrives at its S2 minimum, and then, a nearby S2/S1 conical intersection drives the system down to the S1 state. An S1 → T1 intersystem crossing process then occurs in the vicinity of the S1 minimum, and an energy barrier of 0.30 eV traps the system in the T1 state before it can go to the S0 state. A similar excited-state decay pathway is found from the down configuration, while the central configuration is predicted not to be efficient because of a large energy barrier in the S2 state. Which of these pathways is most probable and in which amounts could be predicted by dynamical simulations.

The inclusion of the pairing cytosine has little influence on the obtained relaxation mechanism, except that the S2 → S1 internal conversion procesn class="Chemical">ses is predicted to be faster and larger energy gaps between the T1 minima and the T1/S0 crossing points are found.

The present work provides an important step toward understanding the excited state properties of chemically modified nucleobases and its response to light, on the way to find better photosensitizers, for example for phototherapies.