胸腺嘧啶紫外损伤与自修复机理的计算机模拟

舒坤贤，马 静，袁 帅，豆育升

（重庆邮电大学生物信息学院，重庆 400065）

0 引言

由于DNA碱基腺嘌呤（adenine，A）、鸟嘌呤（guanine，G）、胞嘧啶（cytosine，C）和胸腺嘧啶（thymine，T）能够有效地吸收紫外光，因此DNA分子在紫外光的照射下易发生电子激发。激发态的DNA分子通过光化反应产生一系列光化产物，从而影响DNA的复制和转录功能，引起遗传突变，导致皮肤癌等疾病［1－3］。然而实验测量到的光化产物的量子产率很低，约为1%［4］，表明绝大多数激发了的DNA在形成损伤前已经通过各种途径衰减至基态。DNA对光辐射去活化的这一超快的生物学过程，对维持生命体的稳定性意义重大，也引发了人们研究碱基激发态性质的强烈兴趣。

随着超快激光实验和量子化学计算方法的发展，碱基激发态无辐射衰减机制的研究取得重大进展［5－8］。实验显示，在水溶液中，受激发的单碱基荧光量子产率约为10－4［9］，因此激发了的碱基主要是通过超快无辐射跃迁失活而回到电子基态。Kohler研究组［10］在2000年第一次精确测量出 DNA碱基1ππ*激发态的寿命为亚皮秒级，并提出碱基的超快无辐射跃迁由电子激发态和基态的圆锥相交（conical intersection，CI）引起。量子化学计算发现所有其他的天然碱基及其衍生物都有着相似的从电子激发态到基态的衰减通道［6］。嘧啶碱基的C5－C6键的扭曲和嘌呤碱基的相关双键扭曲振动是碱基形成圆锥相交的主要构型变化［11－12］。本研究组的动力学模拟研究［13］发现，腺嘌呤的C2和氨基平面外振动是引起分子的电子能量转化为热能量的主要原因。

理解DNA分子中空间结构对激发态衰减的影响时一个极具挑战性的课题，对认识DNA光稳定性和光损伤的机理及其重要。尽管该领域已经有大量的实验和理论研究，但是采用分子动力学方法的研究很少。理解一个化学反应最终需要弄清从反应物到产物的过程中每个原子是如何运动的，然而受限于过于简单的反应坐标，分子整体的运动往往难以顾及。本文中我们采用半经典动力学模型模拟胸腺嘧啶光损伤与自修复反应，在我们的研究中可以包括每一个反应自由度以及能量随时间的变化。模拟结果可以提供碱基体系的超短寿命激发态的形成和失活的详细动力学信息，为皮肤癌等疾病的预防提供一定的理论支持。

1 模型与方法

目前，国际上用于研究光化学反应的动力学过程的半经典模型主要有美国的Martinez研究组［14］提出的AIMS（ab initiomultiple spawning）模型和英国的Robb研究组［15］提出的QM/MM模型。AIMS模型采用从头算的方法处理电子和核运动，能得到较精确的结果，但是计算量很大。QM/MM模型对体系中反应活性部分用量子化学计算，对相对非活性部分用分子力学方法处理，大大降低了计算量。但是因为分子力学不能处理电子激发态结构，所以QM/MM模型在涉及电子激发的光化学反应中受到了限制。

本文所采用的模拟激光诱导光化学反应的计算模型，称为半经典电子－辐射－离子动力学（semiclassically electronic radiation ion dynamics，SERID） 模型［16］。半经典是指电子运动由量子力学方法确定，而用经典力学方法处理辐射场和核运动。这种模拟方法考虑激光场对电子结构的作用，并采用非绝热动力学的近似方法（Ehrenfest方法）和紧束缚近似来建立模型。

在该方法中，电子运动由含时薛定谔（Schrödinger）方程描述为

激光脉冲辐射场的矢势A可以通过含时的Peierls替代耦合到电子哈密顿中:

（3）式中:Urep是离子 -离子排斥势;∂H/∂Xla，∂S/∂Xla，和∂Urep/∂Xla是H，S和U相对于每一个原子位置的梯度。所有的梯度矩阵元可经过数值求解法确定。

SERID模型具有以下优点:

1）在电子 Hamiltonian中通过时间相关的Peierls替代，引入了辐射场的矢量A（x，t），即在模型中引入激光与分子的相互作用。这是SERID模型与其他类似模型相比的独特之处，计算中引入此种相互作用使我们能够研究激光脉冲对光化学反应结果的影响。

2）电子结构采用半经验的DFTB（density functional based tight binding）方法计算，基函数取定域在原子上的非正交基。仅处理价电子，内层电子和核一起构成惰性的原子（离子）实，原子实之间的排斥作用采用参量化处理，因此能量以及力的计算有很高的效率。

3）采用直接动力学计算方案，即不需要预先构建势能面，能量和力的确定采用固定时间间隔的即用计算（on the fly）模式。

这些特点保证了SERID可以进行足够大体系的动力学模拟。我们使用该模型成功地模拟了几个紫外辐射下的生物分子光化学反应，包括偶氮苯异构化反应［17－19］，腺嘌呤无辐射超快失活机理［13］，胸腺嘧啶光致二聚反应［20］和胸腺嘧啶二聚体的光致解离［21］反应。这些模拟结果都与实验及量子化学计算结果符合。

2 结果与讨论

2.1 优化基态构型

在室温、不引入激光条件下，经过200 fs的动力学模拟，得到胸腺嘧啶的基态平衡构型。所得胸腺嘧啶分子为平面构型，其键长、键角、二面角等结构参数与文献报道基本一致，其分子结构和原子编号如图1所示。将平衡构型再进行1 000 fs的模拟，间隔100 fs，取得10个构型作为初始构型，初始构型包括了原子坐标和速度、加速度等动力学信息。

图1 胸腺嘧啶分子结构与原子编号Fig.1 Conformation and atomic labeling for the thyminemolecule

2.2 模拟结果与讨论

选择频率为5.0 eV（相当于248 nm波长），最大半波宽（full width of halfmaxium，fwhm）为50 fs（飞秒）的激光脉冲照射T分子，脉冲频率与SERID方法优化得到的基态平衡状态下HOMO与LUMO的能差相吻合。我们对10个不同初始构型进行模拟脉冲辐射通量密度从10 J·m－2到500 J·m－2，得到200余条轨迹，此处仅报道辐射通量密度为315.57 J·m－2的一条轨迹，其余轨迹的动力学性质基本类似。

图2为胸腺嘧啶激发态失活过程的快照图。初始时刻胸腺嘧啶分子为平面结构（见图2a），激光作用以后，胸腺嘧啶分子被激发，其C5－C6的π键被破坏，C6原子上连接的H6原子开始畸变（见图2b）;H6原子畸变的同时C5原子上所连接的甲基也开始畸变（见图2c，2d）;H6原子畸变程度最大时几乎与分子平面垂直，此时C6原子凸起呈锥状（见图2e）;此后分子恢复平面构型（见图2f），说明已经衰减至基态。

图3描述体系的电子布居数和轨道能量随时间的变化。图3a显示了电子跃迁的细节，激光作用以后，约 1.6个电子从 HOMO激发至 LUMO，表明248 nm紫外光的辐射可以使胸腺嘧啶分子从基态变为电子激发态;在450 fs时发生电子从LUMO至HOMO的跃迁，分子恢复到基态，激发态寿命为350 fs，与激发态寿命是亚皮秒级的文献报道［10］一致。从图 3b看出，初始时刻 LUMO能量为－2.0 eV，激光作用以后，LUMO 能量显著下降，同时HOMO能量略上升。450 fs时HOMO与LUMO能级接近，能隙仅为0.02 eV，发生避免相交，此后轨道能量恢复至其初始值。在避免相交时刻HOMO与LUMO能量兼并，二者之间强烈的耦合作用使得电子发生无辐射跃迁，导致分子返回基态。

图4描述了几个主要结构参数随时间的变化。从图4a可以看出，C5－C6键长初始值为1.39Å，是典型的共轭体系C=C双键键长。激光作用以后，由于胸腺嘧啶分子被激发，π键被破坏，键长延长至平均1.50Å，变为σ键;450 fs以后，由于发生电子无辐射跃迁，分子返回基态，导致C5－C6也恢复至双键键长。图4b和图4c中 C4－C5－C6－H6和C4－C5－C6－N1二面角的变化表征了嘧啶环扭曲的情况。约200 fs以后，C4－C5－C6－H6二面角从开始的180°迅速增大，450 fs时达到最大值310°;C4－C5－C6－N1二面角在450 fs时增大到82°。分子环面在避免相交时刻的强烈扭曲为非绝热衰减提供了驱动力，导致了无辐射跃迁的发生。

图4 C5－C6键长、C4－C5－C6－H6二面角及C4－C5－C6－N1二面角随时间的变化图;Fig 4 Variationswith time of C5－C6 bond、C4－C5－C6－H6 dihedral angle and C4－C5－C6－N1 dihedral angle

图5描述了胸腺嘧啶分子激发－失活过程动能和势能（此处势能包括电子能与核－核排斥势）随时间的变化，从中可以发现能量转换的特征。初始时刻分子处于静止状态，势能为－610 eV，动能为0 eV，激光作用以后，由于电子被激发，电子能迅速上升，势能上升至－602 eV，同时有各原子具有了一定的速度，动能略有增加。450 fs时由于电子衰减至基态，电子能下降为约－605 eV，多余的电子能转变为动能，使得动能增加至约3 eV。这也对图4描述的结构变化做出了解释，失活后由于分子具有动能，因此处于振动态，各个键长或二面角不是一个常数而是在一定范围内波动。在生物体内，分子动能可以传递给溶剂（体液），从而使分子弛豫至振动低能级，确保分子不受损伤。

图5 体系动能和势能随时间的变化Fig 5 Variations with time of kinetic energy and potential energy of thymine.

3 结论

采用半经典动力学方法对激光诱导下胸腺嘧啶分子最低激发态的失活过程进行了模拟。模拟脉冲频率为5.0 eV，fwhm=50 fs，相当于紫外辐射。模拟发现胸腺嘧啶在紫外辐射下会生成电子激发态，激发态的胸腺嘧啶分子发生环扭曲，破坏其原有的平面结构，严重影响碱基的生命功能，且存在产生自由基而引起突变的可能。但是由于激发态寿命极短，仅为350 fs，激发态的胸腺嘧啶分子通过H6原子和C6原子的畸变发生无辐射失活而衰减至基态，失活后释放的电子能转化为分子动能，并通过振动弛豫使其消散于环境，恢复了胸腺嘧啶分子的平面结构，从而最大程度的避免了紫外损伤。

［1］TAYLOR JS.Unraveling the Molecular Pathway from Sunlight to Skin Cancer［J］.Acc Chem Res，1994，27（3）:76－82.

［2］VINK A A，ROZA L.Biological consequences of cyclobutane pyrimidine dimers［J］.JPhotochem Photobiol B，2001，65（2）:101－104.

［3］ BEUKERSR，EKER A，LOHMAN P.50 years thymine dimer［J］.DNA Repair，2008，7（3）:530－543.

［4］ ROCA－SANJUÁN D， OLASO－GONZáLEZ G，GONZáLEZ－RAMíREZ I，et al.Molecular Basis of DNA Photodimerization:Intrinsic Production of Cyclobutane Cytosine Dimers［J］.JAm Chem Soc，2008，130（32）:10768－10779.

［5］ CRESPO－HERNáNDEZC，COHEN B，HARE P，et al.Ultrafast Excited－State Dynamics in Nucleic Acids ［J］.Chem Rev，2004，104（4）:1977－2020.

［6］SHUKLA M，LESZCZYNSKI J.Review Electronic spectra，excited state structures and interactions of nucleic acid bases and base assemblies:a review［J］.JBiomol Struct Dyn，2007，25（1）:93－118.

［7］SAIGUSA H.Excited－state dynamics of isolated nucleic acid bases and their clusters［J］.JPhotochem Photobiol C，2006，7（4）:197－210.

［8］de VRIESM，HOBZA P.Gas－Phase Spectroscopy of Biomolecular Building Blocks［J］.Annu Rev Phys Chem，2007，58:585－612.

［9］ DANIELSM，HAUSWIRTH W.Fluorescence of the Purine and Pyrimidine Bases of the Nucleic Acids in Neutral Aqueous Solution at 300°K ［J］.Science，1971，171:675－677.

［10］ PECOURT J－M L，PEON J，KOHLER B.Ultrafast Internal Conversion of Electronically Excited RNA and DNA Nucleosides in Water［J］.JAm Chem Soc，2000，122（38）:9348－9349.

［11］ ZGIERSKIM Z，PATCHKOVSKIIS，FUJIWARA T，et al.The role of out－of－plane deformations in subpicosecond internal conversion of photoexcited purine bases:Absence of the ultrafast decay channel in propanodeoxyguanosine［J］.Chem Phys Lett，2007，440（1－3）:145－149.

［12］ SERRANO－ANDRÉSL，MERCHÁN M，BORIN A C.A Three－State Model for the Photophysics of Adenine ［J］.Chem Euro J，2006，12（25）:6559－6571.

［13］ LEIY B，YUAN S，DOU Y S，et al.Detailed Dynamics of the Nonradiative Deactivation of Adenine:A Semiclassical Dynamics Study［J］.J Phys Chem A，2008，112（37）:8497－8504.

［14］ BEN－NUN M，MARTíNEZ T J.Ab initio quantum molecular dynamics［J］.Adv Chem Phys，2002，121:439－512.

［15］ BERNARDI F，OLIVUCCIM，ROBB M A.Simulation of MC－SCF results on covalent organic multi－bond reactions:molecularmechanicswith valence bond（MM－VB）［J］.JAm Chem Soc，1992，114（5）:1606－1616.

［16］ DOU Y S，TORRALVA B，ALLEN R J.Semiclassical electron－radiation－ion dynamics（SERID）and cis－trans photoisomerization of butadiene［J］.Mod Optics，2003，50（15－17）:2615－2643.

［17］袁帅，王丹，白明泽，等.偶氮苯调控DNA解链的机理研究［J］.重庆邮电大学学报:自然科学版，2009，21（6）:821－824.

YUAN S，WANG D，BAIM Z，et.al.Mechanism research ofmelting of DNA duplex by photoisomerization of azobenzene［J］.Journal of Chongqing University of Posts and telecommunications（Natural Science Edition），2009，21（6）:821－824.

［18］ YUAN S，WU W，WEN Z，et al.A semiclassical dynamics study of the photoisomerization ofmethyl－substituted azobenzene［J］.Mol Phys，2010，108:3431－3441.

［19］ YUAN S，DOU Y S，WUW F，et al.Why Does trans－Azobenzene Have a Smaller Isomerization Yield forππ*Excitation Than for nπ* Excitation? ［J］.JPhys Chem A，2008，112（51）:13326－13334.

［20］ ZHANGW Y，YUAN S，LI A Y，et al.Photoinduced Thymine Dimerization Studied by Semiclassical Dynamics Simulation［J］.JPhys Chem C，2010，114（12）:5594－5601.

［21］ DOU Y S，XIONG S S，WU W F，et al.Photoinduced dissociation of cyclobutane thymine dimer studied by semiclassical dynamics simulation［J］.J Photoch Photobio B，2010，101（1）:31－36.