基于概念密度泛函理論磷酸酯類反應性物質毒性預測

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

1　引言

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

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

2　理論與計算方法

2.1　理論

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

一階響應函數：可以獲得化學勢 μ和電子密度ρ等。

化學勢：

電子密度：

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

二階響應函數：可以獲得體系硬度、軟度、福井函數等。

整體硬度：

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

整體軟度：

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

親電指數：

局域硬度20：

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

局域軟度：

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

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

福井函數：

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

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

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

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

自由基攻擊福井函數：

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

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

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

2.2　計算方法

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

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

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

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

3　結果與討論

3.1　分子反應性分析

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

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

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

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

3.2　氮原子不同質子化狀態及構象對分子反應性影響

3.2.1質子化狀態對分子反應性影響

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

3.2.2構象對分子反應性影響

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

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

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

3.3　反應性指數的QSPR研究

3.3.1QSPR方程構建

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

方程13 (B3LYP/gas)：

方程14 (B3LYP/water)：

方程15 (MP2/gas)：

方程16 (MP2/water)：

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

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

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

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

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

3.3.2QSPR方程的pLD50預測

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

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

表5　根據CDFT和傳統2D-QSPR方法，方程組對訓練集和測試集化合物pLD50預測結果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.

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

4　結論

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

概念密度泛函方法提出的物理量具有明確的物理意義，能很好描述化學反應性的結構特點，因此越來越受到理論化學與(藥物)應用化學家的廣泛關注與重視。由于磷酸酯類反應性物質的高毒性，使得實驗的危險系數高，數據資料不全或難以獲得。所以采用計算機虛擬評估，應用構建的高精度 QSPR模型，預測未知化合物的毒性，是一個很好的選擇。概念密度泛函理論在材料、藥物化學等領域的QSA(T)R應用，僅在最近幾年才略有文獻報道，并有持續增長的趨勢。本文將最新發展的量子力學概念密度泛函理論方法，應用于這類特殊反應性化合物的毒理性能評估和環境毒理評價，為減少動物的使用量，節約資源，降低環境污染等做一些努力。可望為探索磷酸酯類不可逆抑制劑與乙酰膽堿酯酶的作用機制開拓新思路，并為評估該類化合物的生態和環境毒性提供理論依據。

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