基于量化反應(yīng)性指數(shù)的氨基甲酸酯類農(nóng)藥定量構(gòu)效關(guān)系研究

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

(國(guó)民核生化災(zāi)害防護(hù)國(guó)家重點(diǎn)實(shí)驗(yàn)室，北京 102205)

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

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

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

1 材料與方法

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

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

1.2 量化參數(shù)計(jì)算

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

1.3 量化描述符生成

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

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

1.4 模型構(gòu)建

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

2 結(jié)果與討論

2.1 G/PLS方法QSAR建模結(jié)果

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

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

表2 QSAR預(yù)測(cè)方程Table 2 Predictive QSAR equations

2.2 QSAR預(yù)測(cè)方程討論與分析

G/PLS計(jì)算結(jié)果表明，在設(shè)定不同主成分及方程長(zhǎng)度的條件下，共有以下6個(gè)描述符被用作預(yù)測(cè)方程構(gòu)建：

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

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

delta_f_max1：簡(jiǎn)縮雙描述符最大值，反映分子中原子的親電和親核的凈值的最大值，分子中的各個(gè)原子，實(shí)際上一般都是具有雙重性，即：既有親電性又有親核性。簡(jiǎn)縮雙描述符大于零，則親電性強(qiáng)；簡(jiǎn)縮雙描述符小于零，則親核性強(qiáng)。

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

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

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

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

2.3 一致性模型構(gòu)建

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

在此基礎(chǔ)上構(gòu)建了一致性預(yù)測(cè)模型，并將其作為最終的預(yù)測(cè)模型。利用一致性模型得到的pIC50預(yù)測(cè)值與實(shí)驗(yàn)值之間的散點(diǎn)圖如圖2所示，其中預(yù)測(cè)值與實(shí)驗(yàn)值之間的R2為0.823，均方根誤差 (RMSE)為0.369。

2.4 模型外部測(cè)試

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

表3 外部測(cè)試集預(yù)測(cè)結(jié)果Table 3 Predicted results of the external test set

3 小結(jié)

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

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