基于概念密度泛函理論磷酸酯類反應(yīng)性物質(zhì)毒性預(yù)測

丁曉琴，丁俊杰，李大禹，潘里，裴承新

1　引言

大多數(shù)有機磷酸酯類化合物都是乙酰膽堿酯酶不可逆抑制劑。不可逆抑制劑不同于一般的配體與受體相互作用，它涉及到舊化學鍵的斷裂和新化學鍵的形成，即發(fā)生化學反應(yīng)的過程。因此在構(gòu)建與受體相互作用的化學反應(yīng)性物質(zhì)的定量構(gòu)效關(guān)系(QSAR)或定量構(gòu)性關(guān)系(QSPR)方程時，會遇到一些獨特的困難與問題1,2。目前，QSAR使用的描述參數(shù)大多數(shù)是非反應(yīng)性的，例如傳統(tǒng)的立體參數(shù)、電性參數(shù)、疏水參數(shù)，以及后來不斷發(fā)展的各種分子結(jié)構(gòu)性質(zhì)描述參數(shù)等3-7，這些參數(shù)中的多數(shù)沒有涉及到化合物在發(fā)生化學反應(yīng)過程中的電子轉(zhuǎn)移變化，從而不能有效和精確地描述化合物在體內(nèi)與酶或受體等發(fā)生化學反應(yīng)的作用機制，進而導(dǎo)致構(gòu)建的反應(yīng)性物質(zhì)的 QSAR或QSPR模型預(yù)測性能較差。

隨著量子化學理論、計算方法以及計算機技術(shù)的迅速發(fā)展，特別是近十幾年來概念密度泛函理論(Conceptual Density Functional Theory，CDFT)的發(fā)展，使得化學反應(yīng)性物質(zhì)的毒性虛擬評估成為可能。概念密度泛函理論是密度泛函理論(DFT)的一個重要分支，又稱為密度泛函活性理論或化學密度泛函理論8-11。它是從 DFT理論中把與化學相關(guān)的概念和原理提取出來，通過理論推導(dǎo)，獲得參數(shù)的計算方法12-14。這些參數(shù)主要包括電負性、硬度、軟度、福井函數(shù)、親電性等一系列新型的量子化學描述符指數(shù)。一方面通過分子的化學勢來描述分子在體系中的反應(yīng)性能；另一方面通過分子中各個原子的電荷密度來描述體系中可能發(fā)生化學反應(yīng)的重點部位15-19。概念密度泛函理論方法得到的一系列新型電子結(jié)構(gòu)描述參數(shù)，是一類新的專門用于描述化學反應(yīng)性物質(zhì)的結(jié)構(gòu)指數(shù)。本文應(yīng)用概念密度泛函理論，對磷酸脂類乙酰膽堿酯酶不可逆抑制劑的反應(yīng)性指數(shù)進行系統(tǒng)的理論計算，探索該類描述指數(shù)在構(gòu)建精確 QSPR模型方面的應(yīng)用，進而實現(xiàn)對磷酸酯類反應(yīng)性物質(zhì)毒性的高精度預(yù)測。這方面的研究內(nèi)容在國內(nèi)外還未見相關(guān)的報道。

2　理論與計算方法

2.1　理論

在概念密度泛函理論中，通過兩種微擾對分子體系基態(tài)能量Ε的影響，來研究各種物理化學性質(zhì)的變化。根據(jù)Hohenberg-Kohn定理，兩種微擾一個是體系總電子數(shù)Ν的變化，另一個是外勢ν的變化，分別以dΝ微擾和dν微擾表示。對于分子體系，外勢 ν通常指原子核形成的庫倫勢，外勢的變化可以是幾何構(gòu)型或構(gòu)象的改變，也可以是與其它分子相互作用時相對位置的變化。在微擾的影響下，體系性質(zhì)的改變可以用敏感度系數(shù)(sensitivity coefficient)來表示。我們感興趣的敏感度系數(shù)是體系基態(tài)對微擾dΝ和微擾dr的一階、二階響應(yīng)函數(shù)的變化，甚至三階和更高價的響應(yīng)函數(shù)(其真正的物理和化學意義正在探索研究和發(fā)展過程中)。通過對分子體系的能量泛函Ε[ρ(r)]=Ε[Ν, r]求微商，可以得到各種敏感度系數(shù)，包括化學勢、電負性、硬度、軟度、福井函數(shù)等等。下面對本文研究密切相關(guān)的幾個反應(yīng)性參數(shù)進行簡要說明：

一階響應(yīng)函數(shù)：可以獲得化學勢 μ和電子密度ρ等。

化學勢：

電子密度：

根據(jù)有限差分原理，化學勢近似從計算分子垂直電離勢I和親和勢A中獲得?；瘜W勢μ的負值是電負性χ。

二階響應(yīng)函數(shù)：可以獲得體系硬度、軟度、福井函數(shù)等。

整體硬度：

體系整體硬度(global hardness, η)表示的是當外勢不變，化學勢隨著電子數(shù)目的改變而變化的速率。

整體軟度：

體系整體軟度(global softness, S)與體系整體硬度之間存在倒數(shù)關(guān)系。

親電指數(shù)：

局域硬度20：

局域硬度又稱硬度密度(hardness density)，描述在外勢不變的條件下，體系化學勢隨 r處電子密度改變而改變的難易程度。

局域軟度：

局域軟度用以衡量在外勢不變的情況下，體系某點的電子密度對體系化學勢微擾的敏感程度。整體軟度與局域軟度的關(guān)系S是：S = ∫s(r)dr，而局域軟度則是：

上式局域軟度包含了福井函數(shù)和整體軟度的信息，即包含了相對反應(yīng)物而言的整體反應(yīng)性-硬度或軟度。

福井函數(shù)：

福井函數(shù)用來描述分子內(nèi)的反應(yīng)順序，即同一個分子結(jié)構(gòu)中不同位置的相對反應(yīng)活性信息，而局域軟度可用來研究分子間的反應(yīng)順序。這兩個都是量子化學反應(yīng)性最重要的定量描述參數(shù)。同時，局域軟度還可度量體系化學勢對某區(qū)域外勢微擾的敏感程度，因此有可能包含分子中不同位置相對反應(yīng)性的信息。

福井函數(shù)是一個局域性質(zhì)，即分子中不同的位置有不同的值。因此，福井函數(shù)為研究化學反應(yīng)的區(qū)域選擇性或確定生物體系的活性中心提供了有力的工具。表示在外勢不變的條件下，分子體系的總電子數(shù)發(fā)生改變時，分子中區(qū)域j點的電子數(shù)變化率；或者分子體系的總電荷發(fā)生改變時，區(qū)域j點的電荷變化率。應(yīng)用有限差分近似，福井函數(shù)f(r)可表示為，

親核福井函數(shù)：表示當分子得到1個電子時，第j原子的電荷變化率。

親電福井函數(shù)：表示當分子失去1個電子時，第j原子的電荷變化率。

自由基攻擊福井函數(shù)：

表示當分子的某一原子受到自由基攻擊時的電荷變化率。

上式中 qj(Ν)，qj(Ν + 1)，qj(Ν - 1)分別表示分子為中性、陰離子和陽離子時，分子中原子j所帶的電荷；fj+(r)和 fj-(r)分別為親核進攻指數(shù)和親電進攻指數(shù)，表示分子中的原子j給電子和得電子能力的強弱，而 fj0(r)為自由基攻擊指數(shù)，其值大小反映自由基進攻的能力。

前線分子軌道是控制化學反應(yīng)的難易程度、反應(yīng)區(qū)域及立體選擇性的重要因素：最高占有軌道(HOMO)控制著親電反應(yīng)位置，而最低空軌道(LUMO)是控制著親核反應(yīng)位置。

2.2　計算方法

本文通過對磷酸酯類化合物的結(jié)構(gòu)性質(zhì)分析，選擇了 15個具有大鼠腹腔半數(shù)致死劑量(LD50)活性的各種結(jié)構(gòu)類型化合物作為訓(xùn)練集21，包括有機磷酸酯類殺蟲劑和部分神經(jīng)性毒劑，毒性范圍跨越 5個數(shù)量級以上。首先分別在氣相和水溶液相(CPCM模型)條件下，采用B3LYP/6-311++G(2d,3p)基組對分子進行了結(jié)構(gòu)優(yōu)化。然后我們采用 MP2/6-311++G(2d,3p)及B3LYP/6-311++G(2d,3p)兩個基組，分別在氣相和水溶液條件下，進行分子結(jié)構(gòu)描述參數(shù)計算。其中，NBO性質(zhì)分別進行了中性分子(分子的電子數(shù)為Ν)，得到一個電子(電子數(shù)為Ν + 1)，失去一個電子(電子數(shù)為Ν - 1)條件下的計算。所有訓(xùn)練集和測試集分子的最低能量構(gòu)象的結(jié)構(gòu)與原子編號參見Supporting Information-圖 S1。

計算得到分子的整體描述參數(shù)有：分子的垂直電離勢I、垂直親和勢A、電子化學勢μ、絕對硬度 η、親電性指數(shù) ω以及分子的前線軌道能量等。分子局域描述指數(shù)有：收斂的原子福井函數(shù)、原子NBO凈電荷以及Wiberg鍵級等。上述計算獲得的整體和局域的分子描述參數(shù)和疏水性參數(shù)等原始計算結(jié)果詳見 Supporting Information-圖S2。

采用上述計算的參數(shù)，應(yīng)用逐步回歸統(tǒng)計分析方法，構(gòu)建四個條件(MP2和B3LYP、氣相和溶液相)的線性方程，確定最合理的計算條件，化合物的毒性采用大鼠腹腔半數(shù)致死劑量的負對數(shù)(pLD50)。在此基礎(chǔ)上，進一步采用遺傳最小偏二乘方法(G/PLS)，構(gòu)建線性/非線性QSPR方程，主成分數(shù)分別設(shè)為 1-6，方程的長度設(shè)為 3-10，優(yōu)化并選擇合理的方程組用于測試集的毒性預(yù)測。

本文分子模型的構(gòu)建采用GaussView 5.09軟件，量子化學計算采用Gaussian09軟件22，在曙光TC3600刀片服務(wù)器上完成。相關(guān)統(tǒng)計分析采用Cerius2軟件23完成。另外應(yīng)用ACD/Lab12.024學術(shù)下載版軟件和HyperChem7.025評估程序計算了分子疏水性logP值、分子體積、表面積等參數(shù)。

3　結(jié)果與討論

3.1　分子反應(yīng)性分析

計算結(jié)果表明，四個計算條件得到的結(jié)果，整體趨勢基本一致。為了便于說明，現(xiàn)將其中的一個條件，即：B3LYP/6-311++G(2d,3p)/gas的計算結(jié)果列于表 1-表 3(僅列出 15個訓(xùn)練集化合物的計算結(jié)果，包括了部分叔胺氮游離態(tài)和質(zhì)子化態(tài)的結(jié)果)。

從上述結(jié)果可見：大多數(shù)分子的親電反應(yīng)中心都在磷原子上。但也有例外情況，如一些計算結(jié)果反映出化學反應(yīng)有可能甚至優(yōu)先發(fā)生在碳原子上，這可能與發(fā)生化學反應(yīng)的類型有關(guān)；反應(yīng)中心磷原子上的正電荷和磷的親電進攻反應(yīng)性并沒有明顯的相關(guān)性；分子中的最低鍵級，并不一定就是與酶反應(yīng)中的離去基團，與得失電子后，分子鍵級的變化過程有關(guān)。例如含有 P―S鍵或P―N鍵的分子，在Ν - 1個電子的狀態(tài)下，其鍵級的降低變化就比較大；任何單個反應(yīng)性指數(shù)與毒性似乎沒有明顯的相關(guān)規(guī)律，這可能與其作用機制和動力學過程中的復(fù)雜性有關(guān)，多種因素的共同作用引起乙酰膽堿酯酶的活性抑制、毒性及老化。

表1　訓(xùn)練集分子的pLD50毒性、前線軌道能量和能差、磷原子的NBO電荷和親電反應(yīng)指數(shù)fp-(r)Table 1　The pLD50 value, frontier orbital energy and gap, NBO charge and fp-(r) in training sets.

表2　訓(xùn)練集分子整體反應(yīng)性指數(shù)Table 2　The global reactivity index for the training sets at B3LYP/6-311++G(2d,3p)/gas level.

3.2　氮原子不同質(zhì)子化狀態(tài)及構(gòu)象對分子反應(yīng)性影響

3.2.1質(zhì)子化狀態(tài)對分子反應(yīng)性影響

Amiton、VX、GV和 EDMM 四個分子(化學結(jié)構(gòu)見Supporting Information-S1)，離去基團側(cè)鏈上的乙胺基叔氮原子質(zhì)子化成為胺基陽離子后，反應(yīng)性指數(shù)計算結(jié)果的共同特點是：分子的HOMO和LUMO軌道能量顯著降低，能量差增加；磷原子上的正電性增加；在氣相和溶液相狀態(tài)下，MP2方法計算氮原子的不同質(zhì)子化狀態(tài)都預(yù)測到磷中心的親電反應(yīng)性，而B3LYP方法在氮原子的非質(zhì)子化狀態(tài)，氣相和溶液相都沒有預(yù)測到磷中心的親電反應(yīng)性，氮原子質(zhì)子化后，磷原子的親電性發(fā)生了明顯的轉(zhuǎn)化，由非親電反應(yīng)中心轉(zhuǎn)到較強的親電反應(yīng)中心；同時，分子的垂直電離勢I增加、親和勢 A降低、I-A的能差顯著增加，絕對硬度η增強、親電性ω有所降低，說明質(zhì)子化后，磷中心的親電進攻反應(yīng)性增強。

3.2.2構(gòu)象對分子反應(yīng)性影響

相同分子和質(zhì)子化狀態(tài)，不同構(gòu)象的反應(yīng)性指數(shù)計算結(jié)果表明，對大多數(shù)分子性質(zhì)，例如前線軌道能量、原子電荷、鍵級、分子整體的反應(yīng)性指數(shù)等，計算結(jié)果有微小的差別，相對于其它的誤差，對回歸結(jié)果不會產(chǎn)生本質(zhì)的影響，采用分子的最低能量構(gòu)象計算分子性質(zhì)即可。但 MP2方法有時對構(gòu)象的反應(yīng)過分敏感，在同一個分子中，具有相同性質(zhì)的碳，構(gòu)象位置上稍有差異，對原子的親電反應(yīng)性就產(chǎn)生不可思議的差別，例如在二異丙基氟磷酸酯(Diisopropyl)分子中，兩個P―O―C鍵上的碳原子親電反應(yīng)性應(yīng)該基本類似，但計算結(jié)果顯示，C6具有僅次于磷的第二個親電反應(yīng)中心的性質(zhì)，另一個C8則基本沒有親電反應(yīng)性，結(jié)果見表4。而B3LYP方法，在氣相和溶液相都沒有預(yù)測到C6和C8的親電反應(yīng)性的較大差異。

另一個方面，在大多數(shù)分子中磷中心的親電反應(yīng)性的計算值，MP2方法明顯要高于B3LYP方法計算結(jié)果，在分子的鍵級變化方面也有類似的問題，說明MP2方法計算結(jié)果的不穩(wěn)定性。

表3　關(guān)鍵原子間NBO鍵級及它們的得失1個電子后變化率Table 3　Pivotal NBO bond order and their variability after obtaining 1e or losing 1e.

3.3　反應(yīng)性指數(shù)的QSPR研究

3.3.1QSPR方程構(gòu)建

根據(jù)分子總體反應(yīng)性指數(shù)分析，以及氮原子離子化狀態(tài)、分子構(gòu)象對反應(yīng)性指數(shù)的影響，我們首先建立氮原子為質(zhì)子化狀態(tài)的 QSPR方程。應(yīng)用逐步回歸分析方法，分別構(gòu)建訓(xùn)練集化合物在四組條件下計算的描述符與毒性pLD50的QSPR方程，得到方程 13至方程 16，其中 r2是回歸方程的決定系數(shù)，F(xiàn)test值是顯著性檢驗，XV_r2是交叉驗證相關(guān)系數(shù)，BS_r2是 Bootstrap方法相關(guān)系數(shù)，Νobs表示構(gòu)建方程時訓(xùn)練集的個數(shù)。

方程13 (B3LYP/gas)：

方程14 (B3LYP/water)：

方程15 (MP2/gas)：

方程16 (MP2/water)：

表4　二異丙基氟磷酸酯優(yōu)化結(jié)構(gòu)及四組條件的計算結(jié)果比較Table 4　Optimizing structure, atomic number of diisopropyl phosphofluoridate and comparison of calculation results with four different calculating level.

從前四個方程可以看出，B3LYP和MP2氣相條件的計算結(jié)果要優(yōu)于溶液相條件的結(jié)果，其原因可能是由于QSPR方程的毒性采用LD50，LD50反應(yīng)的是所發(fā)生的反應(yīng)在體內(nèi)蛋白分子環(huán)境中進行，蛋白環(huán)境下的介電常數(shù)要遠小于水環(huán)境，因此氣相的計算結(jié)果更為可靠。且兩個氣相計算結(jié)果中，B3LYP/gas條件的計算結(jié)果更優(yōu)于MP2得到的結(jié)果，其交叉驗證相關(guān)系數(shù)XV_r2= 0.722，為四組條件的最佳結(jié)果。為了深入地比較分析，我們同時采用逐步回歸方法，對相同計算參數(shù)下的 B3LYP/gas方法，氮非質(zhì)子化形式計算的描述符進行了回歸，得到了方程17。

方程17 (B3LYP/gas，氮非質(zhì)子化)：

方程17的決定系數(shù)和交叉驗證相關(guān)系數(shù)等都有明顯的下降，即采用氮的非質(zhì)子化形式，相關(guān)性遠不如采用氮的質(zhì)子化形式。

同時，我們采用G/PLS統(tǒng)計回歸方法，構(gòu)建QSPR方程，通過變化主成分數(shù)和方程長度，合理地選擇方程，并分別采用線性回歸、樣條函數(shù)回歸、線性和樣條函數(shù)的結(jié)合回歸、非線性回歸等多種方法進行統(tǒng)計分析。結(jié)果表明，直接采用線性回歸，不如采用樣條函數(shù)回歸及線性和樣條函數(shù)結(jié)合的回歸獲得的結(jié)果好；另外，因非線性回歸的結(jié)果較差，故沒有采納。

3.3.2QSPR方程的pLD50預(yù)測

根據(jù)上述比較分析，我們選擇方程13作為預(yù)測方程組的QSPR方程之一。同時，選擇G/PLS方程中，滿足 r2＞ 0.89； XV_r2＞ 0.83； BS_r2＞0.89條件的10個方程組，共11個方程構(gòu)成預(yù)測方程組，對含15個化合物的訓(xùn)練集和含6個化合物的測試集進行 pLD50的預(yù)測，其結(jié)果見表 5中Pred.1CDFT的預(yù)測結(jié)果，(詳細參見 Supporting Informtion-S3)。

為了系統(tǒng)比較分析，我們同時采用Cerius2軟件中的方法22，進行傳統(tǒng)2D-QSPR參數(shù)計算、方程組構(gòu)建與毒性預(yù)測，其結(jié)果見上表 5中 Pred.2 2D-QSPR的預(yù)測結(jié)果，(詳細參見 Supporting Informtion-S3，其逐步回歸方程18在表5注釋中列出)。從逐步回歸的2D-QSPR方程看出，r2= 1，且其 Ftest很大，出現(xiàn)顯著的過擬合現(xiàn)象，顯然不合理。從表 5可以看出，基于概念密度泛函方法的預(yù)測結(jié)果，明顯優(yōu)于采用常規(guī)的2D-QSPR方法，除化合物TEP的預(yù)測誤差稍大外，其余三個化合物都相當不錯；而常規(guī)的2D-QSPR方法有三個化合物的預(yù)測誤差較大，只有一個化合物的預(yù)測結(jié)果較好。測試集中24-optfreq和27-optfreq化合物用這兩種方法的預(yù)測結(jié)果具有顯著地差別，CDFT方法預(yù)測結(jié)果毒性相對較高，而常規(guī)的2D-QSPR預(yù)測結(jié)果的毒性較低。但據(jù)文獻報道26，這兩個化合物是實驗測定酶抑制活性logki最高的兩個化合物，這與基于CDFT方法的預(yù)測結(jié)果基本吻合。

表5　根據(jù)CDFT和傳統(tǒng)2D-QSPR方法，方程組對訓(xùn)練集和測試集化合物pLD50預(yù)測結(jié)果Table 5　pLD50 prediction results for the training set and the test set with the established equations average based on CDFT method and conventional 2D-QSPR.

根據(jù)以上研究，我們初步得出以下結(jié)果：離去基團乙胺基上氮的質(zhì)子化形式是與乙酰膽堿酯酶相互結(jié)合的重要區(qū)域；以質(zhì)子化形式進行該類分子的反應(yīng)性指數(shù)理論計算和 QSPR方程構(gòu)建是合理的。氮的質(zhì)子化不僅增加一個受體負電區(qū)相互作用的結(jié)合點，同時增強了以磷原子為中心親電進攻的反應(yīng)性能；從對構(gòu)象的敏感性和 QSPR方程的計算結(jié)果分析，B3LYP/6-311++G(2d,3p)/氣相條件的結(jié)果更能準確地反應(yīng)結(jié)構(gòu)與毒性的差異，而MP2/6-311++G(2d,3p)的氣相和水溶液相條件得到的結(jié)果，對反應(yīng)性能的差異表現(xiàn)過于敏感，使得計算結(jié)果不穩(wěn)定。

4　結(jié)論

通過對訓(xùn)練集和測試集化合物的反應(yīng)性指數(shù)理論計算，QSPR方程組構(gòu)建，以及外部測試集預(yù)測和交叉驗證預(yù)測，我們發(fā)現(xiàn)應(yīng)用基于概念密度泛函理論的化學反應(yīng)性指數(shù)，進行乙酰膽堿酯酶不可逆抑制劑的毒性預(yù)測是完全可行的。研究結(jié)果表明：①大多數(shù)化合物的親電進攻的反應(yīng)中心發(fā)生在磷原子上，某些化合物的乙胺基氮原子上的質(zhì)子化狀態(tài)，能顯著地增強磷原子的親電性；②總的來看，構(gòu)象對反應(yīng)性描述指數(shù)的計算結(jié)果影響較小，但MP2方法有時對構(gòu)象反應(yīng)過分敏感，相似環(huán)境的原子計算值差別較大；③B3LYP/6-311++G(2d,3p)/gas條件得到的結(jié)果優(yōu)于其它條件；④應(yīng)用基于概念密度泛函理論的反應(yīng)性描述指數(shù)構(gòu)建的 QSPR模型，能夠更準確地表達分子結(jié)構(gòu)與毒性之間關(guān)系。

概念密度泛函方法提出的物理量具有明確的物理意義，能很好描述化學反應(yīng)性的結(jié)構(gòu)特點，因此越來越受到理論化學與(藥物)應(yīng)用化學家的廣泛關(guān)注與重視。由于磷酸酯類反應(yīng)性物質(zhì)的高毒性，使得實驗的危險系數(shù)高，數(shù)據(jù)資料不全或難以獲得。所以采用計算機虛擬評估，應(yīng)用構(gòu)建的高精度 QSPR模型，預(yù)測未知化合物的毒性，是一個很好的選擇。概念密度泛函理論在材料、藥物化學等領(lǐng)域的QSA(T)R應(yīng)用，僅在最近幾年才略有文獻報道，并有持續(xù)增長的趨勢。本文將最新發(fā)展的量子力學概念密度泛函理論方法，應(yīng)用于這類特殊反應(yīng)性化合物的毒理性能評估和環(huán)境毒理評價，為減少動物的使用量，節(jié)約資源，降低環(huán)境污染等做一些努力?？赏麨樘剿髁姿狨ヮ惒豢赡嬉种苿┡c乙酰膽堿酯酶的作用機制開拓新思路，并為評估該類化合物的生態(tài)和環(huán)境毒性提供理論依據(jù)。

Supporting lnformation:available free of charge νia the internet at http：//www.whxb.pku.edu.cn.

(1) Mekenyan, O. G.; Veith, G. D. SAR and QSAR in Εnνiron. Res.1994, 2, 129. doi： 10.1080/10629369408028844

(2) Katagi, K. Reν. Εnνiron. Contam. Toxicol. 2002, 175, 79.doi： 10.1007/978-1-4757-4260-2

(3) Karelson, M.; Lobanov, V. S. Chem. Reν. 1996, 96, 1027.doi： 10.1021/cr950202r

(4) Donald, M. M.; Karen, M. B.; Irwin, K.; Richard, E. S. Arch.Toxicol. 2006, 80, 756. doi： 10.1007/s00204-006-0120-2

(5) Ding, J. J.; Ding, X. Q.; Pan, L.; Chen, J. S. Acta. Phys. -Chim. Sin.2014, 30, 2157. [丁俊杰, 丁曉琴, 李大禹, 潘里, 陳冀勝. 物理化學學報, 2014, 30, 2157.] doi： 10.3866/PKU.WHXB201409171

(6) Ding, J. J.; Ding, X. Q.; Zhao, L. F.; Chen, J. S. Acta Pharm. Sin.2005, 40, 340. [丁俊杰, 丁曉琴, 趙立峰, 陳冀勝. 藥學學報,2005, 40, 340.] doi： 10.3321/j.issn：0513-4870.2005.04.011

(7) Katritzky, A. R.; Kuanar, M.; Slavov, S.; Dennis Hall, C. Chem. Reν.2010, 110, 5714. doi： 10.1021/cr900238d

(8) Parr, R. G.; Yang, W. Annu. Reν. Phys. Chern. 1995, 46, 701.doi： 10.1146/annurev.pc.46.100195.003413

(9) John, C. H. J. Am. Chem. Soc. 2010, 132, 7558.doi： 10.1021/ja1030744

(10) Geerlings, P. K.; De Profit, F.; Langenaeker, W. Chem. Reν. 2003,103, 1793. doi： 10.1021/cr990029p

(11) Liu, S. -B. Acta Phys. -Chim. Sin. 2009, 25, 5. [劉述斌. 物理化學學報, 2009, 25, 5.] doi： 10.3866/PKU.WHXB20090332

(12) Bueno, P. R.; Miranda, D.A. Phys. Chem. Chem. Phys. 2017, 19,6184. doi： 10.1039/c6cp02504h.

(13) James, S. M. A.; Junia, M.; Paul, W. A. J. Chem. Theory Comput.2007, 3, 358. doi： 10.1021/ct600164j

(14) Pérez, P.; Yepes, D.; Jaque, P.; Chamorro, E.; Domingo, L. R.;Rojas, R. S.; Toro-Labbé, A. Phys. Chem. Chem. Phys. 2015, 17,10715. doi： 10.1039/c5cp00306g

(15) Domingo. L. R.; Ríos-Gutiérrez. M.; Pérez P. Molecules 2016, 21,748. doi： 10.3390/molecules21060748

(16) Chattaraj, P. K.; Roy, D. R. Chem. Reν. 2007, 107, PR46.doi： 10.1021/cr078014b

(17) Sablon, N.; Proft, F. D.; Ayers, P. W.; Geerlings, P. K. J. Chem.Theory Comput. 2010, 6, 3671. doi： 10.1021/ct1004577

(18) Tim, F.; Sablon, N.; Proft, F. D.; Ayers, P. W.; Geerlings, P. K.J. Chem. Theory Comput. 2008, 4, 1065. doi： 10.1021/ct800027e

(19) Semenyuk, Y. P.; Morozov, P. G.; Burov, O. N.; Kletskii, M. K.;Lisovin, A. V.; Kurbatov, S. V.; Terrier, F. Tetrahedron 2016, 72,2254. doi： 10.1016/j.tet.2016.03.024

(20) Ayers, P. W.; Parr, R. G. J. Chem. Phys. 2008, 128, 184108.doi： 10.1063/1.2918731

(21) http：//www.drugfuture.com/toxic/search.aspx. (accessed March 28,2013).

(22) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb,M. A.; Cheeseman, J. R.; Scalmani, G.; Barone, V.; Mennucci, B.;Petersson, G. A.; et al. Gaussian 09, Revision B.04; Wallingford CT,Pittsburgh, PA： Gaussian Inc., 2009.

(23) Cerius2, Version 4.5; Accelrys Inc.： San Diego, CA 92121, USA,1999.

(24) ACD lab 12.0 software; Advanced Chemistry Development, Inc.：Canada, 2010.

(25) HyperChem7.0 (Beta1.04 for Evaluation copy) Software; Hypercube,Inc.： Gainesville, 2002.

(26) Victor, E. K.; Eugene, N. M.; Anatoly, G. A. QSAR & Comb. Sci.2009, 6-7, 664. doi： 10.1002/qsar.200860117