彭 沁，方業(yè)廣，張騰爍，崔剛龍，方維海

（北京師范大學化學學院,理論及計算光化學教育部重點實驗室,北京100875）

天然堿基腺嘌呤、鳥嘌呤、胞嘧啶、胸腺嘧啶和尿嘧啶是脫氧核糖核酸（DNA）和核糖核酸（RNA）的基本組成單元部分，對核酸序列識別、遺傳信息存儲、DNA聚合酶復制等均有重要作用.它們具有許多優(yōu)異的物理化學性質［1］，如光穩(wěn)定性，可防止紫外線照射下的生物光損傷.為了揭示其光保護的分子機制，在過去十幾年中，大量的實驗和理論研究探索了在不同環(huán)境，如真空、溶液、DNA或RNA中，堿基及堿基對的激發(fā)態(tài)性質和動力學行為［2，3］.現在，普遍認為這些堿基的光穩(wěn)定性歸因于它們可以通過錐形交叉結構，從激發(fā)態(tài)超快地衰減到基態(tài)，最終將能量耗散到環(huán)境中［4～10］.

與天然堿基不同，將堿基上的氧原子進行硫取代生成的硫代堿基具有明顯不同的激發(fā)態(tài)和光物理性質［11，12］.在過去幾年中，實驗和理論研究探索了硫代堿基的激發(fā)態(tài)性質和衰減動力學等［13～23］.Harada等［24］發(fā)現4-硫代胸腺嘧啶三重態(tài)量子產率接近100%；隨后的瞬態(tài)吸收實驗表明，4-硫代胸腺嘧啶在水溶液中存在皮秒時間尺度的系間竄躍過程［25］.Crespo-Hernández等［26，27］發(fā)現4-硫代胸腺嘧啶在離子溶液中存在超快的系間竄躍動力學，持續(xù)時間小于240 fs，解釋了較弱的室溫磷光現象.他們也研究了在不同溶液中2-硫代胸腺嘧啶和2-硫代尿嘧啶亞皮秒的系間竄躍動力學，認為S1暗態(tài)在激發(fā)態(tài)弛豫過程中起著重要作用［28］.鄭旭明課題組［29］利用共振拉曼技術探索了2-硫代尿嘧啶在S2態(tài)上的馳豫動力學.最近，Crespo-Hernández等［30～33］進一步研究了2，4-二硫代胸腺嘧啶，發(fā)現其接近100%的三重態(tài)量子產率和超快的系間竄躍動力學.不難發(fā)現，天然堿基是以超快的內轉換過程為主，但是硫代堿基則是以系間竄躍過程為主要的激發(fā)態(tài)弛豫通道.為了理解其奇特的實驗現象背后的物理本質，一些理論課題組利用高精度的電子結構和非絕熱動力學方法系統(tǒng)地研究了它們的激發(fā)態(tài)和光物理性質［34～47］.

由于氧、硫和硒屬于同一族元素，硒代堿基在過去幾年也引起了一些實驗方面的關注［48～55］，但主要集中在電子的基態(tài).Crespo-Hernández等【56】研究了6-硒代鳥嘌呤（6SeG）在溶液中的激發(fā)態(tài)性質，發(fā)現了單重態(tài)到三重態(tài)有超快的系間竄躍和超高的三重態(tài)量子產率；但是，相比于硫代體系，三重態(tài)的壽命較短.此外，Huang課題組［51，57］成功合成了含硒代堿基的核酸，發(fā)現在DNA中4-硒代胸腺嘧啶和6-硒代鳥嘌呤的吸收光譜發(fā)生明顯的紅移.此外，他們在RNA中也發(fā)現2-硒代尿嘧啶的吸收光譜會發(fā)生明顯的紅移［55］.

在計算方面，Russo等［58，59］用含時密度泛函理論（TD-DFT）方法探索了硒代胸腺嘧啶和鳥嘌呤的光物理性質.Manae和Hazra［60］使用多態(tài)完全活化空間二級微擾理論（MS-CASPT2）方法計算了硫和硒代胸腺嘧啶的吸收光譜.Mai等［61］使用MS-CASPT2和二階代數圖表構造法ADC（2），探索了2-硒代尿嘧啶（2SeU）的激發(fā)態(tài)弛豫路徑，并開展了非絕熱動力學模擬.他們認為與2-硫代尿嘧啶相比，2SeU的三重態(tài)壽命減少了3個數量級，主要原因是缺少第二個三重態(tài)的參與、自旋-軌道耦合的增強以及T1/S0勢能面交叉結構和T1極小能量結構的較小能量差.最近，本課題組［62］使用高精度MS-CASPT2方法，研究了6SeG的激發(fā)態(tài)弛豫路徑，發(fā)現增強的T1/S0旋-軌耦合導致了較短的三重態(tài)壽命.隨后繼續(xù)使用MS-CASPT2方法，分別研究了1個和2個硒原子取代的2SeU，4SeU和24SeU的激發(fā)態(tài)性質和光物理機理，發(fā)現不同位點的硒取代會顯著影響垂直和絕熱激發(fā)能、激發(fā)態(tài)性質和弛豫路徑［63］.

DNA環(huán)境非常復雜，如氫鍵作用可能對堿基的光物理性質產生較大的影響，但至今大多數關于硒代堿基的理論計算研究并沒有較為合理地考慮DNA環(huán)境的影響.為了考慮DNA環(huán)境中互補堿基對的氫鍵對堿基的光物理機理的影響，本文使用量子力學/分子力學組合方法和高精度的電子結構方法，即QM（MS-CASPT2//CASSCF）/MM方法，研究和比較了DNA環(huán)境中2-硒代胸腺嘧啶和腺嘌呤堿基對（2SeT-A）以及4-硒代胸腺嘧啶和腺嘌呤堿基對（4SeT-A）的激發(fā)態(tài)性質和光物理過程.

1 計算方法

1.1 計算模型

根據實驗研究提供的結構信息，利用AMBER16程序包中的NAB模塊建立了DNA序列為5′-CTTCTTGTCCG-3′∶3′-GAAGAACAGGC-5′的雙鏈結構［64］，并放置于80 nm×80 nm×80 nm的立方體水盒子中.然后，加入20個鈉離子，使得體系變?yōu)殡娭行?在建立的初始體系基礎上，在前1000步中，整個DNA序列保持凍結，優(yōu)化水分子及鈉離子；在剩下的2000步里對整個體系能量進行最小化處理，不施加任何限制.然后加熱20 ps，使體系從0 K升溫至300 K；最后，進行1 ns的恒溫和恒壓的分子動力學模擬.在分子力場的結構優(yōu)化和動力學模擬中，DNA結構和鈉離子使用Amber ff99sb力場描述［65］，而水分子使用TIP3P模型描述［66］.隨后，在分子動力學模擬得到的穩(wěn)定結構基礎上，將其中一個單鏈的第五個胸腺嘧啶堿基的2和4號位置的氧原子分別替換為硒原子生成含有2-硒和4-硒代胸腺嘧啶堿基（2SeT或4SeT）的DNA雙鏈結構，其中與硒取代胸腺嘧啶堿基配對的是互補鏈的腺嘌呤堿基A.

1.2 QM/MM計算

所有電子結構計算均在量子力學/分子力學（QM/MM）方法框架下進行［67～70］.QM區(qū)域包含2SeT或4SeT及與之配對A堿基（2SeT-A或4SeT-A）.MM區(qū)域則由剩余的DNA結構及所有的離子和水分子構成.對于QM/MM極小能量結構和交叉點結構的優(yōu)化，QM區(qū)域的電子結構方法選用完全活化空間自洽場（CASSCF）方法［71］.在CASSCF計算中，單重態(tài)和三重態(tài)均采用5態(tài)平均計算，每個電子態(tài)的權重為0.2.在單重態(tài)計算中平均了S0～S4這5個電子態(tài)，而在三重態(tài)中平均了T1～T5這5個電子態(tài).由于該方法不能充分考慮動態(tài)電子的相關效應，因此，在QM/MM單點能量計算的時候，QM區(qū)域選用多態(tài)完全活化空間自洽場二級微擾理論（MS-CASPT2）方法［72］.同樣地，在MS-CASPT2計算中均采用5態(tài)平均，態(tài)平均計算參數與CASSCF計算一致.同時，采用0.2原子單位的虛能級移動來避免入侵組態(tài)問題［73］.根據本課題組之前的經驗和測試及近期Zobel等［74］的測試結果，電離和親核勢（IPEA）校正設置為零［75］.同時使用Cholesky分解技術快速而準確地計算雙電子積分［76］，截斷閾值設置為10-4.在CASSCF和MS-CASPT2計算中，活化空間選用14電子11軌道（本文支持信息中圖S1和圖S2）.C，H，O和N原子選用cc-pVDZ基組，重原子Se則選用較大的aug-cc-pVDZ基組［77，78］.旋-軌耦合常數利用平均場近似方案計算［79～81］，表示為

式中：ψI和ψJ是單重態(tài)和三重態(tài)的電子波函數；Hsox，Hsoy和Hsoz是旋-軌算符的x，y和z分量.所有的QM/MM計算使用與TINKER6.3.2軟件包連接的MOLCAS8.0軟件包.

2 結果與討論

2.1 基態(tài)的穩(wěn)定結構和光譜性質

采用QM（CASSCF）/MM方法首先優(yōu)化了2SeT-A和4SeT-A在S0態(tài)上的極小能量結構（圖1，S0-MIN）.盡管這些結構都具有類似平面的對稱性，但是也存在一些明顯的結構差異.當O原子變?yōu)镾e原子時，鍵長明顯變長，如2SeT-A中的C2―Se7鍵長伸長到0.184 nm.此外，在2SeT-A和4SeT-A的S0-MIN結構中，與C―O和C―Se鍵直接相連的化學鍵的長度也會發(fā)生明顯變化.

Fig.1 QM(CASSCF)/MMoptimized minima of 2SeT-A and 4SeT-A base pair in the S0 states

Franck-Condon區(qū)域的激發(fā)態(tài)性質對于理解2SeT-A和4SeT-A的激發(fā)態(tài)弛豫動力學具有重要作用.在優(yōu)化得到的S0極小能量結構處，QM（MS-CASPT2）/MM計算表明2SeT-A和4SeT-A的S1態(tài)是具有1nπ*性質的暗態(tài)，而S2態(tài)是具有1ππ*性質的明態(tài).但是，它們的1nπ*和1ππ*態(tài)的電子組態(tài)涉及的前線軌道不同.對于2SeT-A，1nπ*態(tài)對應Se7原子的n電子激發(fā)到嘧啶和C2-Se7雙鍵π*軌道.相比之下，4SeT-A的1nπ*態(tài)激發(fā)Se8原子的n電子到嘧啶和C4-Se8雙鍵的π*軌道（圖2）.這些差異導致2SeT-A和4SeT-A的1nπ*態(tài)的垂直激發(fā)能不同，其中2SeT-A為336.7 kJ/mol，而4SeT-A降低到243.1 kJ/mol（表1）.有關1ππ*態(tài)，π*軌道與1nπ*態(tài)中的相同，而π軌道不同.2SeT-A中的π軌道主要由Se7原子的p軌道組成；而4SeT-A中是由Se8原子的p軌道組成（圖2）.同樣地，1ππ*態(tài)的垂直激發(fā)能也略有不同，其中2SeT-A為378.2 kJ/mol，而4SeT-A為347.3 kJ/mol（表1）.在2SeT-A和4SeT-A中，具有1ππ*性質的更高的電子激發(fā)單重態(tài)能量約為460.0 kJ/mol，分別比最低的1ππ*態(tài)高約84.0和113.0 kJ/mol，實驗所用激發(fā)波長難以到達.因此，本文主要研究從最低的1ππ*態(tài)出發(fā)的激發(fā)態(tài)弛豫機理.

Fig.2 Molecular orbitals relevant to the lowest two excited singlet states at the Franck-Condon points of 2SeT-A and 4SeT-A

Table 1 QM(MS-CASPT2)/MM calculated vertical excitation energies(kJ/mol)at the S0 minima of 2SeT-A and 4SeT-A

2.2 激發(fā)態(tài)極小能量結構

除了S0態(tài)的穩(wěn)定結構，也用QM（CASSCF）/MM方法優(yōu)化了S2，S1，T2和T1激發(fā)態(tài)的極小能量結構.對于2SeT-A，分別優(yōu)化出了C―Se鍵沿著嘧啶環(huán)平面向上翹（-U）和向下翹（-D）兩套激發(fā)態(tài)結構（本文支持信息圖S3）.因為這兩套激發(fā)態(tài)極小能量結構對應的激發(fā)態(tài)能量比較接近，因此文中主要討論其中C―Se鍵沿著嘧啶環(huán)平面向上翹的結構（加個后綴-U）.圖3和圖4分別為2SeT-A-U的激發(fā)態(tài)結構示意圖及S2，S1，T2和T1激發(fā)態(tài)極小能量結構相對于S0穩(wěn)定結構鍵長的變化值.可以發(fā)現，這些結構的變化主要位于2SeT和4SeT部分，而A堿基的貢獻幾乎可以忽略不計.這些幾何結構特征和上述討論的激發(fā)態(tài)電子結構特征是一致的.與S0態(tài)的穩(wěn)定結構相比，S2態(tài)的C2―Se7鍵長變化最大，伸長了0.039 nm；其次是C6―N1，N1―C2，N3―C4和C4―C5鍵長，分別變化了約0.002，0.001，-0.001和0.001 nm.在S1，T2和T1態(tài)的極小能量結構中，也是C2―Se7鍵長變化最明顯，分別變化了0.023，0.025和0.031 nm，略小于S2態(tài)的C2―Se7鍵長；但是C6―N1，N1―C2，C2―N3，N3―C4和C4―C5鍵長變化相對較大（圖4）.QM（MS-CASPT2）/MM計算表明，2SeT-A-U和2SeT-A-D的S2，S1，T2和T1的絕熱激發(fā)能分別為313.8和312.5、257.3和268.2、257.7和254.4、242.7和256.5 kJ/mol.2SeT-A-D的激發(fā)態(tài)極小能量結構及對應鍵長的變化見圖S4和圖S5（本文支持信息）.

Fig.3 QM(CASSCF)/MMoptimized minima of 2SeT-A-U base pair in the S1,S2,T1 and T2 states

Fig.4 Bond-length changes of QM(CASSCF)/MM optimized excited-state minimum-energy structures of 2SeT-A-U relative to corresponding ground-state structures

對于4SeT-A，除了U和D兩套激發(fā)態(tài)極小能量結構外，還優(yōu)化得到了C―Se鍵與嘧啶環(huán)基本共平面的結構（-P）.同樣地，這3套結構對應的能量均比較接近（本文支持信息圖S13），因此主要討論其中C―Se鍵沿著嘧啶環(huán)平面向上翹的結構（加個后綴-U）.圖5和圖6分別為4SeT-A-U的激發(fā)態(tài)平衡結構及S2，S1，T2和T1相關激發(fā)態(tài)的平衡結構相對于S0平衡結構對應鍵長的變化值.同樣可以發(fā)現，這些激發(fā)態(tài)結構的變化主要發(fā)生于4SeT，而A堿基部分的結構仍然變化較小.與2SeT-A-U的激發(fā)態(tài)極小能量結構不同，4SeT-A-U相對于S0穩(wěn)定結構變化最大的是C4―Se8鍵.其中，S2態(tài)的極小能量結構的C4―Se8鍵長變化最大，增加了0.042 nm；其次為C4―C5，C5―C6和C6―N1鍵長，分別變化了-0.002，0.001和-0.001 nm.在S1，T2和T1態(tài)的結構中，C4―Se8鍵長也是變化最大的，分別為0.015，0.021和0.018 nm，但是小于S2態(tài)的C4―Se8鍵長變化；而N1―C2，C2―N3，C3―N4和C6―N1鍵長相對變化也較大（圖6）.對于4SeT-A-U（4SeT-A-D，4SeT-A-P），QM（MS-CASPT2）/MM計算的S2，S1，T2和T1態(tài)的絕熱激發(fā)能分別為299.6（289.1，297.1），210.5（211.3，212.1），204.2（216.7，206.3）和199.6（208.8，209.6）kJ/mol.對應的4SeT-A-D和4SeT-A-P的激發(fā)態(tài)結構及相關鍵長變化見圖S6～圖S9（本文支持信息）.

Fig.5 QM(CASSCF)/MMoptimized minima of 4SeT-A-U base pair in the S1,S2,T1 and T2 states

Fig.6 Bond-length changes of QM(CASSCF)/MM optimized excited-state minimum-energy structures of 4SeT-A-U relative to corresponding ground-state structures

2.3 交叉點結構

眾所周知，不同勢能面之間的交叉結構對于激發(fā)態(tài)弛豫過程有著重要作用，是連接不同勢能面之間的“橋梁”.對于2SeT-A，在S2-MIN-U和S2-MIN-D的極小能量結構處，S2和S1的能量已經比較接近，QM（MS-CASPT2）/MM計算相對能量分別為313.8/296.2和312.5/297.9 kJ/mol.因此，在這兩個極小能量結構附近，處于S2態(tài)的體系可以有效地衰減到S1態(tài).此外，還優(yōu)化得到了T1/S0的交叉結構T1S0（圖7）.在該結構處，C2―Se7鍵長被拉長到0.246 nm，Se7-C2-N1-N3二面角變?yōu)?25.1°.在QM（MS-CASPT2）/MM計算級別，T1和S0的能量為265.3和252.3 kJ/mol.對于4SeT-A，找到了S2/S1和T1/S0的兩個交叉結構S2S1和T1S0（圖7）.它們的C4―Se8鍵分別拉長到0.278和0.247 nm，而Se8-C6-C2-N1二面角分別為175.8°和146.9°.在S2S1結構，QM（MS-CASPT2）/MM計算的S2和S1態(tài)的能量分別為374.5和359.8 kJ/mol，而在T1S0結構中，計算的T1和S0能量分別為228.4和214.6 kJ/mol.

Fig.7 QM(CASSCF)/MMoptimized intersection structures of 2SeT-A and 4SeT-A

2.4 激發(fā)態(tài)的失活機理

根據優(yōu)化得到的極小能量結構和交叉結構，并結合計算得到的線性內插內坐標（LIIC）路徑，分析得到了2SeT-A和4SeT-A在DNA環(huán)境中的主要激發(fā)態(tài)失活機理.對于2SeT-A體系，當光照到達Franck-Condon區(qū)域的S2態(tài)（FC S2）后，存在2條不同的弛豫路徑，分別到達S2勢能面中的兩個極小能量結構S2-MIN-U和S2-MIN-D.對于路徑Ⅰ，如圖8所示，體系首先弛豫到S2-MIN-U.在該結構處，由于S2和S1態(tài)的能量比較接近，能量差僅為17.6 kJ/mol；因此，可以通過內轉換有效地弛豫到S1態(tài)，然后弛豫到S1-MIN-U，該過程沒有能壘.在S1態(tài)上面，由于S1，T2和T1態(tài)具有較小的能差，因此S1態(tài)到T2和T1態(tài)的系間竄躍過程是允許發(fā)生的.但是，較大的S1/T1旋-軌耦合（395.0 cm-1）使得S1→T1的系間竄躍過程比較有利.到達T1-MIN-U后，2SeT-A可以通過T1/S0交叉點衰減到S0態(tài)，該過程要克服22.6 kJ/mol的能壘.在該結構處，較大的T1/S0旋-軌耦合（432.6 cm-1）將有助于T1→S0的系間竄躍過程，但是較大的能壘使T1態(tài)有可測量的壽命.如圖S10（本文支持信息）所示，路徑Ⅱ與路徑Ⅰ類似.體系首先無能壘弛豫到S2-MIN-D，然后繼續(xù)衰減到S1-MIN-D，最后到T1-MIN-D.不同之處在于，在路徑Ⅱ中，從T1-MIN-D到T1/S0交叉點需要克服73.6 kJ/mol的能壘，使T1態(tài)的壽命更長.

Fig.8 Main excited relaxation path of 2SeT-A in DNA calculated by QM(MS-CASPT2)/MM(the conformation of C─Se bond up along the pyrimidine ring plane)

對于4SeT-A，光激發(fā)到S2態(tài)Franck-Condon區(qū)域后，可以有3個不同的路徑.對于路徑Ⅰ，分子首先弛豫到S2-MIN-U，但是從S2-MIN-U到S2/S1圓錐交叉點需克服60.7 kJ/mol的能量.在S1態(tài)上面，存在大量S1，T2和T1態(tài)近簡并的區(qū)域，如圖9所示.沿著這些路徑，具有較大的S1/T1旋-軌耦合（410.0 cm-1）的S1→T1系間竄躍過程更為有利；而S1→T2系間竄躍過程，由于較小的S1/T2旋-軌耦合（39.1 cm-1）則相對較慢.這些現象符合EI-Sayed規(guī)則［82］.因為S1和T2態(tài)具有相同的nπ*電子結構，而S1和T1態(tài)則具有不同的nπ*和ππ*電子結構.因此，在S1→T1無輻射躍遷過程中，總的角動量可以守恒，有利于系間竄躍發(fā)生.到達T1-MIN-U后，考慮到T1/S0交叉結構處的較大旋-軌耦合（373.0 cm-1），4SeT-A可以進一步衰減到S0態(tài)，但是該過程要克服28.9 kJ/mol的能壘.因此，體系會在T1態(tài)上停留一段時間.路徑Ⅱ和Ⅲ（本文支持信息圖S11和圖S12）與路徑Ⅰ的激發(fā)態(tài)弛豫過程接近，因此，3條路徑都有可能對4SeT-A在DNA環(huán)境中的激發(fā)態(tài)失活做出貢獻，具體的重要性要進一步的非絕熱動力學模擬才能完全確定.

Fig.9 Main excited relaxation path of 2SeT-A in DNA calculated by QM(MS-CASPT2)/MM(the conformation of C—Se bond up along the pyrimidine ring plane)

3 結 論

使用QM（MS-CASPT2//CASSCF）/MM方法，系統(tǒng)研究了在DNA環(huán)境中2-硒和4-硒代胸腺嘧啶和腺嘌呤堿基對2SeT-A和4SeT-A的激發(fā)態(tài)性質和光物理過程.盡管2SeT-A和4SeT-A具有相似的S2（ππ*），S1（nπ*），T2（nπ*）和T1（ππ*）電子激發(fā)態(tài)，但是由于涉及不同的n和π軌道，它們的激發(fā)態(tài)行為存在較大差異，導致2SeT-A的垂直和絕熱激發(fā)能都比4SeT-A高.此外，無論2SeT-A還是4SeT-A在DNA環(huán)境中均存在多個具有不同構象的激發(fā)態(tài)極小能量結構和激發(fā)態(tài)失活途徑.盡管這些激發(fā)態(tài)弛豫路徑有相似之處，基本是按照S2→S1→T1的光物理過程失活到最低的三重T1激發(fā)態(tài)，使T1態(tài)具有可測量的壽命.但是具體的S2→S1內轉換過程在2SeT-A和4SeT-A中還是有顯著區(qū)別.前者基本是個沒有能壘的過程，后者則需要克服60.7 kJ/mol的能壘.總體上可以看出，2SeT-A和4SeT-A在DNA環(huán)境中的激發(fā)態(tài)弛豫路徑主要是系間竄躍到最低三重態(tài).計算工作也表明不同位置的硒取代胸腺嘧啶在DNA環(huán)境中具有顯著不同的激發(fā)態(tài)性質和光物理機理，這些差異將有助于理解硒取代堿基激發(fā)態(tài)性質和光物理機理的復雜環(huán)境效應.

支持信息見http：//www.cjcu.jlu.edu.cn/CN/10.7503/cjcu20210117.