The Behavior of Triplet Thymine in a Model B-DNA Strand. Energetics and Spin Density Localization Revealed by ab initio Molecular Dynamics Simulations

Among the naturally occurring nucleobases, thymine presents the lowest triplet state, hence it represents a hotspot for energy transfer and photosensitization. In turn, the population of the triplet state may lead to thymine dimerization and hence to the production of toxic DNA lesions and has been the subject of intensive theoretical and experimental investigations. Relying on QM/MM molecular dynamics simulations, we have sought to situate the energy of the lowest triplet state of thymine embedded in a B-DNA environment. The energy gap varies between 305 and 329 kJ mol − 1 when a single thymine is treated at the quantum chemistry level, depending on its position in the model double-stranded 16-bp oligonucleotide. The energy of triplet state decreases up to 300 kJ mol − 1 , due to polarization effects, when we consider coupled stacked nucleobases up to the inclusion of four nucleobases. Our value lies in good agreement with the energy inferred experimentally by Miranda and coworkers (270 kJ mol − 1 ), and our theoretical exploration opens the door to investigations toward other more complex and bio-logically relevant environments, such as thymines embedded in nucleosome core particles. Our investigations also provide a reference for further studies using semi-empirical approaches such as density functional-based tight-binding, allowing to further rationalize sequence effects.


INTRODUCTION
Thymine plays a central role in DNA photochemistry. Indeed, direct absorption of UV-light leads to thymine dimerization, producing either cyclobutane pyrimidine dimers (CPD) or 6-4 photoproduct (64-PP) as the most common lesions (1)(2)(3). The former damages are also toxic when they accumulate, as they are either repair-resistant or highly mutagenic and are correlated to the development of malignant skin cancers as the result of unprotected exposure to sunlight (4)(5)(6)(7). Furthermore, thymine is also an hot spot in DNA photosensitization (8)(9)(10). Indeed, its triplet state, the lowest among the natural nucleobases, can be populated via triplet-triplet energy transfer from the excited sensitizer. This process can also be facilitated by the additional stabilization of thymine triplet 3 T due to π-stacking, and the related multichromophoric coupling in the DNA helical structure. Thymine photochemistry has been explored as early as 1960s (11), but the presence of 3 T (12) has triggered renewed attention, notably for delineating the mechanism of triplet-triplet energy transfer with external photosensitizers (13), or with modified nucleobases and oxidative lesions recognized as Trojan Horses (9,14). It has also been shown that the population of 3 T opens rather efficient photochemical channels for dimerization (15,16), even if there is no experimental evidences for the formation of 64-PP mediated by 3 T. The energy of 3 T in a DNA strand has been inferred experimentally by Bosca et al. by assessing the energy transfer efficiency of different photosensitizers, such as norfloxacin and its N4'-acetyl derivative, whose energies are precisely known (17).
The assumption being that only those sensitizers whose triplet state is higher than the one of 3 T embedded in the DNA environment will be able to induce energy transfer. This approach has provided a reference value of 270 kJ mol −1 , which confirms the stabilization of 3 T in DNA and provides a macroscopic average in a complex environment. Yet this value encompasses a manifold of different nucleobases whose direct surroundings differs largely. Indeed, the average value masks many crucial aspects that may modify the local energy alignment, and hence may influence the transfer process. These include the nature of the adjacent bases and sequence effects, the helical embedding and its structural variability, the differential screening induced by the DNA backbone and the presence of additional charges due to salt concentration or other charged species. An important and recent experimental study along 64 possible sequences has delineated some factors tuning direct and triplet sensitized CPD formation in DNA : the effect of flanking bases can increase up to 12-fold CPD formation (18). Also, in biological media, the situation may be complicated even further by the presence of proteins interacting with DNA, which may induce structural and electronic constraints. As a matter of fact, sequence dependence is known to impact energy transfer in DNA oligonucleotides (19). The presence of adjacent thymines may be particularly relevant in stabilizing 3 T due to the formation of excitonic or excimeric intermediates presenting a delocalization of the spin density over multiple bases (20). In turn, the presence of delocalized states, may alter the efficiency of triplet-triplet energy transfer, and hence the photosensitization efficiency, but also the rate of 3 T → 1 S 0 decay, influencing the subsequent DNA photochemistry. A rather large panel of experimental data are available to assess the decay channels of 3 T leading either to the population of the singlet ground state or to the formation of photodamages, for isolated nucleotides or in single-stranded DNA thymines (21,22). In parallel, several computational studies have considered 3 T embedded in B-DNA, to assess characteristic energies and lifetimes for triplet-triplet energy transfer, either with photosensitizers (23) or along thymine stacks (24,25). In 2011, the characteristic time for triplet migration was estimated at 6.4 ns by semi-empirical quantum chemical methods coupled to classical molecular-dynamics simulations (26). The photochemistry of isolated thymine has also been intensively studied by molecular modeling using either static or dynamics approaches. Notably, it has allowed to rationalize the quantum yield of the population of the lowest triplet state (27) or the evolution of Watson-Crick pairing (28). The evolution of the long-lived 3ππ* state has been shown to have a characteristic time which may encompass the nanosecond range, as it has also been confirmed by the study of the decay in DNA single strand (22,29), for thymine and their analogous (30,31). Furthermore, the lifetime has also shown a dependency on the possible formation of excimeric states characterized by the delocalization of the triplet spin density over several stacked nucleobases. However, the precise amount of such delocalization remains to be quantified. In this contribution, we rely on quantum mechanics/molecular mechanics molecular dynamics (QM/MM-MD) simulations to assess the energy and nature of 3 T on a representative 16-bp double-stranded oligonucleotide harboring a 4 thymine TTTT motif in its center. The 16-bp sequence was chosen as it features two alternating 6bp dG-dC blocks, conferring rigidity to the structure and being representative of experimental choices. In addition, four nucleobases are usually considered as the upper scale limit for spin density delocalization (32). By tuning our QM/MM partition to include either only one (1T) or the four consecutive thymines (4T), we may systematically assess delocalization and polarization effects on the energy of the triplet state. Of note, the choice of running QM/MM-MD simulations also provides us with the possibility to assess the influence of DNA flexibility on the dispersion of the triplet energy levels, hence bridging the gap with the macroscopic and averaged experimental determinations.

MATERIALS AND METHODS
Preliminary classic MD simulations. Classical all-atom molecular dynamics (MD) simulations have been performed using the Amber18 package (33) while the set-up has been performed with the Ambertools suite. Classical MD is used to provide a representative starting structure for the DNA duplex prior to perform the QM/MM-MD. The sequence of the chosen 16-bp duplex is represented in Fig. 1. The ff99bsc1 force field was employed for the nucleic acid (34), since it notably includes the so-called Barcelona's correction for the DNA backbone (35,36). The system was solvated in a box of 7858 TIP3P water (37) molecules. Potassium (K + ) counterions were added to neutralize the total charge of the box. Long-range electrostatic interactions were computed using Particle Mesh Ewald (PME) algorithm (38,39). After 1 ns of equilibration, a production run of 100 ns was performed to generate a representative starting structure, which has been chosen based on a cluster analysis of the MD trajectory.
QM/MM molecular dynamics simulations. QM/MM-MD simulations were performed using the Orca (40,41) suite of programs, interfaced with the Amber18 package. The QM partition was selected to encompass either one of the four thymine nucleobases (T7, T8, T9 or T10) or the whole four stacked thymines in order to assess the delocalization character of the triplet state (see Fig. 1). The QM subsystem consists in the nucleobases, with a link atom placed on the dangling N-glycosidic bond. QM/MM-MD trajectories were run for a total of 1 ns with a timestep of 1.0 fs. The resulting QM/MM and QM-only energies have Figure 1. Sequence of the 16-bp double-stranded oligonucleotide. The QM part encompasses either one nucleobase here T8 (a) "minimal 1T system" or the T7T8T9T10 stack (b) "4T system" in the middle of the oligonucleotide. (c) The simulation box for the solvated 16-bp duplex is shown, DNA is represented in cartoon, while K + cations are in yellow van der Waals spheres. The consecutive TTTT strand is also evidenced in ball and stick representation as shown in the zoom inlay. been averaged over the last 500 ps. Note that 5 replicas have been run to increase the statistical sampling. Our chosen global time length is meaningful to provide a good analysis of the behavior of the triplet state in DNA, which can reach the ns scale. The singlet and triplet states have been obtained at density functional theory (DFT) and its time-dependent (TD-DFT) extension level using the M06-2X (42) functional and the Pople's double-ζ 6-31G(d, p) basis set. This level of theory was benchmarked against four other different functionals, and also against the def2-DVZP basis set. This choice may, indeed, provide both an accurate representation of the electronic structure and a sufficiently fast convergence allowing for a reasonable sampling of the phase space. The centroid of the spin density C S has been obtained as the norm of the center of spin vector which is calculated by integrating over all the volume dV the product of the position vector R with the spin density ρ S : From a practical standpoint, C S has been obtained by a numerical integration of the spin density cube files for the selected snapshots.
MD and QM/MM-MD trajectories have been analyzed and visualized using the VMD code (43).

RESULTS
First, we compare the overall averaged energies (ΔE) between 1 S 0 and 3 T 1 along the central TTTT stack obtained from five replicas for each of the thymine in the 1T model (see Table 1). These energies range between 288 and 346 kJ mol −1 along the 20 trajectories, and are subject to important fluctuations for the same QM part. Comparison of the singlet-to-triplet energies (ΔE) suggests that the formation of the triplet state may be more favorable on T9 and T10 than on T8 and T7. Indeed, localization of spin density on the latter bases would lead to an energy penalty of 18 kJ mol −1 .
However, the difference is not systematic along the trajectories, and falls in the range of the accuracy of our energy estimates. Conversely, for a given QM region, the variations are significant and may reach up to 43 kJ mol −1 . Hence, one can hypothesize a scenario in which photosensitization, and triplet energy transfer, is possible due to dynamic and vibrational fluctuations leading to a temporary favorable alignment of the energy levels. Note also that a consistent, although moderate, stabilization of the triplet energy compared with the one of isolated thymine is observed already at this level of theory.
Despite having access to the entire QM/MM-MD trajectories, the identification of the structural factors breaking the symmetry of the system and leading to the preference for the localization of the triplet on T9 and T10, are not self-evident. Indeed, one has to take into account that the average values mask rather important deviations and oscillations and hence, should not be overinterpreted. Furthermore, the high variability of the energy value is also due to a subtle coupling between the DNA vibrational degrees of freedom, the solvent reorganization, and the counterions movements. As a matter of fact, ions have been recognized as important for triggering DNA chemistry (44,45) and photochemistry (46) and in many instance they play a noninnocent role. However, the detailed analysis of the ensemble of these factors, and their cross-talks, goes beyond the scope of the present contribution, and would probably require longer sampling time and sophisticated, or machine-learning-based, analytical tools. Furthermore, as shown in the Supporting Information, the analysis of DNA structural parameters, performed using Curves+ code (47), does not pinpoint any clear differences between the behavior of the four thymines. Indeed, the Watson-Crick network is mostly conserved, in agreement with previous studies that have reported that the 1 S 0 → 3 T 1 excitation of thymine does not affect the base pairing. Unsurprisingly, the overall bending angle is not changed upon excitation, varying from 6 to 18°(see Supporting Information).
Interestingly, when including the entire TTTT strand (4T model) in the QM partition, we observe a nonnegligible stabilization of 3 T whose ΔE drops to an average of 299 kJ.mol -1 . This effect is indeed compatible with the enhanced treatment of multichromophoric coupling and polarization effects, which are not included in our force field based on fixed-charge and hence in the MM embedding. It is also in line with experimental evidences pointing to the decrease of the triplet energy in the DNA environment (17). However, if the decrease of the triplet energy is undeniable, a finer analysis of the QM/MM-MD simulation reveals a more contrasted behavior. Indeed, the 4T model experiences the most prominent oscillations, globally spanning an interval of 99.6 kJ mol −1 on the ensemble of the five replicas. As a matter of fact, while trajectory number 1 provides a relative low energy (268 kJ mol −1 ), trajectory 5 yields an average ΔE of 316 kJ mol −1 , that represents a value that is comparable to the ones obtained for the 1T model. Globally our results are in reasonable agreement with the experimental determination of the triplet energy in thymine, the value of 270 kJ mol −1 being comprised in our confidence interval. However, our data also offer a more complex scenario, in which the energy of the triplet state is strongly dependent on the precise, and instantaneous, arrangement of the DNA and solvent structural parameters, leading to the coexistence of higher and lower energy situations, which thus can be accessible for triplet energy transfer or not. This aspect can have a strong influence in dictating the overall characteristic time for the energy transfer, but also the final quantum yield, especially in presence of competitive channels, such as those due to triplet quenchers. While we have clearly underlined the role of DNA, and most specific the effects of π-stacked chromophores in modulating the energy alignment between singlet and triplet states, the question of a localized or delocalized state, as well as of spin density mobility, remains unanswered.
To tackle this issue, we report the spin density calculated over equidistant snapshots extracted from the QM/MM-MD simulation. An illustration of some representative snapshots is also provided in Fig. 2, while the corresponding movie can be seen in the Supporting Information. It appears evident that the spin density is, for the majority of the time, highly localized over only one thymine nucleobase. Also, coherently with the results of the 1T model, a marked preference for T9 is observed. Interestingly enough, while partial delocalization between two consecutive chromophores was observed only very rarely, it is clear that πstacking is mostly acting by enhancing polarization rather than by inducing long-range spread excitonic states. If the spin density is localized, it is nonetheless extremely mobile; indeed, during the time spanned by our QM/MM-MD, we may observe the population of virtually all the thymine nucleobases. Furthermore, the triplet transfer seems to proceed via a noncoherent hop-based mechanism rather than by coherent energy transfer, which should require persistent delocalized states.
To assess more quantitatively the localization of the spin density, we report in Fig. 3 the evolution of the position of its centroid C S projected on the DNA vertical axis. We can first point out a coherent behavior between the different replicas, again evidencing a highly localized spin-density. Indeed, in most of the cases, C S is localized in close proximity to the center of mass of T9, the small deviation presented may also be seen as an evidence of spin polarization. All the five replicas also agree in evidencing a high preference for spin localization on the central thymines rather than at the periphery. This result, which may appear in disagreement with the 1T model can however be easily rationalized considering the border effect of the QM/MM partition, which disfavor the peripheral chromophores. The strand orientation favors the thymines T9 and T10 over T8 and T7, which indeed correspond to the experimental evidences for the thymines at 3'-end being more prone to stabilize their triplet states (18).
However, and coherently with what observed from the representative movie provided as ESI, we may still evidence a rather high spin mobility for all the replicas and especially for trajectories 4 and 5. Indeed, while we can see that the spin density may localize on virtually all the thymines, we also evidence, especially for trajectory 4, the establishment, at around 800 ps, of an oscillatory behavior with the unpaired electrons rapidly swapping between T9 and T8. Interestingly, in this case the increased distance of C S with the center of mass of the corresponding base could also suggest a stronger amount of delocalization, or at least an increased polarization.
Our results may delineate a coherent picture of the behavior of the triplet state of thymine embedded in B-DNA and of the role of the environment. While π-stacking appears fundamental in stabilizing the triplet state, and hence in opening the way to photosensitization, the spin density appears localized, and spin polarization, rather than excitonic delocalization, should be responsible for the altered photophysical properties.
Of note, we should also underline that the reduced role of delocalization, as evidenced in our work, could also result from the fact that using MD simulations, we may account for thermal disorder effects. On the contrary, idealized and ordered model systems, having a perfect symmetry, are, obviously, nonphysically over emphasizing delocalization.

DISCUSSION AND CONCLUSIONS
In this contribution and thanks to the use of extended QM/MM-MD simulations, we have succeeded in analyzing the energy alignment and the spin localization of the lowest triplet state of thymine embedded in a B-DNA environment. In particular, while we have shown that the triplet state is usually highly localized, with the participation of only one nucleobase, we have also clearly shown that the spin polarization induced by the thymine environment, and correctly captured only by the larger QM/MM partition, is fundamental to provide a better description of the Interestingly, we have also shown that the triplet energy is undergoing rather important variations, especially in the case of the 4T model, as evidenced by the dispersion of the corresponding ΔE values. This observation opens important questions regarding the photosensitization process. Indeed, the coexistence of low and high energy triplet states may also point to a more complicated photosensitization mechanisms, which could strongly affect both the characteristic times and the quantum yields. Experimental determinations, which asses only the average macroscopic energy difference present clear limitations, since they stay blind to subtle local effects susceptible to play a key role in the further photophysical outcome. Furthermore, the variability of the energy alignment, together with the observed high mobility of the spin density, also necessitate to reconsider sequence effects and the possible presence of hot spots that may lead to triplet accumulation in close spatial proximity, and hence eventually to DNA cluster lesions.
These aspects may be seen even more crucial and dramatic in more complex environment as compared with our solvated, rather short, DNA double strand. In particular, the mapping of the energy of triplet states in nucleosomal environment is particularly intriguing and could allow to underline the presence of specific sinks or spin traps.
Even though we have provided the first atomistic and electronic-resolved determination of the triplet energy of thymine triplet state, some important aspects are still missing. The use of semi-empirical methodologies, such as density functional tight binding (DFTB) may allow to significantly increase the statistical sampling, hence providing more correlation from which the role of precise structural and dynamic factors could be more easily inferred. Furthermore, DFTB may allow a more systematic exploration of sequence effects, which have been shown to play an important role (18), but also to increase the size of the QM partition to avoid spurious side effects. Finally, efficient semiempirical methods will be necessary to bring the study to the nucleosomal realm (48) without an insurmountable computational burden. However, DFTB, and more in general semi-empirical methods, requires a tailored parameterization against either experimental results or ab initio methods. Due to the difficulty of obtaining nucleobase-resolved experimental determinations, we believe that our work may also serve as a solid base for the benchmark of the less expensive methods that we preview to use in the near future. Such methods will offer a sound cross-talk with experimental trends which have been delineated recently, without the limitation of considering a unique double-stranded model oligomer.
Acknowledgements-The authors thank the SYSPROD project and AXELERA Pôle de Compétitivité for financial support (PSMN Data Center). AM also thanks ANR (Agence Nationale de la Recherche) and CGI (Commissariatà l'Investissement d'Avenir) for their financial support of this work through Labex SEAM (Science and Engineering for Advanced Materials and devices) ANR 11 LABX 086, ANR 11 IDEX 05 02. This article is part of a special issue to celebrate the scientific achievements of Jean Cadet. The authors are indebted to the scientific legacy brought by him and to the considerable advancements he has performed allowing a significant improvement of the understanding of DNA lesions from a chemical and biological point of view.

SUPPORTING INFORMATION
Additional supporting information may be found online in the Supporting Information section at the end of the article:  Figure S1. Distribution of the global bending calculated for each QM/MM-MD replica of the 1T model in the triplet state, including T7, T8, T9 and T10. Figure S2. Whiskers plots for the Shear, Stretch and Stagger parameters obtained for the ensemble of the QM/MM-MD of the 5 replicas in the triplet state. Note that the parameters have been calculated for the nucleobase in the QM partition. Figure S3. Whiskers plots for the Buckle, Propel and Opening parameters obtained for the ensemble of the QM/MM-MD of the 5 replicas in the triplet state. Note that the parameters have been calculated for the nucleobase in the QM partition. Table S1. Position of the Spin density Centroid on a representative snapshot obtained with different density functionals. Energies of transition S0 ->T (in kJ.mol -1 ) were also recomputed.