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

（重慶郵電大學生物信息學院，重慶 400065）

0 引言

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

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

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

1 模型與方法

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

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

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

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

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

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

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

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

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

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

2 結果與討論

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

在室溫、不引入激光條件下，經過200 fs的動力學模擬，得到胸腺嘧啶的基態(tài)平衡構型。所得胸腺嘧啶分子為平面構型，其鍵長、鍵角、二面角等結構參數(shù)與文獻報道基本一致，其分子結構和原子編號如圖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方法優(yōu)化得到的基態(tài)平衡狀態(tài)下HOMO與LUMO的能差相吻合。我們對10個不同初始構型進行模擬脈沖輻射通量密度從10 J·m－2到500 J·m－2，得到200余條軌跡，此處僅報道輻射通量密度為315.57 J·m－2的一條軌跡，其余軌跡的動力學性質基本類似。

圖2為胸腺嘧啶激發(fā)態(tài)失活過程的快照圖。初始時刻胸腺嘧啶分子為平面結構（見圖2a），激光作用以后，胸腺嘧啶分子被激發(fā)，其C5－C6的π鍵被破壞，C6原子上連接的H6原子開始畸變（見圖2b）;H6原子畸變的同時C5原子上所連接的甲基也開始畸變（見圖2c，2d）;H6原子畸變程度最大時幾乎與分子平面垂直，此時C6原子凸起呈錐狀（見圖2e）;此后分子恢復平面構型（見圖2f），說明已經衰減至基態(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的躍遷，分子恢復到基態(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ā)生避免相交，此后軌道能量恢復至其初始值。在避免相交時刻HOMO與LUMO能量兼并，二者之間強烈的耦合作用使得電子發(fā)生無輻射躍遷，導致分子返回基態(tài)。

圖4描述了幾個主要結構參數(shù)隨時間的變化。從圖4a可以看出，C5－C6鍵長初始值為1.39?，是典型的共軛體系C=C雙鍵鍵長。激光作用以后，由于胸腺嘧啶分子被激發(fā)，π鍵被破壞，鍵長延長至平均1.50?，變?yōu)棣益I;450 fs以后，由于發(fā)生電子無輻射躍遷，分子返回基態(tài)，導致C5－C6也恢復至雙鍵鍵長。圖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)面在避免相交時刻的強烈扭曲為非絕熱衰減提供了驅動力，導致了無輻射躍遷的發(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)能量轉換的特征。初始時刻分子處于靜止狀態(tài)，勢能為－610 eV，動能為0 eV，激光作用以后，由于電子被激發(fā)，電子能迅速上升，勢能上升至－602 eV，同時有各原子具有了一定的速度，動能略有增加。450 fs時由于電子衰減至基態(tài)，電子能下降為約－605 eV，多余的電子能轉變?yōu)閯幽?，使得動能增加至約3 eV。這也對圖4描述的結構變化做出了解釋，失活后由于分子具有動能，因此處于振動態(tài)，各個鍵長或二面角不是一個常數(shù)而是在一定范圍內波動。在生物體內，分子動能可以傳遞給溶劑（體液），從而使分子弛豫至振動低能級，確保分子不受損傷。

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

3 結論

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

［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］袁帥，王丹，白明澤，等.偶氮苯調控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.