胸腺嘧啶紫外損傷與自修復(fù)機理的計算機模擬

舒坤賢，馬 靜，袁 帥，豆育升

（重慶郵電大學(xué)生物信息學(xué)院，重慶 400065）

0 引言

由于DNA堿基腺嘌呤（adenine，A）、鳥嘌呤（guanine，G）、胞嘧啶（cytosine，C）和胸腺嘧啶（thymine，T）能夠有效地吸收紫外光，因此DNA分子在紫外光的照射下易發(fā)生電子激發(fā)。激發(fā)態(tài)的DNA分子通過光化反應(yīng)產(chǎn)生一系列光化產(chǎn)物，從而影響DNA的復(fù)制和轉(zhuǎn)錄功能，引起遺傳突變，導(dǎo)致皮膚癌等疾病［1－3］。然而實驗測量到的光化產(chǎn)物的量子產(chǎn)率很低，約為1%［4］，表明絕大多數(shù)激發(fā)了的DNA在形成損傷前已經(jīng)通過各種途徑衰減至基態(tài)。DNA對光輻射去活化的這一超快的生物學(xué)過程，對維持生命體的穩(wěn)定性意義重大，也引發(fā)了人們研究堿基激發(fā)態(tài)性質(zhì)的強烈興趣。

隨著超快激光實驗和量子化學(xué)計算方法的發(fā)展，堿基激發(fā)態(tài)無輻射衰減機制的研究取得重大進展［5－8］。實驗顯示，在水溶液中，受激發(fā)的單堿基熒光量子產(chǎn)率約為10－4［9］，因此激發(fā)了的堿基主要是通過超快無輻射躍遷失活而回到電子基態(tài)。Kohler研究組［10］在2000年第一次精確測量出 DNA堿基1ππ*激發(fā)態(tài)的壽命為亞皮秒級，并提出堿基的超快無輻射躍遷由電子激發(fā)態(tài)和基態(tài)的圓錐相交（conical intersection，CI）引起。量子化學(xué)計算發(fā)現(xiàn)所有其他的天然堿基及其衍生物都有著相似的從電子激發(fā)態(tài)到基態(tài)的衰減通道［6］。嘧啶堿基的C5－C6鍵的扭曲和嘌呤堿基的相關(guān)雙鍵扭曲振動是堿基形成圓錐相交的主要構(gòu)型變化［11－12］。本研究組的動力學(xué)模擬研究［13］發(fā)現(xiàn)，腺嘌呤的C2和氨基平面外振動是引起分子的電子能量轉(zhuǎn)化為熱能量的主要原因。

理解DNA分子中空間結(jié)構(gòu)對激發(fā)態(tài)衰減的影響時一個極具挑戰(zhàn)性的課題，對認識DNA光穩(wěn)定性和光損傷的機理及其重要。盡管該領(lǐng)域已經(jīng)有大量的實驗和理論研究，但是采用分子動力學(xué)方法的研究很少。理解一個化學(xué)反應(yīng)最終需要弄清從反應(yīng)物到產(chǎn)物的過程中每個原子是如何運動的，然而受限于過于簡單的反應(yīng)坐標，分子整體的運動往往難以顧及。本文中我們采用半經(jīng)典動力學(xué)模型模擬胸腺嘧啶光損傷與自修復(fù)反應(yīng)，在我們的研究中可以包括每一個反應(yīng)自由度以及能量隨時間的變化。模擬結(jié)果可以提供堿基體系的超短壽命激發(fā)態(tài)的形成和失活的詳細動力學(xué)信息，為皮膚癌等疾病的預(yù)防提供一定的理論支持。

1 模型與方法

目前，國際上用于研究光化學(xué)反應(yīng)的動力學(xué)過程的半經(jīng)典模型主要有美國的Martinez研究組［14］提出的AIMS（ab initiomultiple spawning）模型和英國的Robb研究組［15］提出的QM/MM模型。AIMS模型采用從頭算的方法處理電子和核運動，能得到較精確的結(jié)果，但是計算量很大。QM/MM模型對體系中反應(yīng)活性部分用量子化學(xué)計算，對相對非活性部分用分子力學(xué)方法處理，大大降低了計算量。但是因為分子力學(xué)不能處理電子激發(fā)態(tài)結(jié)構(gòu)，所以QM/MM模型在涉及電子激發(fā)的光化學(xué)反應(yīng)中受到了限制。

本文所采用的模擬激光誘導(dǎo)光化學(xué)反應(yīng)的計算模型，稱為半經(jīng)典電子－輻射－離子動力學(xué)（semiclassically electronic radiation ion dynamics，SERID） 模型［16］。半經(jīng)典是指電子運動由量子力學(xué)方法確定，而用經(jīng)典力學(xué)方法處理輻射場和核運動。這種模擬方法考慮激光場對電子結(jié)構(gòu)的作用，并采用非絕熱動力學(xué)的近似方法（Ehrenfest方法）和緊束縛近似來建立模型。

在該方法中，電子運動由含時薛定諤（Schr?dinger）方程描述為

激光脈沖輻射場的矢勢A可以通過含時的Peierls替代耦合到電子哈密頓中:

（3）式中:Urep是離子 -離子排斥勢;?H/?Xla，?S/?Xla，和?Urep/?Xla是H，S和U相對于每一個原子位置的梯度。所有的梯度矩陣元可經(jīng)過數(shù)值求解法確定。

SERID模型具有以下優(yōu)點:

1）在電子 Hamiltonian中通過時間相關(guān)的Peierls替代，引入了輻射場的矢量A（x，t），即在模型中引入激光與分子的相互作用。這是SERID模型與其他類似模型相比的獨特之處，計算中引入此種相互作用使我們能夠研究激光脈沖對光化學(xué)反應(yīng)結(jié)果的影響。

2）電子結(jié)構(gòu)采用半經(jīng)驗的DFTB（density functional based tight binding）方法計算，基函數(shù)取定域在原子上的非正交基。僅處理價電子，內(nèi)層電子和核一起構(gòu)成惰性的原子（離子）實，原子實之間的排斥作用采用參量化處理，因此能量以及力的計算有很高的效率。

3）采用直接動力學(xué)計算方案，即不需要預(yù)先構(gòu)建勢能面，能量和力的確定采用固定時間間隔的即用計算（on the fly）模式。

這些特點保證了SERID可以進行足夠大體系的動力學(xué)模擬。我們使用該模型成功地模擬了幾個紫外輻射下的生物分子光化學(xué)反應(yīng)，包括偶氮苯異構(gòu)化反應(yīng)［17－19］，腺嘌呤無輻射超快失活機理［13］，胸腺嘧啶光致二聚反應(yīng)［20］和胸腺嘧啶二聚體的光致解離［21］反應(yīng)。這些模擬結(jié)果都與實驗及量子化學(xué)計算結(jié)果符合。

2 結(jié)果與討論

2.1 優(yōu)化基態(tài)構(gòu)型

在室溫、不引入激光條件下，經(jīng)過200 fs的動力學(xué)模擬，得到胸腺嘧啶的基態(tài)平衡構(gòu)型。所得胸腺嘧啶分子為平面構(gòu)型，其鍵長、鍵角、二面角等結(jié)構(gòu)參數(shù)與文獻報道基本一致，其分子結(jié)構(gòu)和原子編號如圖1所示。將平衡構(gòu)型再進行1 000 fs的模擬，間隔100 fs，取得10個構(gòu)型作為初始構(gòu)型，初始構(gòu)型包括了原子坐標和速度、加速度等動力學(xué)信息。

圖1 胸腺嘧啶分子結(jié)構(gòu)與原子編號Fig.1 Conformation and atomic labeling for the thyminemolecule

2.2 模擬結(jié)果與討論

選擇頻率為5.0 eV（相當(dāng)于248 nm波長），最大半波寬（full width of halfmaxium，fwhm）為50 fs（飛秒）的激光脈沖照射T分子，脈沖頻率與SERID方法優(yōu)化得到的基態(tài)平衡狀態(tài)下HOMO與LUMO的能差相吻合。我們對10個不同初始構(gòu)型進行模擬脈沖輻射通量密度從10 J·m－2到500 J·m－2，得到200余條軌跡，此處僅報道輻射通量密度為315.57 J·m－2的一條軌跡，其余軌跡的動力學(xué)性質(zhì)基本類似。

圖2為胸腺嘧啶激發(fā)態(tài)失活過程的快照圖。初始時刻胸腺嘧啶分子為平面結(jié)構(gòu)（見圖2a），激光作用以后，胸腺嘧啶分子被激發(fā)，其C5－C6的π鍵被破壞，C6原子上連接的H6原子開始畸變（見圖2b）;H6原子畸變的同時C5原子上所連接的甲基也開始畸變（見圖2c，2d）;H6原子畸變程度最大時幾乎與分子平面垂直，此時C6原子凸起呈錐狀（見圖2e）;此后分子恢復(fù)平面構(gòu)型（見圖2f），說明已經(jīng)衰減至基態(tài)。

圖3描述體系的電子布居數(shù)和軌道能量隨時間的變化。圖3a顯示了電子躍遷的細節(jié)，激光作用以后，約 1.6個電子從 HOMO激發(fā)至 LUMO，表明248 nm紫外光的輻射可以使胸腺嘧啶分子從基態(tài)變?yōu)殡娮蛹ぐl(fā)態(tài);在450 fs時發(fā)生電子從LUMO至HOMO的躍遷，分子恢復(fù)到基態(tài)，激發(fā)態(tài)壽命為350 fs，與激發(fā)態(tài)壽命是亞皮秒級的文獻報道［10］一致。從圖 3b看出，初始時刻 LUMO能量為－2.0 eV，激光作用以后，LUMO 能量顯著下降，同時HOMO能量略上升。450 fs時HOMO與LUMO能級接近，能隙僅為0.02 eV，發(fā)生避免相交，此后軌道能量恢復(fù)至其初始值。在避免相交時刻HOMO與LUMO能量兼并，二者之間強烈的耦合作用使得電子發(fā)生無輻射躍遷，導(dǎo)致分子返回基態(tài)。

圖4描述了幾個主要結(jié)構(gòu)參數(shù)隨時間的變化。從圖4a可以看出，C5－C6鍵長初始值為1.39?，是典型的共軛體系C=C雙鍵鍵長。激光作用以后，由于胸腺嘧啶分子被激發(fā)，π鍵被破壞，鍵長延長至平均1.50?，變?yōu)棣益I;450 fs以后，由于發(fā)生電子無輻射躍遷，分子返回基態(tài)，導(dǎo)致C5－C6也恢復(fù)至雙鍵鍵長。圖4b和圖4c中 C4－C5－C6－H6和C4－C5－C6－N1二面角的變化表征了嘧啶環(huán)扭曲的情況。約200 fs以后，C4－C5－C6－H6二面角從開始的180°迅速增大，450 fs時達到最大值310°;C4－C5－C6－N1二面角在450 fs時增大到82°。分子環(huán)面在避免相交時刻的強烈扭曲為非絕熱衰減提供了驅(qū)動力，導(dǎo)致了無輻射躍遷的發(fā)生。

圖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描述了胸腺嘧啶分子激發(fā)－失活過程動能和勢能（此處勢能包括電子能與核－核排斥勢）隨時間的變化，從中可以發(fā)現(xiàn)能量轉(zhuǎn)換的特征。初始時刻分子處于靜止?fàn)顟B(tài)，勢能為－610 eV，動能為0 eV，激光作用以后，由于電子被激發(fā)，電子能迅速上升，勢能上升至－602 eV，同時有各原子具有了一定的速度，動能略有增加。450 fs時由于電子衰減至基態(tài)，電子能下降為約－605 eV，多余的電子能轉(zhuǎn)變?yōu)閯幽?，使得動能增加至約3 eV。這也對圖4描述的結(jié)構(gòu)變化做出了解釋，失活后由于分子具有動能，因此處于振動態(tài)，各個鍵長或二面角不是一個常數(shù)而是在一定范圍內(nèi)波動。在生物體內(nèi)，分子動能可以傳遞給溶劑（體液），從而使分子弛豫至振動低能級，確保分子不受損傷。

圖5 體系動能和勢能隨時間的變化Fig 5 Variations with time of kinetic energy and potential energy of thymine.

3 結(jié)論

采用半經(jīng)典動力學(xué)方法對激光誘導(dǎo)下胸腺嘧啶分子最低激發(fā)態(tài)的失活過程進行了模擬。模擬脈沖頻率為5.0 eV，fwhm=50 fs，相當(dāng)于紫外輻射。模擬發(fā)現(xiàn)胸腺嘧啶在紫外輻射下會生成電子激發(fā)態(tài)，激發(fā)態(tài)的胸腺嘧啶分子發(fā)生環(huán)扭曲，破壞其原有的平面結(jié)構(gòu)，嚴重影響堿基的生命功能，且存在產(chǎn)生自由基而引起突變的可能。但是由于激發(fā)態(tài)壽命極短，僅為350 fs，激發(fā)態(tài)的胸腺嘧啶分子通過H6原子和C6原子的畸變發(fā)生無輻射失活而衰減至基態(tài)，失活后釋放的電子能轉(zhuǎn)化為分子動能，并通過振動弛豫使其消散于環(huán)境，恢復(fù)了胸腺嘧啶分子的平面結(jié)構(gòu)，從而最大程度的避免了紫外損傷。

［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，F(xiàn)UJIWARA 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］袁帥，王丹，白明澤，等.偶氮苯調(diào)控DNA解鏈的機理研究［J］.重慶郵電大學(xué)學(xué)報:自然科學(xué)版，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.