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

(国民核生化灾害防护国家重点实验室，北京 102205)

氨基甲酸酯类农药 (carbamates, CMs) 由于具有成本低廉且高效广谱等特点，已被广泛用作农业杀虫剂和除草剂[1]。然而此类农药的不合理滥用，导致其大量残留于生态环境中，给整个生态系统的安全造成了严重威胁[2-6]。因此，建立其生物活性与急性毒性的评估机制对于预防和降低此类农药的危害和潜在风险具有重要意义。

氨基甲酸酯的活性 (毒性) 主要源自于其对生物体内乙酰胆碱酯酶 (acetylcholinesterase, AChE)的可逆共价抑制作用[7]，毒理学中用于评估此类抑制作用效力的重要指标之一是酶活性半数抑制浓度 (the half maximal inhibitory concentration, IC50)。然而，传统上依赖于体内/体外实验的IC50值测试不仅耗时耗力，而且越来越难以适应大数据时代下大规模新型农药的快速评估，同时也有悖于3R原则[8]和21世纪毒性测试的愿景[9]。随着先进的人工智能算法的快速发展以及化合物数据库的不断完善[10-11]，计算机模拟预测方法开始出现并日益走向成熟[12-15]。其中定量构效关系 (quantitative structure-activity relationship, QSAR) 经过近60年的持续改进和跨学科突破，已成为目前化合物的理化性质和生物活性建模的最常用方法之一[16-18]。在QSAR研究中，由于量子化学描述符[19-21]及其中基于概念密度泛函理论 (conceptual density functional theory,CDFT) 的反应性指数[22-24]具有不依赖实验、无统计误差、物理意义明确、可解释性强、能够精确描述分子结构、电子结构及反应性等独特优势，在QSAR模型中的应用逐渐增加，且取得了良好的效果[25]。

本文利用数理意义明确的量化描述符以及线性建模方法，选取目前常用的氨基甲酸酯类农药分子及其他氨基甲酸酯类乙酰胆碱酯酶抑制剂分子作为研究对象，构建了预测性和解释性兼备的氨基甲酸酯IC50预测模型，从而为其农药的活性预测和风险评估提供指导。

1 材料与方法

1.1 数据收集与预处理

严格按照Fourches等[26]提出的QSAR建模数据集预处理流程，精选了146个结构多样且具有确定的对苍蝇乙酰胆碱酯酶IC50值的氨基甲酸酯分子作为训练集。将IC50值转换为以10为底的负对数形式pIC50，使其近似服从正态分布，所得活性范围大于4个数量级。所有原始数据均来源于Pubchem数据库，且经过多次人工仔细检查，尽可能地降低了数据的错误率。

1.2 量化参数计算

综合考虑计算稳定性和计算耗时成本[27]，确定量化参数的计算条件为B3LYP泛函、6-311+G(2d, p) 基组、基于密度的溶剂化模型 (solvation model based on density, SMD) 及water环境，即B3LYP/6-311+G(2d, p)/SMD/water。首先，对于数据集中的每一个分子，均利用GaussView 6中的GMMX3.0模块对其结构进行系统的构象搜索和筛选；然后，在B3LYP/6-311G(2d, p)/SMD/water的计算条件下做进一步的几何结构优化和振动频率分析，得到不含虚频的最低能量构象；最后，对最低能量构象的分子结构在相同计算条件下进行单点能的计算，以及其中性分子 (分子的电子数为N)、得到1个电子 (电子数为N+1) 和失去1个电子 (电子数为N−1) 3种条件下的自然键轨道(natural bond orbital, NBO) 分析。上述计算过程均在软件Gaussian16中完成。

1.3 量化描述符生成

基于量子化学及CDFT理论，共生成模型构建所需的53个描述符。其中分子全局描述参数有：分子的垂直电离势(I)、垂直亲和势(A)、电子化学势(μ)、绝对硬度(η)、亲电性指数(ω)以及分子的前线轨道能量等；分子局部描述参数有：收敛的原子福井函数、原子NBO净电荷以及Wiberg键级及其最弱键级的变化率，如得一电子下最易断裂化学键的键级变化率和失一电子下最易断裂化学键的键级变化率等。

分子的垂直电离势I、垂直亲和势A由所研究分子的中性分子、得到1个电子分子和失去1个电子分子的基态能量或分子的前线轨道能量计算获得。而电子化学势μ、绝对硬度η、亲电性指数ω等则是基于垂直电离势、垂直亲和势数学推导公式获得。简缩亲电福井函数、简缩亲核福井函数、双描述符指数、净亲电指数 (multiphilicity descriptor) 等局域反应性指数则是由Gaussian程序计算结果的原子NBO净电荷计算获得[28]。Wiberg键级由Gaussian程序计算输出文件的Wiberg键级矩阵 (Wiberg bond index matrix in the NAO) 中获得，而得失电子最易断裂化学键的键级变化率，则是依据Wiberg键级数值经公式计算获得。描述符生成过程均在本课题组自主研发的量化描述符提取软件Quantum V1.0[29]中完成。

1.4 模型构建

将训练集中化合物的量化描述符与其pIC50值一一对应进行整合，得到下一步用于模型训练的数据集。设置因变量为pIC50，自变量为53个量化描述符，选用线性建模方法——遗传/偏最小二乘法 (genetic/partial least square, G/PLS) 进行QSAR模型构建。其中G/PLS中种群数设为100，循环迭代 (进化) 次数设为5 000，分别设定主成分数为2～5，方程长度为4～10 (其中一项为常数项，即方程中描述符数量为3～9个)。此外，为确保模型具有良好拟合优度的同时具有较好的预测性能，本文还选取留一法 (leave-one-out, LOO)和自举法 (bootstrapping, BS) 对模型性能进行了严格的交叉验证。以上模型构建及相关统计学分析均在Cerius2软件上进行。

2 结果与讨论

2.1 G/PLS方法QSAR建模结果

表1 不同长度及主成分数的QSAR方程统计参数Table 1 Statistic parameters of QSAR equations from different lengths and primary components

考虑到方程长度增加，所涉及到的描述符数量增加，使得模型复杂且可解释性变差。在方程统计学性能接近的情况下，选取方程长度最短的方程作为下一步的预测方程，共得到如表2所示的4个QSAR方程，各方程中pIC50的预测值与实验值之间的散点图如图1 所示。

表2 QSAR预测方程Table 2 Predictive QSAR equations

2.2 QSAR预测方程讨论与分析

G/PLS计算结果表明，在设定不同主成分及方程长度的条件下，共有以下6个描述符被用作预测方程构建：

Electrophilicity_index_D：采用基于前线轨道理论推导而不是基于基态能量推导的亲电指数。反映分子得失电子后的动态反应性能。

fn_min1：比较分子亲核简缩福井函数最小值，反映在研究体系中，不同分子中的各个原子最小的亲核性能。

delta_f_max1：简缩双描述符最大值，反映分子中原子的亲电和亲核的净值的最大值，分子中的各个原子，实际上一般都是具有双重性，即：既有亲电性又有亲核性。简缩双描述符大于零，则亲电性强；简缩双描述符小于零，则亲核性强。

f0_max1：比较分子自由基简缩福井函数最大值，反映在研究体系中，不同的分子中的各个原子的最大的自由基反应性能。

Max_W1：最大的净亲电指数，反映分子中亲电能力最强的原子的化学反应性能。

Quadrupole：中性分子四级距，反映二维空间中的电荷分布情况。

由分析方程1～4可见：Max_W1及delta_f_max1与pIC50呈正相关，表明氨基甲酸酯对AChE活性的抑制能力可能与其分子中亲电性最强原子的反应性有关，且反应性越强，其抑制能力就越强；反映分子中原子亲核性强弱的描述符fn_min1则与pIC50之间存在负相关关系，即亲核性越强，pIC50反而越小，对AChE活性的抑制能力就越弱。此外，f0_max1与pIC50之间呈现显著的正相关关系，说明氨基甲酸酯对AChE活性的抑制能力还与分子中各原子的最大自由基反应性能有关，且抑制能力随着自由基反应性能的增强而增强。在分子的全局反应性层次，表征分子整体亲电性强弱的亲电性指数Electrophilicity_index_D及表征中性分子极性的Quadrupole则与pIC50呈负相关，由此可推断出，分子整体的亲电性和极性越强，其对AChE活性的抑制能力就越弱。

2.3 一致性模型构建

为平衡各方程之间的预测误差，将各个预测方程进行算数平均，得到以下预测方程：

在此基础上构建了一致性预测模型，并将其作为最终的预测模型。利用一致性模型得到的pIC50预测值与实验值之间的散点图如图2所示，其中预测值与实验值之间的R2为0.823，均方根误差 (RMSE)为0.369。

2.4 模型外部测试

为确保所取得的模型具有良好的预测能力，本文对未参与模型训练的12个氨基甲酸酯分子进行pIC50的外部预测，其预测结果见表3。可以看出，模型的外部预测残差均在一个数量级以内，表明所得到的一致性模型具有较好的预测能力。

表3 外部测试集预测结果Table 3 Predicted results of the external test set

3 小结

本研究严格遵循经济合作与发展组织 (OECD)关于QSAR模型构建的五项原则，构建了氨基甲酸酯类农药及其他氨基甲酸酯类乙酰胆碱酯酶抑制剂分子对苍蝇乙酰胆碱酯酶抑制活性的定量构效关系模型。所有预测方程均经过严格的留一法和自举法交叉验证，其实验值与预测值之间的决定系数R2均在0.8以上，交叉验证R2均在0.7以上。良好的统计学参数结果表明，本文所建立的预测模型不仅具有良好的拟合优度，也具备较好的预测性能。此外，采用量子化学描述符及其中基于CDFT提出的反应性指数进行目标化合物分子表征，一方面更为准确地反映了目标分子化学反应性的结构特征，另一方面也大大增强了所得模型的解释性。研究结果表明：氨基甲酸酯对AChE的抑制活性可能与其分子亲电性及原子自由基反应性有关，且随着二者的增强而增强；此外，分子整体的亲电性和极性也可能影响氨基甲酸酯的AChE抑制活性，且对AChE活性的抑制能力随二者的增强而减弱。

基于量化反应性指数构建定量构效关系模型，对氨基甲酸酯类乙酰胆碱酯酶可逆抑制剂进行毒理评估具有科学可行性，且可大大减少用于毒理学试验的动物数量，符合3R原则和对21世纪毒性测试的美好愿景。研究结果有望为探索和发现氨基甲酸酯类农药与乙酰胆碱酯酶的作用机制提供新思路，可在一定程度上指导设计并合成安全高效的此类新农药；同时，所构建模型具有良好的预测性能，可有效预估此类危险化合物的生物活性，为其风险评估和监管决策提供理论依据，从而有助于降低此类农药使用过程中对生态系统造成的潜在风险和危害。