武汉大学学报(理学版) 2016, Vol. 62 Issue (3): 293-298
0

文章信息

陈艳 , 王思华 , 李靖 . 2016
CHEN Yan, WANG Sihua, LI Jing . 2016
昆虫酚氧化酶抑制剂的活性预测模型
Prediction Model of the Activity of Insect Henoloxidase Inhibitors
武汉大学学报(理学版), 2016, 62(3): 293-298
Journal of Wuhan University(Natural Science Edition), 2016, 62(3): 293-298
http://dx.doi.org/10.14188/j.1671-8836.2016.03.014

文章历史

收稿日期:2014-11-27
昆虫酚氧化酶抑制剂的活性预测模型
陈艳, 王思华, 李靖    
徐州工程学院 化学化工学院, 江苏 徐州 221111
摘要: 采用电性拓扑状态指数(En)表征昆虫酚氧化酶(PO)抑制剂的分子结构,通过最佳变量子集回归的方法建立了57种PO抑制剂抑制活性(pIC50)的多元线性回归模型,非交叉相关系数和交叉相关系数分别为0.920和0.908,经Jackknife和变异膨胀因子(VIF)检验具有良好的稳定性和预测能力.该模型显示影响PO抑制剂抑制活性的主要因素是-OH,-O-和C=O等分子结构片段.以模型中的3个参数E13E14E16为人工神经网络输入层,设定3∶6∶1的网络结构构建人工神经网络的BP算法模型,相关系数达到0.988.结果表明,与多元线性回归模型相比,BP人工神经网络模型的相关性和预测能力均有较大的提高.
关键词: 电性拓扑状态指数     酚氧化酶抑制剂     人工神经网络     定量结构-活性相关    
Prediction Model of the Activity of Insect Henoloxidase Inhibitors
CHEN Yan, WANG Sihua, LI Jing    
School of Chemistry and Chemical Engineering, Xuzhou Institute of Technology, Xuzhou 221111, Jiangsu,China
Abstract: Electrotopological state index (En)was used to describe the molecular structure of 57 insect henoloxidase inhibitors in this paper. The multiple liner regression(MLR) model was constructed by leaps-and-bounds regression (LBR), and the traditional correlation coefficient and the cross-validation correlation coefficient are 0.920 and 0.908, respectively.The model exhibited excellent stability and predictability evaluated by Jackknife and variance inflation factors (VIF) method. The model shows that the dominant influencing factors of inhibited activity are the molecular structure fragments: -OH, -O-, and C=O in PO inhibitors.The three structural parameters E13, E14, E16were used as the input neurons of artificial neural network,and a 3∶6∶1 network architecture is employed. A satisfied model can be constructed with the back-propagation algorithm, the correlation coefficient R2was 0.988. The results show that the back-propagation algorithm of artificial neural network model have higher correlation and predictability than multiple liner regression model. 
Key words: electrotopological state index     insect henoloxidase inhibitor     artificial neural network     QSAR    
0 引 言

随着科学技术的不断进步,害虫的防治理念已从“杀死”过渡到“控制”,杀死害虫的化学农药长期使用会带来许多问题,例如由于昆虫产生抗药性会导致农药使用剂量不断增加,在粮食、水果、蔬菜中产生残留,所以农药的研发发展到了生物合理设计(biorationaldesign)的阶段.

酚氧化酶(phenoloxidase,PO)是生物体内一类含铜的氧化酶,广泛存在于无脊椎动物、脊椎动物、植物、细菌和真菌等生物体内[1].PO在昆虫的正常发育过程中具有重要的生理功能,是生物体内黑色素合成的关键酶,在昆虫体内参与表皮的硬化、黑化过程,对外来侵染物的免疫防御反应以及伤口的愈合反应,对昆虫的变态发育和免疫系统都起着重要的作用[2].因此,研究PO在昆虫体中的作用及调控方式、在免疫系统中的功能以及抑制机理将有助于增进人们对它的认识深度,进而开发以该酶为靶标的新型农药来控制或减轻害虫的危害[3].

本文使用电性拓扑状态指数(En)为分子描述符表征了57个昆虫酚氧化酶抑制剂的分子结构,建立了PO抑制剂抑制活性pIC50(IC50为酶活性被抑制一半时抑制剂的浓度)的两种预测模型:一是采用最佳变量子集回归的方法确立的三元线性回归方程,二是采用人工神经网络(artificial neural network,ANN)[4]的预测模型,即以最佳变量子集方法筛选出的3个结构参数作为输入层单元,采用3∶6∶1的网络结构建立的BP模型,并对两种模型的相关性、稳定性和预测能力进行了对比.

1 研究方法 1.1 数据来源

57个PO抑制剂化合物的母体结构见图1[5],其分子结构及相应的pIC50[5]列于表1.

图 1 PO抑制剂化合物的分子结构 Figure 1 The molecular structure of PO inhibitors
表1 PO抑制剂抑制活性的实验值与预测值 Table 1 Predicted valuesvs experimental values of PO inhibitors
序号取代基E13E14 E16 pIC50
Exp.[5]Prep.1Err.1Prep.2Err.2
1te0.0000.00010.0062.232.840.612.22-0.01
2tr2-OH8.8790.00010.0672.382.490.112.28-0.10
3va3-OH8.7850.00010.0482.022.500.482.210.19
4tr4-OH8.7360.00010.0362.512.510.002.18-0.33
5tr2,4-OH17.6570.00010.0972.312.16-0.152.500.19
6te3,4-OH17.5960.00010.0772.052.170.122.02-0.03
7tr2,5-OH17.6910.00010.1082.982.15-0.832.82-0.16
8va2-OMe0.0004.90510.3060.972.101.130.81-0.16
9tr3-OMe0.0004.89810.2191.372.130.761.790.42
10tr4-OMe0.0004.89810.1652.722.15-0.572.47-0.25
11te2-OH;3-OMe9.2334.77910.2802.241.78-0.462.260.02
12tr2-OH;4-OMe9.1114.81910.2262.701.80-0.902.68-0.02
13va3-OH;4-OMe9.1394.77910.2061.601.810.211.600.00
14tr4-OH;3-OMe9.0904.77910.2491.301.800.501.29-0.01
15tr4-Propyl-i0.0000.00010.2883.072.74-0.332.97-0.10
16te4-Butyl0.0000.00010.3323.132.72-0.413.12-0.01
17tr4-Butyl-t0.0000.00010.3673.142.71-0.433.250.11
18va8.3850.00010.2011.852.460.612.200.35
19tr2-OH17.3060.00010.2621.472.110.641.720.25
20tr3-OH17.1970.00010.2421.792.120.331.74-0.05
21te4-OH17.1400.00010.2311.792.130.341.790.00
22tr2,4-OH26.1020.00010.2911.521.780.261.50-0.02
23va2,5-OH26.1440.00010.3031.491.770.281.45-0.04
24tr3,4-OH26.0260.00010.2721.531.790.261.590.06
25tr3-OH;4-NH217.4090.00010.3131.592.090.501.910.32
26te4-OH;3-NH217.3700.00010.3272.422.09-0.331.91-0.51
27tr3-OMe8.5494.84210.4141.881.75-0.131.85-0.03
28va4-OMe8.5084.85610.3602.261.77-0.492.290.03
29tr3-OH;4-OMe17.6734.73710.4011.321.440.121.360.04
30tr4-OH;3-OMe17.6574.72210.4431.161.430.271.12-0.04
31te4-Propyl8.6060.00010.4702.762.36-0.402.67-0.09
32tr4-Propyl-i8.6140.00010.4822.722.35-0.372.69-0.03
33va4-Butyl8.6520.00010.5262.982.34-0.642.78-0.20
34tr4-Butyl-t8.6790.00010.5622.682.32-0.362.860.18
35tr4-Hexyl8.7220.00010.6112.932.30-0.632.970.04
36te0.0000.0000.0006.086.370.296.330.25
37tr2-OH9.2930.0000.0005.666.030.375.920.26
38va3-OH9.0890.0000.0006.286.04-0.246.15-0.13
39tr4-OH8.9690.0000.0006.366.04-0.326.28-0.08
40tr2,4-OH18.3040.0000.0005.985.70-0.285.990.01
41te3,4-OH18.1320.0000.0006.545.71-0.836.540.00
42tr2,5-OH18.4080.0000.0005.655.700.055.650.00
43va2-OMe0.0005.1200.0004.675.711.044.65-0.02
44tr3-OMe0.0005.0500.0005.725.720.005.70-0.02
45tr4-OMe0.0005.0110.0006.225.73-0.496.270.05
46te2-OH;3-OMe9.6464.9310.0005.685.38-0.305.820.14
47tr3-OH;4-OMe9.4424.8920.0006.595.40-1.195.97-0.62
48va4-OH9.3234.9310.0005.485.39-0.095.950.47
49tr2,4-OMe0.00010.2560.0004.815.050.245.190.38
50tr3,4-OMe0.00010.2550.0005.575.05-0.525.19-0.38
51te2,4,5-OMe0.00015.5870.0004.004.360.364.000.00
52tr4-Propyl-i0.0000.0000.0006.466.37-0.096.33-0.13
53va4-Butyl0.0000.0000.0006.556.37-0.186.33-0.22
54tr4-Butyl-t0.0000.0000.0006.256.370.126.330.08
55tr3,5-Butyl-i; 2-OH10.4510.0000.0004.145.991.854.09-0.05
56te4-N(CH3)20.0000.0000.0006.446.37-0.076.33-0.11
57tr-0.0000.0000.0006.196.370.186.330.14
注:tr(训练集),va(验证集),te(测试集),Exp.(实验值),Prep.(预测值),Err.(误差)
1.2 电性拓扑状态指数(En)的建构方法

电性拓扑状态指数是由Kier和Hall所建立的表征分子中每个非氢原子(n)的结构参数[6],由两部分组成,一是非氢原子本身的原子结构及局部拓扑环境,称为固有状态值(本征值),该数值与非氢原子所处的具体分子环境无关;二是反映该原子受分子中其他非氢原子扰动程度的本征值增量(ΔIn).

分子中每个非氢原子的固有状态值被定义为[6]

(1)

式中,Nn为原子n的电子层数;δn=σnhn ,其中δn,σn为原子n的支化度及形成σ键的键数,hn为原子n之间键合的氢原子数;δνn原为原子n的Kier[7]点价,后被Hall等[6, 8]修正为:

(2)

式中,πnkn分别为π键数及孤对电子数.

Kier和Hall在文献[6]中对常见化合物中各种不同的结构片段指定了不同类型,共有79种.本文57个PO抑制剂化合物所涉及的原子类型及其固有状态值见表2.

表2 PO抑制剂分子中原子类型的符号及固有状态值 Table 2 The symbols and intrinsic state values of atoms type in PO inhibitors
编号结构片段符号InEn编号结构片段符号InEn
1-CH3sCH32.000E19—OHsOH6.000E13
2-CH2ssCH21.500E210—O—ssO3.500E14
3>CH—sssCH1.333E311O=dO7.000E16
4>C<ssssC1.250E412—NH2sNH24.000E17
5—CH=sdCH2.000E613—NH—ssNH2.500E18
6—CH=aaCH2.000E714>N—sssN2.000E20
7>C=ssdC1.667E815—N=sdN3.000E23
8>C=aaCa1.667E916S=dS3.667E33
注:s(单键),d(双键),a(芳香环中的共轭键)

分子中每个成键原子都不是孤立的,要受到周围其他原子的干扰,而使其电性发生改变,此改变量即为增量:

(3)

(3)式中,rnj为原子n与原子j之间相隔的最短路径数加1,是对分子中原子n以外的其他原子求和.

原子n的电性拓扑状态指数按下式计算:

(4)

(4)式中,t为具有相同固有状态值的原子的个数.

本文在Matlab环境下,利用程序[9, 10]计算了57个PO抑制剂分子16个电性拓扑状态指数.

2 pIC50的多元线性回归模型 2.1 模型的建立

将16种电性拓扑状态指数作为结构参数与57个PO抑制剂分子的抑制活性pIC50输入Minitab软件系统,经最佳变量子集回归建立的多元线性回归(multiple liner regression,MLR)模型见表3.

表3 pIC50En的最佳变量子集回归结果 Table 3 Leaps-and-bounds regression for PI50 as a function of En
序号RR2R2adjR2CVSF变量
10.9340.8710.8680.8610.694372.771E16
20.9500.9020.8990.8920.610249.914E16 ,E14
30.9590.9200.9160.9080.557203.908E16 ,E14,E13
40.9610.9240.9180.8930.548158.442E16 ,E14,E13,E3
n(样本数),R(相关系数),R2(判断系数),R2CV(逐一剔除法的交叉验证系数),S(估计标准误差),F(Fischer检验值).

表3可见,随着模型中变量数增多,R2逐渐增大、相应的S逐渐减小.但其交叉验证的R2CV先是增大至0.908,随后下降至0.893,说明第3个数学模型具有最好的稳健性和预测能力.该模型为:

用(5)式计算所得的预测值(Prep.1)列于表1,平均误差为0.422.

2.2 稳健性检验 2.2.1 逐一剔除法交叉验证

逐一剔除法的交叉验证系数RCV2是目前较为广泛使用的模型验证方法之一,从表3可见,所选最优模型((5)式)的交叉验证系数R2CV为0.908,符合目前公认的R2CV≥0.5的标准[11],同时该模型的R2CVR2(0.920)略小,交叉验证标准偏差SCV(0.597)比S(0.557)略大,说明该模型的稳定性和预测能力都比较理想.另外,所建模型的R2adjR2CV相差为0.008,远小于0.3[11],说明该模型没有过拟合,不存在不相关的其他变量或离域点.

2.2.2 Jackknife验证

采用Jackknife法[12]对所建模型进行稳健性检验,由于样本超过30个,所以本文采用逐组剔除法(把化合物按个位数从0到9分组,共分成10组,每次从57个化合物中剔除1组,剩余的9组根据(5)式进行回归)得到的10个相关系数(R2)的分布是:0.917以下2次,0.917~0.923之间6次,0.923以上2次,呈良好的正态分布,平均值为0.921,与所建模型((5)式)的0.920非常接近,说明该模型中不存在异常数据和随机相关,具有总体可接受的稳健性.

2.2.3 模型中自变量的相关性检验

用变异膨胀因子(variance inflation factors,VIF)评价模型中各自变量间的多重相关性,VIF的定义式为[13]

(6)

(6)式中,R2为模型中一个自变量与余下自变量的相关系数.当VIF=1,表明模型中各自变量间不存在相关性;VIF<5时,各自变量间相关性很弱,该模型是稳定可以接受的;VIF>5则存在明显的相关性,所建模型不能用于估算和预测[13].所建模型中E13E14E16的VIF值分别为1.241,1.148和1.152,说明所建模型中各自变量间不存在自相关,模型具有可接受的稳定性.

2.3 构效分析

进入模型的E13是羟基中氧原子的电拓扑指数、E14是-O-中氧原子的拓扑指数、E16是C=O中氧原子的拓扑指数,这3个指数均和氧原子有关.氧原子电负性高且含有孤对电子,可以和昆虫酚氧化酶中心周围的基团如—SH,—NH2,—OH中的活性氢以氢键结合,在活性中心周围形成空间位阻从而阻止底物和活性中心作用而发挥抑制剂的功效,所以所建模型揭示了影响昆虫酚氧化酶抑制剂抑制活性的主要因素和抑制剂的抑制机理.

3 PO抑制剂pIC50的BP-ANN模型

人工神经网络(ANN)利用计算机来模拟生物神经网络的某些结构和功能,具有强大的非线性处理、自组织协调和容错能力,已广泛应用于定量构效关系(QSAR)中[14].其中应用最为广泛的是基于误差反向传播(Back-propagation,BP)算法的三层神经网络,即输入层、隐蔽层和输出层.本文采用Matlab提供的神经网络工具箱中的BP算法进行拟合.输入层单元选为最佳变量子集回归得到的最佳变量组合E13E14E16.根据Andrea规则寻找最佳隐蔽层的单元数H

(7)

其中N,M分别是样本数和网络总权重.M被定义为:

(8)

(8)式中:IHQ分别是输入层、隐蔽层和输出层的单元数.本文的I=3,Q=1,N=57,可得4.98<H< 6.13,取H=6.

至此,本文采用3∶6∶1的网络结构建立模型,其中输出层为PO抑制剂的pIC50.为了进一步避免过拟合、过训练,将数据集分成3组:训练集、验证集和测试集(依次在表1中使用tr,va,te上标标注),各集的化合物数分别是34,11,12.在训练过程中,当验证集误差开始上升便立即停止训练,既防止了网络的过训练、过拟合,又减少了训练时间,由此建立的模型为:训练集R=0.994,验证集R=0.994,测试集R=0.996,均与总体的相关系数(R=0.994)比较接近,说明所建模型是稳定的,不存在过训练、过拟合现象.该模型给出的预测值与实验值非常吻合,列于表1(Prep.2),平均误差为0.139.

用两种模型得到的预测值和实验值的相关性见图2,可以看出,用BP-ANN模型得到的预测值更接近实验值,进一步说明了BP-ANN模型具有更强的预测能力.

图 2 用两种模型预测的pIC50的预测值和实验值的相关性 Figure 2 Predicted values vs experimental values of pIC50 by two models
4 结 论

1) 电性拓扑状态指数En较好地表征分子的结构特征,有效解释了影响昆虫酚氧化酶抑制剂抑制活性的本质因素;

2) 从多元线性回归模型的结构参数来看,影响昆虫酚氧化酶抑制剂抑制活性的主要结构片段是—OH,—O—和C=O;

3) 所建多元线性回归模型中各变量的VIF值均小于1.5,说明它们是相互独立的,不存在多重共线性.模型经Jackknife法、逐一剔除法交互检验证明具有良好的稳定型和预测能力;

4) 文献建立了八元线性回归模型[5],预测值和实验值的相对误差为0.252,本研究采用人工神经网络的建模方法,仅用了3个变量作为输入层单元,预测平均误差为0.139,说明本研究方法及预测结果优于文献.

综上,本研究通过两种预测模型的建立,进一步揭示了昆虫酚氧化酶抑制剂的抑制机理,增进了对该酶的认识深度,可望为开发以该酶为靶标的新型农药提供理论依据.

参考文献
[1] 徐亚玲, 李文楚. 昆虫酚氧化酶作用机制的研究进展[J]. 安徽农业科学, 2010 ,38 (27) : 14844 –14846. XU L L, LI W C. Research progress in activation mechanisms of phenoloxide in insects[J]. Journal of Anhui Agricultural Sciences, 2010, 38 (27) : 14844 –14846.
[2] 王树栋. 四种植物源化合物对甜菜夜蛾酚氧化酶的影响及酶免疫学研究[D].泰安:山东农业大学, 2010. WANG S D.Study on the Effects of Four Plant-Derived Compounds on Phenoloxidase of Spodotera exigua ( hübner) and Immunology of the Enzyme[D]. Taian:Shandong Agricultural University, 2010(Ch).
[3] 刘守柱, 杜学林, 戴明勋. 昆虫酚氧化酶特性及其生理学功能[J]. 安徽农业科学, 2010 ,38 (34) : 19249 –19251. LIU S Z, DU X L, DAI M X. Characteristics and physiological functions of henoloxidase from insect[J]. Journal of Anhui Agricultural Sciences, 2010, 38 (34) : 19249 –19251.
[4] 黄俊, 余刚, 王溢磊, 等. 用量子化学方法预测PCDFs的Ah受体结合能力(1)——神经网络模型[J]. 南京理工大学学报, 2003 ,27 (6) : 709 –714. HUANG J, YU G, WANG Y L, et al. Predicting Ah receptor binding data of PCDFs using quantum chemical method (1): Artificial neural network approach[J]. Journal of Nanjing University of Science and Technology, 2003, 27 (6) : 709 –714.
[5] 莫凌云, 谢天宝, 印春生. 57个昆虫酚氧化酶抑制剂的QSAR研究[J]. 桂林理工大学学报, 2011 ,31 (4) : 574 –578. MO L Y, XIE T B, YIN C S. QSAR of 57 insect phenoloxidase inhibitors[J]. Journal of Guilin University of Technology, 2011, 31 (4) : 574 –578.
[6] HALL L H, KIER L B. Molecular similarity based on novel atom-type electrotopological state indices[J]. Journal of Chemical Information and Computer Science, 1995, 35 (6) : 1074 –1080.
[7] KIER L B, HALL L H. Molecular Connectivity in Structure Activity Analysis. England: Research studies press.[M] 1986 .
[8] HALL L H, KIER L B. Electrotopological state indices for atom types: A novel combination of electronic, topological and valence state information[J]. Journal of Chemical Information and Computer Science, 1995, 35 (6) : 1039 –1045.
[9] 胡黔楠, 梁逸曾, 王亚丽, 等. 直观队列命名法的基本原理及其在矩阵与拓扑指数计算中的应用[J]. 计算机与应用化学, 2003 ,20 (4) : 386 –390. HU Q N, LIANG Y Z, WANG Y L, et al. The basic principles of heruisticqueue notation and its applications in calculation of matrix and topological index for topogicalgraghs[J]. Computers and Applied Chemistry, 2003, 20 (4) : 386 –390.
[10] 张婷, 梁逸曾, 赵晨曦, 等. 基于分子结构预测气相色谱程序升温保留指数[J]. 分析化学, 2006 ,34 (11) : 1607 –1610. ZHANG T, LIANG Y Z, ZHAO C X, et al. Prediction of temperature-programmed retention indices from molecule structures[J]. Chinese Journal of Analytical Chemistry, 2006, 34 (11) : 1607 –1610.
[11] SWAMINA S, FUREY W, PLCHER J, et al. Crystalstructure of staphylococcal enterotoxin B, a superan-tigen[J]. Nature, 1992, 359 (6389) : 801 –806.
[12] DIETRICH W S, DREYER N D, HANSCH C. Confidence intervalestimators for parameters associated with quantitative structure/activity relationship[J]. J Med Chem, 1980, 23 : 1201 –1205.
[13] 冯长君. 手性有机酸保留指数的手性指数及原子类型电拓扑指数模型[J]. 物理化学学报, 2010 ,26 (1) : 193 –198. FENG C J. Mathematical model for retention indices, chiral index and electrotopological state indices for atom types of chiral organic acids[J]. Acta Physico-Chimica Sinica, 2010, 26 (1) : 193 –198.
[14] 赵筱萍, 范骁辉, 余杰, 等. 一类基于组效关系神经网络模型的中药药效预测方法[J]. 中国中药杂志, 2004 ,29 (11) : 1082 –1085. ZHAO X P, FAN X H, YU J, et al. A method for predicting activity of traditional Chinese medicine based on quantitative composition-activity relationship of neural network model[J]. Journal of Chinese MateriaMedica, 2004, 29 (11) : 1082 –1085.