Intersystem Crossing Pathways in the Noncanonical Nucleobase 2-Thiouracil: A Time-Dependent Picture.

The deactivation mechanism after ultraviolet irradiation of 2-thiouracil has been investigated using nonadiabatic dynamics simulations at the MS-CASPT2 level of theory. It is found that after excitation the S2 quickly relaxes to S1, and from there intersystem crossing takes place to both T2 and T1 with a time constant of 400 fs and a triplet yield above 80%, in very good agreement with recent femtosecond experiments in solution. Both indirect S1 → T2 → T1 and direct S1 → T1 pathways contribute to intersystem crossing, with the former being predominant. The results contribute to the understanding of how some noncanonical nucleobases respond to harmful ultraviolet light, which could be relevant for prospective photochemotherapeutic applications.

Noncanonical nucleobases, i.e., DNA or RNA nucleobases with chemical modifications, are gaining increased attention as researchers are realizing that these modifications provide additional possibilities to encode genetic information.[1] One important chemical alteration is thionation, where the oxygen atom of a carbonyl group is substituted by a sulfur atom, yielding a thiobase. Such thiobases have been isolated from transfer-RNAs[2−4] from all domains of life,[5] but they are also commonly found in artificial DNA and RNA systems. They prevent mispairing in DNA-based nanotechnology;[6,7] are employed as fluorescent markers;[8−10] and—given their excellent binding affinity to biological targets—are also used clinically as cytotoxic agents,[11,12] as anti-inflammatory drugs,[13] or as antithyroid drugs.[14,15] Particularly due to their medicinal relevance in photochemotherapeutic applications,[16−21] unravelling the photophysical and photochemical properties of thiobases is decisive. The ability to produce high triplet yields is pivotal to generate singlet oxygen and other reactive species responsible for the cytotoxicity. This ability seems to be linked to the heavy-atom effect of the sulfur atom, because canonical nucleobases exhibit ultrafast (i.e., few picosecond) internal conversion back to the electronic ground state,[22,23] while thiobases do not.[24] Understanding the photochemistry of noncanonical nucleobases in general is also important in the context of molecular evolution in prebiotic settings, where only the most photostable molecules survived the harsh ultraviolet irradiation on the early earth.[22,25,26]

A number of recent studies focused on the excited-state dynamics of thiouracils, but a consensus regarding the mechanism to populate triplet states has not been reached. Time-resolved experimental studies show ultrafast near-unity triplet yield for 2-thiouracil (2TU), 2-thiothymine, 2-thiouridine, and 2,4-dithiouracil.[21,27−31] Stationary quantum-chemical calculations are available only for 2TU, but different relaxation pathways are proposed to explain intersystem crossing (ISC).[32−34] After populating the lowest bright state S2, all studies agree that internal conversion to the S1 should take place. As a subsequent deactivation pathway, two routes for ISC have been proposed: one direct S1 → T1 and involving mostly planar geometries[32,33] and one indirect S1 → T2 → T1 with geometries showing a large degree of pyramidalization.[34] Because stationary calculations cannot predict the relative importance of these routes and might even overlook additional pathways, a time-dependent study is mandatory.

Here, we investigate the dynamics of the excited states of 2TU (Figure ), with the aim of providing a clear-cut deactivation mechanism that can also help to enlighten the general photophysics of thiobases and other noncanonical nucleobases, which is a timely topic.[21,24,30,35−42] Accurate CASPT2-based surface hopping molecular dynamics simulations including singlet and triplet states have been performed, employing a local development version of the SHARC (Surface Hopping including ARbitrary Couplings) software suite.[43−45] Ab initio methods regularly employed for excited-state dynamics simulations, such as CASSCF or TD-DFT, often give a qualitatively inaccurate description of the potential energy landscape.[46−50] Hence, here we rely on the MS-CASPT2(12,9)/cc-pVDZ level of theory (multistate complete active space perturbation theory) for all underlying on-the-fly computations of energies, numerical gradients, and relevant couplings between the electronic states. While this is not the first application of MS-CASPT2 in excited-state dynamics simulations (pioneered by the group of Martínez[51]), here we employ the largest MS-CASPT2 setup to date, to the best of our knowledge. The usage of the MS variant of CASPT2 for dynamics is necessary, because the single-state variant fails to correctly describe excited-state crossings.[51,52] Although MS-CASPT2 can create state-oversplitting artifacts,[53,54] the good agreement between single-state and multistate CASPT2 static calculations for 2TU[34] suggests that these artifacts are not present here. As will be shown below, this computationally expensive approach achieves very good agreement with the available experimental ISC time scales.[21,30,42] The orbitals included in the active space (Figure S1) and other details of the calculations as well as a discussion regarding the quality and stability of the computations can be found in the Supporting Information.

Figure 1

Structure of 2-thiouracil and atom numbering.

To provide initial conditions for the dynamics simulations, first an absorption spectrum (see Figure S2) was computed from 2000 geometries sampled from the ground-state equilibrium Wigner distribution. The single absorption band obtained shows a maximum at 4.2 eV, almost exclusively due to excitation to the S2 state,[34] which is a ππ* excitation from the sulfur atom to one π* antibonding orbital of the pyrimidine ring, denoted as πSπ6* henceforth (see orbital labels in Figure S1). The lower-lying S1 state is of nSπ2* character and does not contribute much to the absorption. Higher-lying excited states were not computed for the absorption spectrum because we are focused on the dynamics after excitation to S2. Accordingly, the dynamics simulations are initiated in the bright S2 state and consider the interaction with all states below, i.e., the S1 and S0 as well as the three lowest triplet states, T1, T2, and T3 (3πSπ2*, 3nSπ2*, 3π5π6* at the Franck–Condon geometry). A total of 44 trajectories with a simulation time of 1 ps have been considered; the nuclear dynamics was simulated with a 0.5 fs time step, and the electronic wave function was propagated with a 0.02 fs time step using the “local diabatization” formalism.[55,56]

Figure a shows the time-evolution of the excited-state populations (thin lines) and the results of a global fit described below (thick lines). For all analysis purposes, the populations of T2 and T3 are merged, because the T3 population is only 2% on average. As can be seen, after excitation to the S2 state, internal conversion to the S1 state is basically completed within 200 fs. Simultaneously, ISC to all three triplet states proceeds on a time scale of a few hundred femtoseconds, whereas hardly any population of the ground state is observed within 1 ps. The triplet population at 1 ps is 80%, while most of the remaining population still resides in S1. This dynamical behavior is different from the canonical nucleobase uracil, where ISC is smaller and ground-state recovery is ultrafast.[22,23]

Figure 2

(a) Excited-state populations from 44 trajectories (thin lines) and global fits (thick lines) to the populations. (b) Net population transfer (numbers of trajectories) among the electronic states. (c) Fitted time constants (in femtoseconds) from the global fit. The T3 population was absorbed into the T2 population in all panels. The thickness of the lines in panels b and c relates to the amount of transferred population.

To identify the most relevant deactivation pathways and branching ratios, the net population transfer among the excited states is plotted in Figure b. Three different ISC routes are revealed: S2 → T2,3, S1 → T2,3, and S1 → T1. The S1 → T2 route (ignoring the T3 population henceforth), proposed in ref (34), accounts for about 40% of the ISC events and is thus the main ISC pathway. However, routes S1 → T1 and S2 → T2, suggested in earlier works,[32,33] account for 35% and 25% of ISC, respectively.

Time scales associated with the excited-state populations in Figure a could in principle be derived from a global fit to a kinetic model based on the transition graph in Figure b. However, the graph is relatively complex because it contains parallel pathways. In terms of chemical reaction network theory,[57] the stochiometry matrix corresponding to the transfer graph exhibits a two-dimensional nullspace. Intuitively, this means that several different combinations of rate constants can lead to the same overall population dynamics. Stated differently, for the given transition graph, the populations do not contain enough information to uniquely estimate all individual rate constants. Therefore, attempting to extract the individual transition rates for the three ISC routes by fitting will not be successful because the different ISC rate constants would be strongly correlated and large fitting errors and essentially arbitrary ISC rates would be obtained.

Instead, for the purpose of obtaining time scales, a simplified kinetic model was devised, where the S1 → T2 → T1 pathway is merged into the competing S1 → T1 pathway. A global fit based on this model (see the Supporting Information for more details) yields the time constants presented in Figure c and the corresponding fitted functions shown as thick lines in Figure a, which closely follow the actual populations. The error estimates given in Figure c were obtained with the bootstrapping method.[58] They describe mainly the uncertainty in the populations due to incomplete sampling (i.e., due to the small number of trajectories); systematic errors due to the level of theory and the semiclassical nature of the dynamics method are not included. Consequently, the errors serve as a check as to whether the number of trajectories is sufficient. According to the simplified kinetic model, internal conversion from S2 → S1 proceeds with a time constant of 59 fs. The S1 → S0 deactivation is a much slower process that is not well described here because there are too few hops to the ground state. In addition, two ISC time constants were obtained, one of 540 ± 100 fs for S1 → T1 and another of 250 ± 60 fs for S2 → T2. Both values can be combined to obtain an average ISC time constant of about 400 fs, which is in very good agreement with the values of 340–360 fs reported in aqueous and acetonitrile solution using transient absorption spectroscopy,[21,30] and also with the 330 fs obtained for 2-thiothymine.[28] Our ISC time constant also agrees reasonably well with the value of 775 fs from very recent multiphoton time-resolved photoelectron spectroscopy experiments,[42] considering that our simulations do not explicitly account for the multiphoton ionization probe. It has been shown that the inclusion of the time-dependent ionization probabilities in the simulations can easily change the time scales by a factor of 2.[59] The kinetic model can also be used to extrapolate the populations to times longer than 1 ps. Including the errors associated with the time constants, a total triplet yield of 90 ± 10% has been obtained, in agreement with the experimental value of 75 ± 20%[30] and also with the ones of 90 ± 10%[28] or 100 ± 5%[29] for the related 2-thiothymine.

With the aim of identifying important geometrical changes along the excited-state dynamics, the time-evolution of selected internal coordinates is displayed in Figure . After excitation in the Franck–Condon region, relaxation to two different S2 minima is possible according to stationary calculations.[32−34] These minima differ in structure and the occupation of the antibonding orbitals, either more localized at the C6 atom[32−34] or C2 atom[34] (recall Figure ), and therefore are labeled as 1πSπ2* and 1πSπ6* minima, respectively. Which of the minima is actually populated is easily solved by the dynamics simulations. The length of the C5=C6 bond length shown in Figure a undoubtedly reveals that most trajectories exhibit a bond length of 1.43 Å during the first 50 fs, which corresponds to the planar 1πSπ6* minimum. In contrast, only two trajectories moved into the pyramidalized minimum (with bond length 1.38 Å) of the 1πSπ2* state, being trapped there for several hundreds of femtoseconds, because there exist sizable barriers to leave this minimum. Accordingly, the S2/S1 (1πSπ2*/1nSπ2*) conical intersection presented in ref (34) is not operative. In contrast, because the 1πSπ6* minimum is located very close to another S2/S1 conical intersection involving the (1πSπ6*/1nSπ2*) states, it effectively funnels the population toward the S1 state, leading to an S2 lifetime of 59 fs. Subsequently, the bond length distribution shifts to smaller values, and for the rest of the simulation, the trajectories spend most of the time around 1.38 Å, which corresponds to the 1nSπ2* pyramidalized minimum in the S1 state (as well as to other minima, see below).

Figure 3

Time-dependent evolution of selected internal coordinates in the ensemble of 44 trajectories. The definition of the chosen internal coordinates is schematically depicted on the right (see also Supporting Information). Panel e shows the evolution of the dipole moment of the active state (contour plot) and the ground-state dipole moment μ(S0) (dashed line).

At the operative S2/S1 crossing, both the pyrimidine ring and the thiocarbonyl group are planar, but at the S1 minimum the S atom is displaced out of the ring plane. This strong pyramidalization at the C2 atom is characteristic of the 1nSπ2* minimum and can be quantified by the pyramidalization angles p7213 (defined as the angle between the S7=C2 bond and the C2–N1–N3 plane, see Figure b, and related to the planarity of the thiocarbonyl group) and p7513 (which describes how much the S atom is displaced from the molecular plane, see Figure c). Both pyramidalization angles are initially distributed around zero, but within 100 fs an equilibrium value of p7213 = ±40° is reached, whereas it takes more than 400 fs for p7513 to stabilize. Inspection of the trajectories (see also Figure S3 in the Supporting Information) confirms that the pyramidalization mechanism involves two stages. First, the system relaxes by puckering the C2 atom out of the ring plane, thereby stabilizing the energy of the 1nSπ2* state by about 0.35 eV. Because only the C2 atom needs to be moved, but not the heavy sulfur atom, puckering takes only tens of femtoseconds after transition to S1 (Figure b). In a second stage, the pyramidalizations at the nitrogen atoms N1 and N3 are removed and the planarity of the ring recovered while the pyramidalization at C2 is conserved. Because this second stage necessitates heavy-atom motion, it is far slower, taking several hundred femtoseconds (Figure c); this stage leads only to a marginal energetic stabilization (0.03 eV).

After relaxation to the pyramidalized S1 minimum, ISC proceeds efficiently, promoted by small energy gaps between the three states S1, T2, and T1 as well as sizable spin–orbit couplings (up to 150 cm–1).[34] On the basis of the trajectories, the “wagging” mode of the S atom has been identified as an important tuning mode that modulates the gap between the 1nSπ2* and 3πSπ2* states. This mode is quantified with the angle α721, sketched in Figure d. This angle affects the energy of the 3πSπ2* state more than that of 1nSπ2* and 3nSπ2* (Figure S4 in the Supporting Information). Thus, at angles of about 102° and 128°, the 3πSπ2* state crosses with the 1nSπ2* state, at energies only slightly above the minimum energy of the S1 minimum energy. At these two angles one can also see that the two triplet characters are exchanged, suggesting a nearby T2/T1 conical intersection that favors very efficient internal conversion between both triplet states. These two angles are regularly reached by the trajectories (Figure d), providing ample opportunity for ISC and subsequent internal conversion to T1 to take place.

After the transition to T1, the system equilibrates between two distinct local T1 minima, termed “pyramidalized” and “boat” in ref (34). In Figure a, most trajectories exhibit values of the C5=C6 bond length of around 1.38 Å, which is typical for the pyramidalized T1 minimum as well as the S1 minimum. A small fraction of the ensemble shows values around 1.49 Å, corresponding to a population of the boat T1 minimum. Pyramidalization angles p7213 (Figure b) around ±40° indicate the pyramidalized T1 minimum and values close to zero indicate a population of the boat T1 minimum. Detailed examination of the individual trajectories showed that the average population ratio in the two T1 minima is approximately 2:1, in favor of the pyramidalized minimum.

Internal conversion to the ground state is minimal, in line with experimental evidence of high ISC yields and phosphorescence.[30,60] In agreement with this, none of the reported S1/S0 conical intersections[32−34] is observed.

Finally, Figure e presents the temporal evolution of the permanent dipole moment of the active state. After excitation to S2, its value is about 6.4 D, whereas after relaxation to the lower states it quickly changes to 3.8 D and then does not change notably. This implies that the involved excited states (S1, T1, and T2) have very similar dipole moments, which are also very close to the value of the ground state (4.2 D, see also ref (34)). This fact has important consequences for the photochemical behavior of 2TU in polar solvents, where more polar states tend to be stabilized while unpolar states increase in energy. In the present case, a polar solvent would decrease the energy of S2 (πSπ6*), further enhancing S2 → S1 internal conversion, and would leave the remaining states, and hence the remaining dynamics, mostly unchanged. This could offer an explanation for the good agreement of our simulations and time-resolved experiments in solution,[30] even though the simulations neglect solvation effects. This hypothesis is also supported by the very similar absorption maxima reported for different solvents in the literature (see Table S1 in the Supporting Information). In this regard, 2TU seems to be a rare exception from the observation that solvation strongly affects the dynamics of nucleobases.[22,23,26] Especially in uracil, the reported ISC yields are strongly solvent-dependent and vary from a few percent in water to about 50% in aprotic, nonpolar solvents.[61,62] The SHARC method, using CASSCF potentials, yields for uracil in the gas phase an ISC yield of 67% (extrapolated from the data reported in ref (63)).

In conclusion, nonadiabatic dynamics simulations show that the mechanism of ISC is very efficient in the noncanonical nucleobase 2-thiouracil. The global sketch shown in Figure summarizes the deactivation mechanism. After excitation to the first bright state S2 of 1πSπ* character, the system quickly moves to a mostly planar conical intersection with the S1 (1nSπ*) state and decays with a time constant of 59 fs to S1. Internal conversion from S2 to S1 mainly involves bond length inversion followed by pyramidalization. Almost no ground-state relaxation was observed; instead, extremely efficient ISC from S1 to T2 and T1 occurs, with an average time constant of 400 fs, in agreement with time-resolved experiments in solution.[21,30] The most important ISC channel is the indirect S1 → T2 → T1 one, with S1 → T1 and S2 → T2 playing smaller roles. ISC involves strong pyramidalization of the thiocarbonyl group and is promoted by the wagging mode of the same group, which efficiently modulates the gap between 1nSπ* (S1) and 3πSπ* states. Our extrapolated triplet yield is 90 ± 10%, which agrees nicely with the near-unity yields of several experimental studies.[28−30]

Figure 4

Summary of the excited-state dynamics of 2-thiouracil, showing excitation (), internal conversion in the singlet states (59 fs time constant), and intersystem crossing (400 fs).

These results not only are relevant for 2TU but also evidence how DNA and RNA modifications strongly influence their response to ultraviolet light. For example, thionation can abate the photoprotective abilities inherent to the canonical nucleobases, unless relaxation from the reactive triplet states to the ground state takes place sufficiently fast. However, the ability to populate efficiently triplet states can be exploited to generate singlet oxygen with therapeutic aims, as has been experimentally shown in thiouracils.[21,29]