田家豪, 王亮亮, 潘 里, 丁俊杰, 丁曉琴

(國民核生化災害防護國家重點實驗室，北京 102205)

氨基甲酸酯類農(nóng)藥 (carbamates, CMs) 由于具有成本低廉且高效廣譜等特點，已被廣泛用作農(nóng)業(yè)殺蟲劑和除草劑[1]。然而此類農(nóng)藥的不合理濫用，導致其大量殘留于生態(tài)環(huán)境中，給整個生態(tài)系統(tǒng)的安全造成了嚴重威脅[2-6]。因此，建立其生物活性與急性毒性的評估機制對于預防和降低此類農(nóng)藥的危害和潛在風險具有重要意義。

氨基甲酸酯的活性 (毒性) 主要源自于其對生物體內(nèi)乙酰膽堿酯酶 (acetylcholinesterase, AChE)的可逆共價抑制作用[7]，毒理學中用于評估此類抑制作用效力的重要指標之一是酶活性半數(shù)抑制濃度 (the half maximal inhibitory concentration, IC50)。然而，傳統(tǒng)上依賴于體內(nèi)/體外實驗的IC50值測試不僅耗時耗力，而且越來越難以適應大數(shù)據(jù)時代下大規(guī)模新型農(nóng)藥的快速評估，同時也有悖于3R原則[8]和21世紀毒性測試的愿景[9]。隨著先進的人工智能算法的快速發(fā)展以及化合物數(shù)據(jù)庫的不斷完善[10-11]，計算機模擬預測方法開始出現(xiàn)并日益走向成熟[12-15]。其中定量構效關系 (quantitative structure-activity relationship, QSAR) 經(jīng)過近60年的持續(xù)改進和跨學科突破，已成為目前化合物的理化性質(zhì)和生物活性建模的最常用方法之一[16-18]。在QSAR研究中，由于量子化學描述符[19-21]及其中基于概念密度泛函理論 (conceptual density functional theory,CDFT) 的反應性指數(shù)[22-24]具有不依賴實驗、無統(tǒng)計誤差、物理意義明確、可解釋性強、能夠精確描述分子結構、電子結構及反應性等獨特優(yōu)勢，在QSAR模型中的應用逐漸增加，且取得了良好的效果[25]。

本文利用數(shù)理意義明確的量化描述符以及線性建模方法，選取目前常用的氨基甲酸酯類農(nóng)藥分子及其他氨基甲酸酯類乙酰膽堿酯酶抑制劑分子作為研究對象，構建了預測性和解釋性兼?zhèn)涞陌被姿狨C50預測模型，從而為其農(nóng)藥的活性預測和風險評估提供指導。

1 材料與方法

1.1 數(shù)據(jù)收集與預處理

嚴格按照Fourches等[26]提出的QSAR建模數(shù)據(jù)集預處理流程，精選了146個結構多樣且具有確定的對蒼蠅乙酰膽堿酯酶IC50值的氨基甲酸酯分子作為訓練集。將IC50值轉換為以10為底的負對數(shù)形式pIC50，使其近似服從正態(tài)分布，所得活性范圍大于4個數(shù)量級。所有原始數(shù)據(jù)均來源于Pubchem數(shù)據(jù)庫，且經(jīng)過多次人工仔細檢查，盡可能地降低了數(shù)據(jù)的錯誤率。

1.2 量化參數(shù)計算

綜合考慮計算穩(wěn)定性和計算耗時成本[27]，確定量化參數(shù)的計算條件為B3LYP泛函、6-311+G(2d, p) 基組、基于密度的溶劑化模型 (solvation model based on density, SMD) 及water環(huán)境，即B3LYP/6-311+G(2d, p)/SMD/water。首先，對于數(shù)據(jù)集中的每一個分子，均利用GaussView 6中的GMMX3.0模塊對其結構進行系統(tǒng)的構象搜索和篩選；然后，在B3LYP/6-311G(2d, p)/SMD/water的計算條件下做進一步的幾何結構優(yōu)化和振動頻率分析，得到不含虛頻的最低能量構象；最后，對最低能量構象的分子結構在相同計算條件下進行單點能的計算，以及其中性分子 (分子的電子數(shù)為N)、得到1個電子 (電子數(shù)為N+1) 和失去1個電子 (電子數(shù)為N?1) 3種條件下的自然鍵軌道(natural bond orbital, NBO) 分析。上述計算過程均在軟件Gaussian16中完成。

1.3 量化描述符生成

基于量子化學及CDFT理論，共生成模型構建所需的53個描述符。其中分子全局描述參數(shù)有：分子的垂直電離勢(I)、垂直親和勢(A)、電子化學勢(μ)、絕對硬度(η)、親電性指數(shù)(ω)以及分子的前線軌道能量等；分子局部描述參數(shù)有：收斂的原子福井函數(shù)、原子NBO凈電荷以及Wiberg鍵級及其最弱鍵級的變化率，如得一電子下最易斷裂化學鍵的鍵級變化率和失一電子下最易斷裂化學鍵的鍵級變化率等。

分子的垂直電離勢I、垂直親和勢A由所研究分子的中性分子、得到1個電子分子和失去1個電子分子的基態(tài)能量或分子的前線軌道能量計算獲得。而電子化學勢μ、絕對硬度η、親電性指數(shù)ω等則是基于垂直電離勢、垂直親和勢數(shù)學推導公式獲得。簡縮親電福井函數(shù)、簡縮親核福井函數(shù)、雙描述符指數(shù)、凈親電指數(shù) (multiphilicity descriptor) 等局域反應性指數(shù)則是由Gaussian程序計算結果的原子NBO凈電荷計算獲得[28]。Wiberg鍵級由Gaussian程序計算輸出文件的Wiberg鍵級矩陣 (Wiberg bond index matrix in the NAO) 中獲得，而得失電子最易斷裂化學鍵的鍵級變化率，則是依據(jù)Wiberg鍵級數(shù)值經(jīng)公式計算獲得。描述符生成過程均在本課題組自主研發(fā)的量化描述符提取軟件Quantum V1.0[29]中完成。

1.4 模型構建

將訓練集中化合物的量化描述符與其pIC50值一一對應進行整合，得到下一步用于模型訓練的數(shù)據(jù)集。設置因變量為pIC50，自變量為53個量化描述符，選用線性建模方法——遺傳/偏最小二乘法 (genetic/partial least square, G/PLS) 進行QSAR模型構建。其中G/PLS中種群數(shù)設為100，循環(huán)迭代 (進化) 次數(shù)設為5 000，分別設定主成分數(shù)為2～5，方程長度為4～10 (其中一項為常數(shù)項，即方程中描述符數(shù)量為3～9個)。此外，為確保模型具有良好擬合優(yōu)度的同時具有較好的預測性能，本文還選取留一法 (leave-one-out, LOO)和自舉法 (bootstrapping, BS) 對模型性能進行了嚴格的交叉驗證。以上模型構建及相關統(tǒng)計學分析均在Cerius2軟件上進行。

2 結果與討論

2.1 G/PLS方法QSAR建模結果

表1 不同長度及主成分數(shù)的QSAR方程統(tǒng)計參數(shù)Table 1 Statistic parameters of QSAR equations from different lengths and primary components

考慮到方程長度增加，所涉及到的描述符數(shù)量增加，使得模型復雜且可解釋性變差。在方程統(tǒng)計學性能接近的情況下，選取方程長度最短的方程作為下一步的預測方程，共得到如表2所示的4個QSAR方程，各方程中pIC50的預測值與實驗值之間的散點圖如圖1 所示。

表2 QSAR預測方程Table 2 Predictive QSAR equations

2.2 QSAR預測方程討論與分析

G/PLS計算結果表明，在設定不同主成分及方程長度的條件下，共有以下6個描述符被用作預測方程構建：

Electrophilicity_index_D：采用基于前線軌道理論推導而不是基于基態(tài)能量推導的親電指數(shù)。反映分子得失電子后的動態(tài)反應性能。

fn_min1：比較分子親核簡縮福井函數(shù)最小值，反映在研究體系中，不同分子中的各個原子最小的親核性能。

delta_f_max1：簡縮雙描述符最大值，反映分子中原子的親電和親核的凈值的最大值，分子中的各個原子，實際上一般都是具有雙重性，即：既有親電性又有親核性。簡縮雙描述符大于零，則親電性強；簡縮雙描述符小于零，則親核性強。

f0_max1：比較分子自由基簡縮福井函數(shù)最大值，反映在研究體系中，不同的分子中的各個原子的最大的自由基反應性能。

Max_W1：最大的凈親電指數(shù)，反映分子中親電能力最強的原子的化學反應性能。

Quadrupole：中性分子四級距，反映二維空間中的電荷分布情況。

由分析方程1～4可見：Max_W1及delta_f_max1與pIC50呈正相關，表明氨基甲酸酯對AChE活性的抑制能力可能與其分子中親電性最強原子的反應性有關，且反應性越強，其抑制能力就越強；反映分子中原子親核性強弱的描述符fn_min1則與pIC50之間存在負相關關系，即親核性越強，pIC50反而越小，對AChE活性的抑制能力就越弱。此外，f0_max1與pIC50之間呈現(xiàn)顯著的正相關關系，說明氨基甲酸酯對AChE活性的抑制能力還與分子中各原子的最大自由基反應性能有關，且抑制能力隨著自由基反應性能的增強而增強。在分子的全局反應性層次，表征分子整體親電性強弱的親電性指數(shù)Electrophilicity_index_D及表征中性分子極性的Quadrupole則與pIC50呈負相關，由此可推斷出，分子整體的親電性和極性越強，其對AChE活性的抑制能力就越弱。

2.3 一致性模型構建

為平衡各方程之間的預測誤差，將各個預測方程進行算數(shù)平均，得到以下預測方程：

在此基礎上構建了一致性預測模型，并將其作為最終的預測模型。利用一致性模型得到的pIC50預測值與實驗值之間的散點圖如圖2所示，其中預測值與實驗值之間的R2為0.823，均方根誤差 (RMSE)為0.369。

2.4 模型外部測試

為確保所取得的模型具有良好的預測能力，本文對未參與模型訓練的12個氨基甲酸酯分子進行pIC50的外部預測，其預測結果見表3?？梢钥闯?，模型的外部預測殘差均在一個數(shù)量級以內(nèi)，表明所得到的一致性模型具有較好的預測能力。

表3 外部測試集預測結果Table 3 Predicted results of the external test set

3 小結

本研究嚴格遵循經(jīng)濟合作與發(fā)展組織 (OECD)關于QSAR模型構建的五項原則，構建了氨基甲酸酯類農(nóng)藥及其他氨基甲酸酯類乙酰膽堿酯酶抑制劑分子對蒼蠅乙酰膽堿酯酶抑制活性的定量構效關系模型。所有預測方程均經(jīng)過嚴格的留一法和自舉法交叉驗證，其實驗值與預測值之間的決定系數(shù)R2均在0.8以上，交叉驗證R2均在0.7以上。良好的統(tǒng)計學參數(shù)結果表明，本文所建立的預測模型不僅具有良好的擬合優(yōu)度，也具備較好的預測性能。此外，采用量子化學描述符及其中基于CDFT提出的反應性指數(shù)進行目標化合物分子表征，一方面更為準確地反映了目標分子化學反應性的結構特征，另一方面也大大增強了所得模型的解釋性。研究結果表明：氨基甲酸酯對AChE的抑制活性可能與其分子親電性及原子自由基反應性有關，且隨著二者的增強而增強；此外，分子整體的親電性和極性也可能影響氨基甲酸酯的AChE抑制活性，且對AChE活性的抑制能力隨二者的增強而減弱。

基于量化反應性指數(shù)構建定量構效關系模型，對氨基甲酸酯類乙酰膽堿酯酶可逆抑制劑進行毒理評估具有科學可行性，且可大大減少用于毒理學試驗的動物數(shù)量，符合3R原則和對21世紀毒性測試的美好愿景。研究結果有望為探索和發(fā)現(xiàn)氨基甲酸酯類農(nóng)藥與乙酰膽堿酯酶的作用機制提供新思路，可在一定程度上指導設計并合成安全高效的此類新農(nóng)藥；同時，所構建模型具有良好的預測性能，可有效預估此類危險化合物的生物活性，為其風險評估和監(jiān)管決策提供理論依據(jù)，從而有助于降低此類農(nóng)藥使用過程中對生態(tài)系統(tǒng)造成的潛在風險和危害。