双酚A类似物雌激素干扰效应的定量结构-活性关系模型

杨先海，刘会会，杨倩,3，刘济宁,

1. 环境保护部南京环境科学研究所，南京 210042 2. 南京理工大学 环境与生物工程学院江苏省化工污染控制与资源化高校重点实验室，南京 210094 3. 南京工业大学生物与制药工程学院，南京 210009



双酚A类似物雌激素干扰效应的定量结构-活性关系模型

杨先海1，刘会会2，杨倩1,3，刘济宁1,

1. 环境保护部南京环境科学研究所，南京 210042 2. 南京理工大学 环境与生物工程学院江苏省化工污染控制与资源化高校重点实验室，南京 210094 3. 南京工业大学生物与制药工程学院，南京 210009

双酚A类似物；内分泌干扰物；雌激素效应；内分泌干扰效应；定量结构-活性关系

Received 31 March 2016 accepted 5 May 2016

双酚A(bisphenol A, BPA)是全球生产量最大的化学原料之一，被广泛应用于日用品的生产[1]。然而，大量的研究表明，BPA具有较强的内分泌干扰等毒性效应[2-3]，可干扰激素受体[4-5]，进而影响胚胎/胎儿发育[6]，损害雄性生殖健康[7]等。考虑到BPA的潜在危害效应，部分国家已经开始限制使用含BPA的塑料制品，尤其是对婴幼儿等敏感人群使用的制品，许多国家/地区直接采取了禁用政策[8-9]。为了应对管制措施，一些双酚A类似物，如双酚AF(bisphenol AF, BPAF)、双酚S(bisphenol S, BPS)等被开发用来替代含BPA的工业化学品的生产。由于其与双酚A结构相似，双酚A类似物是否也具有潜在内分泌干扰效应就受到越来越多的关注。

研究表明，双酚A类似物不仅具有急、慢性毒性效应[10]，也具有内分泌干扰效应[11-13]，且部分双酚A类似物的雌激素活性和抗雄性激素活性比BPA更强。例如，Matsushima等[14]的体外细胞实验结果显示，BPAF与雌激素受体α(ERα)的亲和力相当于BPA的20倍，与雌激素受体β(ERβ)的亲和力相当于BPA的48倍。Kitamura等[11]也通过荧光素酶报告基因研究了19种BPA类似物在MCF-7细胞中的雌激素活性，发现这类物质在10-8～10-4mol·L-1就能刺激MCF-7细胞的增殖；此外，还发现除四溴双酚A和四氯双酚A外，其他几种物质均具有较强的抗雄性激素活性。Ji等[15]的研究表明，当把雄性和雌性斑马鱼暴露于0.5 μg·L-1～ 50 μg·L-1BPS的水样中21 d后，雌性和雄性斑马鱼血浆中的雌二醇含量显著升高，而雄性斑马鱼血浆中的睾酮含量则显著降低。因此，双酚A类似物的潜在内分干扰效应也应得到足够重视。目前，尽管已有研究表明部分双酚A类似物具有雌激素活性和抗雄性激素活性等内分泌干扰效应，但是主要是体外试验结果为主，考虑到双酚A类似物种类较多，难以对双酚A类似物一一开展整体动物试验研究。为解决该困境，一方面需要新的测试策略，以快速填补缺失的数据，另一方面需要通过设定优先级，筛查高优先级的双酚A类似物，减少进行整体动物测试的化学品数量。

计算毒理学技术在设定优先级及填补数据缺失方面具有一定的优势。在双酚A类似物的内分泌干扰效应研究方面，针对双酚A类似物分子结构与其雌激素干扰效应关系的研究还较少[16-19]。现有的研究主要涉及作用机制研究和构建双酚A类似物的雌激素干扰效应预测模型。Delfosse等[17]分析了BPA及其2种类似物与雌激素受体的作用模式。Zhuang等[18]采用分子动力学模拟方法，研究了BPA及其7种类似物与雌激素受体alpha(ERα)的作用机制，发现静电作用对双酚A类似物的雌激素效应具有重要作用。Coleman等[16]测试了24种双酚A类似物的雌激素效应，并采用比较分子力场法(CoMFA)、比较分子相似性指数(CoMSIA)方法和分子全息QSAR方法(HQSAR)构建了预测模型。该研究未基于构建的模型预测其他双酚A类似物的雌激素效应数据。Ng等[19]基于18个化合物，采用分子对接方法，计算了小分子与雌激素受体的结合能力(得分)，建立了化合物雌激素效应数据与打分之间的关系，并用构建的关系填补了其他27种化合物(其中双酚A类似物17种)的雌激素活性数据。本研究将采用定量结构-活性关系(QSAR)方法，基于24个有雌激素效应的数据，构建双酚A类似物雌激素效应的预测模型，并填补其他双酚A类似物缺失的雌激素效应数据。

1 材料与方法 (Materials and methods)

1.1 数据集

数据集I包含24个双酚A类似物的雌激素效应值。实验数据来源于Coleman等[16]的测定结果。实验数据是采用人乳腺癌细胞(MCF-7细胞)增殖实验获取。实验测定的指标是半数效应浓度(EC50)。此外，本文还收集了另外22个无试验测试数据的双酚A类似物(数据集II)，将用本研究构建的模型预测这部分化合物的雌激素效应值。所有化合物信息及其效应值列于表1。

QSAR建模的指标是取对数的相对效应势(logRP)，其定义为：

(1)

其中，EC50,E2和EC50,BPA分别是17β-雌二醇(E2)和双酚A类似物的半数效应浓度。

1.2 分子结构优化

首先采用ChemBioOffice 2010软件生成初始的化合物分子结构，进而根据上述分子结构生成Gaussian输入文件，然后运用Gaussian 09优化化合物分子结构[20]。具体而言，采用Gaussian 09程序包中的密度泛函理论(DFT)中的B3LYP方法，在6-31+g(d,p)基组水平上进行结构优化，获得其稳定的分子构象。对优化好的几何构型进行频率分析，以确保体系无虚频。在上述几何优化和DFT的计算中都使用自洽场(SCRF)和积分方程形式极化连续介质模型(IEFPCM)来考虑溶剂(水)效应[21-22]。

1.3 分子描述符计算

其次，基于Gaussian 09优化的结构，使用Dragon 6.0软件计算描述符[23]。Dragon 6.0包括29大类4 885种分子描述符，是一个全面且高效的分子描述符计算软件。基于如下标准，对计算的4 885种描述符进行了初步筛选：去掉常数和近似常数的描述符，去掉至少有一个缺失值的描述符。最终，描述符集包含914种描述符。

表1 双酚A及其类似物信息、雌激素效应数据

注：*验证集化合物。

Note: * Compounds selected as the external validation set.

1.4 QSAR模型构建与表征

(2)

此外，还计算了均方根误差(RMSE)和模型的平均绝对误差(MAE)，以评估预测方法的可靠性。其中，MAE计算式如下[26]：

(3)

其中，logRP (实验值)和logRP(预测值)分别是实验测定和预测的logRP值，n是化合物数。

采用欧几里德距离图和基于杠杆值(leverage)的Williams图定义模型的应用域。欧几里德距离图使用AMBIT Discovery v0.04 (http://ambit.sourceforge.net/download_ ambitdiscovery.html)软件计算。

Williams图是由标准残差(δ*)和杠杆值(以hi表示，i代表不同的化合物)定义的一种QSAR模型应用域[24]。其中，δ*使用下式计算：

(4)

训练集化合物的H值计算式为：

H = X(XTX)-1XT

(5)

其中，X为n × m的矩阵，n为训练集化合物的个数，m为模型中预测变量的个数，X矩阵表征了模型中训练集化合物的描述符空间。

训练集和验证集中化合物的hi值可由式6计算得到：

hi= xiT(XTX)-1xi

(6)

其中，xi为训练集和验证集中第i个化合物预测变量的行向量。

警戒值(h*)定义式为：

h*= 3(m + 1)/n

(7)

其中，m为预测变量的个数。

由于使用SPSS逐步回归时，可能得到一系列模型，本研究基于上述的模型表征参数、方差膨胀因子(VIF)值及描述符个数要求等指标，选取最优模型。

2 结果与讨论(Results and discussion)

2.1 最优模型及其表征结果

最优模型为：

logRP = -10.7-2.53Mor17u+64.1Vs,max

(8)

2.2 应用域表征

在OECD关于QSAR构建与验证导则中[24]，明确指出要对所建模型进行应用域表征。模型应用域表征结果如图1和图2所示。模型定义的欧几里德距离小于0.748。训练集和验证集中化合物都在模型应用域中，说明模型训练集具有较好的代表性。此外，在Williams图中(图2)，若化合物的标准残差δ*落在(-3.0, +3.0)以外时，认为该点是离群点。如图可见，所有化合物的标准残差δ*均落于(-3.0, +3.0)以内，且所有化合物的杠杆值均小于警戒值h*，说明模型没有离群点。因此，该模型能够用于预测应用域内其他双酚A类似物的抗雄激素效应值。

图1 基于欧几里德距离方法表征的模型应用域Fig. 1 Characterization of application domain for developed model based on the Euclidean distance

图2 模型Williams图Fig. 2 Williams plot of developed model

图3 静电势分布图Fig. 3 Electrostatic potentials on molecular surface

图4 数据集II中化合物的适用性评估Fig. 4 Applicability evaluation of the compounds in the data set II

2.3 描述符意义与机理解释

所建模型包含2个描述符，即Mor17u和Vs,max，Mor17u是3D-MoRSE类描述符，Vs,max是量化描述符。其中，Mor17u是未加权的3D-MoRSE描述符；3D-MoRSE(3D molecule representation of structures based on electron diffraction)描述符是基于电子衍射表征的分子三维结构描述符[29]，表征了官能团的影响。与双酚A相比，分子结构中引入吸电子取代基团会使Mor17u值增加，而引入供电子取代基团会使Mor17u值减小。例如，当双酚A羟基邻位被溴或氯原子取代时，四溴双酚A和四氯双酚A的Mor17u值从双酚A的-1.072，增加至-0.559和-0.319；而当双酚A羟基邻位被1或2个甲基取代时，2,2-双(4-羟基-3,5-二甲基苯基)丙烷和4,4'-(1-甲基亚乙基)双(2-甲基苯酚)的Mor17u值从双酚A的-1.072，减小至-1.273和-1.479。方程中Mor17u为负系数(r = -0.80, P< 0.0001)，说明化合物logRP值随Mor17u值的减小而增加。因此，当双酚A羟基邻位引入吸电子取代基团会减弱分子的雌激素效应，而引入供电子取代基团则会增强分子的雌激素效应。

图5 基于模型预测的数据集II中化合物logRP值Fig. 5 Predicted logRP values of compounds in data set II based on developed model

Vs,max是分子表面上的最正静电势。图3显示了分子静电势分布，对于含酚羟基的化合物(如图3中双酚A)，分子中静电势最正的区域是酚羟基对应的部分。前人的研究表明，与雌激素受体形成氢键相互作用、静电作用是小分子表现雌激素效应的关键[18,21-22]。对双酚A类似物而言，酚羟基是其与雌激素受体形成氢键相互作用、静电作用的核心基团[18]。因此，该描述符表征了化合物分子与雌激素受体之间的氢键相互作用、静电作用。从图4可以看出，酚羟基邻位取代基也会影响分子表面正静电势的分布。与双酚A相比，分子结构中引入吸电子取代基团会使分子表面正静电势减小。在方程中Vs,max为正系数(r = 0.75, P< 0.0001)，说明化合物logRP值随Vs,max值的增加而增加。因此，当双酚A酚羟基邻位引入吸电子取代基团会减弱分子的雌激素效应。综上，模型选出的2个变量都解释了取代基对双酚A类似物雌激素效应的影响。

2.4 预测其他双酚A类似物的雌激素效应值

本文采用构建的模型预测其他22个双酚A类似物的雌激素效应数据。首先计算了模型对这22个化合物的适用性。计算结果如图4所示，有3个化合物欧几里德距离大于0.748，处于模型应用域外，3个化合物分别是2,2-二苯基丙烷、二苯基甲烷和双酚A双(二苯基磷酸酯)。由于模型训练集化合物分子结构中都含有酚羟基，因而2,2-二苯基丙烷和二苯基甲烷分子结构中不含酚羟基基团是导致其处于模型应用域外的原因。对双酚A双(二苯基磷酸酯)而言，其分子结构中含2个二苯基磷酸酯基，与训练集化合物分子结构差异较大。综上，模型对其余19个化合物适用。

模型对数据集II中化合物的预测值列于表1。此外，前人的研究结果表明，在该实验体系下，17β-雌二醇(E2)的半数效应浓度(EC50)约为10-11mol·L-1。根据E2的EC50值和式(1)，计算了数据集I中化合物实测的EC50值和数据集II中化合物预测的EC50值(表1)。图5显示了预测的效应值分布范围及其对应的化合物数量。结果表明，数据集II中化合物logRP值介于-8 ～ -2之间，而数据集I(建模数据集)其logRP值介于-9 ～ -3之间。数据集II中仅有3个双酚A类似物的预测logRP值介于-3 ～ -2，超出了建模数据集化合物的logRP值范围。但是由于这3个化合物在模型应用域内，说明预测值可信。

从logRP值的定义可以看出，其值越大，表示对应化合物的雌激素效应越强。从表1可知，有14个双酚A类似物的logRP值大于双酚A。根据EC50值可以看出，数据集II中有7个化合物的EC50值仅是双酚A类似物EC50值的二十分之一。在将来的研究中，需进一步试验确认其雌激素效应。

对于不在模型应用域内的3个化合物，模型预测表明，2,2-二苯基丙烷、二苯基甲烷和双酚A双(二苯基磷酸酯)的logRP分别为-7.02、-7.56和-3.03，即表明预测的2,2-二苯基丙烷和二苯基甲烷雌激素效应较弱，而双酚A双(二苯基磷酸酯)雌激素效应较强。根据上述分析，分子结构中不含酚羟基基团是导致2,2-二苯基丙烷和二苯基甲烷处于模型应用域外的原因。由于分子结构中酚羟基基团的缺失，导致2,2-二苯基丙烷和二苯基甲烷不能与雌激素受体形成氢键相互作用。因此，2,2-二苯基丙烷和二苯基甲烷雌激素效应较弱的预测结果可信度较高。对于双酚A双(二苯基磷酸酯)而言，由于其分子结构与训练集化合物差异最大，且分子中不含酚羟基等关键基团。在模型训练集中，同为酯类的双酚A双(二甲基酯)的雌激素效应较弱(logRP = -7.36)，因此，双酚A双(二苯基磷酸酯)的雌激素效应预测数据不确定性较大。

2.5 模型比较

现有双酚A类似物相关的预测模型情况列于表2。与表2中的模型相比，仅本研究所建预测模型严格按照OECD提出的QSAR模型构建与验证准则进行了模型拟合优度、稳健性、预测能力及应用域的表征，且还进行了机理解释。因此，本研究所建预测模型能够较好地预测应用域内其他双酚A类似物的雌激素效应数据。

表2 不同模型之间的比较

注：alogRP是取对数的相对效应势，其定义为式(1)；logRBA是取对数的相对亲和势，其定义为：

其中，IC50,E2和IC50,BPAs分别是17β-雌二醇(E2)和双酚A类似物的半数抑制浓度。

c原文未计算验证集表征参数，表中数据是根据原文验证集的实验效应数据和预测效应数据计算(原文表5)。

Note:alogRP is the logarithm of the relative competing potency and is defined by equation 1; logRBA is logarithm of relative binding affinity, which is defined as:

where IC50,E2and IC50,BPAsare the concentrations of 17β-estradiol and bisphenols at 50% inhibition of 17β-estradiol binding from estrogen receptor, respectively.

cThe statisticsis calculated by employing the observed and predicted data of the compounds in validation set (Table 5 in the original paper) becuase of the statistics for validation set is not calculated in the original paper.

[1] Eladak S, Grisin T, Moison D, et al. A new chapter in the bisphenol A story: Bisphenol S and bisphenol F are not safe alternatives to this compound [J]. Fertility and Sterility, 2015, 103(1): 11-21

[2] Rochester J R. Bisphenol A and human health: A review of the literature [J]. Reproductive Toxicology, 2013, 42: 132-155

[3] Glausiusz J. Toxicology: The plastics puzzle [J]. Nature, 2014, 508(7496): 306-308

[4] Sui Y, Ai N, Park S H, et al. Bisphenol A and its analogues activate human pregnane X receptor [J]. Environmental Health Perspectives, 2012, 120(3): 399-405

[5] Delfosse V, Grimaldi M, le Maire A, et al. Nuclear receptor profiling of bisphenol-A and its halogenated analogues [J]. Vitamins & Hormones, 2014, 94: 229-251

[6] Virginie R F, René H, Gabriel L. Effects of endocrine disruptors on the human fetal testis [J]. Annales Dendocrinologie, 2014, 75(2): 54-57

[7] Li D K, Zhou Z J, Miao M, et al. Urine bisphenol-A (BPA) level in relation to semen quality [J]. Fertility and Sterility, 2011, 95(2): 625-630

[8] Commission Directive 2011/8/EU of 28 January 2011 amending Directive 2002/72/EC as regards the restriction of use of Bisphenol A in plastic infant feeding bottles [EB/OL]. (2011-01-28) [2015-12-17]. http://faolex.fao.org/docs/pdf/eur100741.pdf

[9] Health Canada. Archived — Order Amending Schedule I to the Hazardous Products Act (bisphenol A) [EB/OL]. (2010-03-31) [2015-12-17]. http://www.gazette.gc.ca/rp-pr/p2/2010/2010-03-31/html/sor-dors53-eng.html

[10] Chen M Y, Ike M, Fujita M. Acute toxicity, mutagenicity, and estrogenicity of bisphenol-A and other bisphenols [J]. Environmental Toxicology, 2002, 17(1): 80-86

[11] Kitamura S, Suzuki T, Sanoh S, et al. Comparative study of the endocrine-disrupting activity of bisphenol A and 19 related compounds [J]. Toxicological Sciences, 2005, 84(2): 249-259

[12] Riu A, Grimaldi M, le Maire A, et al. Peroxisome proliferator-activated receptor γ is a target for halogenated analogs of bisphenol A [J]. Environmental Health Perspectives, 2011, 119(9): 1227-1232

[13] Molina-Molina J M, Amaya E, Grimaldi M, et al. In vitro study on the agonistic and antagonistic activities of bisphenol-S and other bisphenol-A congeners and derivatives via nuclear receptors [J]. Toxicology and Applied Pharmacology, 2013, 272(1): 127-136

[14] Matsushima A, Liu X, Okada H, et al. Bisphenol AF is a full agonist for the estrogen receptor ERα but a highly specific antagonist for ERβ [J]. Environmental Health Perspectives, 2010, 118: 1267

[15] Ji K, Hong S, Kho Y, et al. Effects of bisphenol S exposure on endocrine functions and reproduction of zebrafish [J]. Environmental Science & Technology, 2013, 47(15): 8793-8800

[16] Coleman K P, Toscano Jr W A, Wiese T E. QSAR models of the in vitro estrogen activity of bisphenol A analogs [J]. QSAR & Combinatorial Science, 2003, 22: 78-88

[17] Delfosse V, Grimaldi M, Pons J L, et al. Structural and mechanistic insights into bisphenols action provide guidelines for risk assessment and discovery of bisphenol A substitutes [J]. Proceedings of the National Academy of Sciences of the United States of America, 2012, 109(37): 14930-14935

[18] Zhuang S L, Zhang C L, Liu W P. Atomic insights into distinct hormonal activities of bisphenol A analogues toward PPARγ and ERα receptors [J]. Chemical Research in Toxicology, 2014, 7(10): 1769-1779

[19] Ng H W, Shu M, Luo H, et al. Estrogenic activity data extraction and in silico prediction show the endocrine disruption potential of bisphenol A replacement compounds [J]. Chemical Research in Toxicology, 2015, 28(9): 1784-1795

[20] Frisch M J, Trucks G W, Schlegel H B, et al. Gaussian Inc. Wallingford CT [CP]. 2009

[21] Li F, Li X H, Shao J P, et al. Estrogenic activity of anthraquinone derivatives: In vitro and in silico studies [J]. Chemical Research in Toxicology, 2010, 23(8): 1349-1355

[22] 李斐. 部分有机污染物雌激素效应和甲状腺激素效应的计算模拟与验证[D]. 大连: 大连理工大学, 2010: 56-59

Li F. Computational simulation and validation on the estrogenic and thyroid hormone activities of some organic chemicals [D]. Dalian: Dalian University of Technology, 2010: 56-59 (in Chinese)

[23] Talete srl, DRAGON (Software for Molecular Descriptor Calculation) Version 6.0-2012 [CP]. http://www.talete.mi.it/

[24] Organisation for Economic Co-operation and Development (OECD). Guidance document on the validation of (quantitative) structure activity relationships ((Q)SAR) models, Technical Report for OECD Environment, Health and Safety Publications Series on Testing and Assessment No. 69 [R]. Paris: OECD, 2007

[25] Schüürmann G, Ebert R U, Chen J W, et al. External validation and prediction employing the predictive squared correlation coefficient — Test set activity mean vs training set activity mean [J]. Journal of Chemical Information and Modeling, 2008, 48(11): 2140-2145

[26] Hyndman R J, Koehler A B. Another look at measures of forecast accuracy [J]. International Journal of Forecasting, 2006, 22(4): 679-688

[27] Golbraikh A, Shen M, Xiao Z, et al. Rational selection of training and test sets for the development of validated QSAR models [J]. Journal of Computer-Aided Molecular Design, 2003, 17(2-4): 241-253

[28] Kaneko H, Funatsu K. Criterion for evaluating the predictive ability of nonlinear regression models without cross-validation [J]. Journal of Chemical Information and Modeling, 2013, 53(9): 2341-2348

[29] Schuur J, Gasteiger J. Infrared spectra simulation of substituted benzene derivatives on the basis of a 3D structure representation [J]. Analytical Chemistry, 1997, 69(13): 2398-2405

[30] Wiese T E, Kral L G, Dennis K E, et al. Optimization of estrogen growth response in MCF-7 cells [J]. In Vitro Cellular & Developmental Biology - Animal, 1992, 28A(9-10): 595-602

Predicting Estrogenic Activity of Bisphenols Using Quantitative Structure-Activity Relationships

Yang Xianhai1, Liu Huihui2, Yang Qian1,3, Liu Jining1,*

1. Nanjing Institute of Environmental Science, Ministry of Environmental Protection, Nanjing 210042, China 2. Jiangsu Key Laboratory of Chemical Pollution Control and Resources Reuse, School of Environmental and Biological Engineering, Nanjing University of Science and Technology, Nanjing 210094, China 3. The College of Biotechnology and Pharmaceutical Engineering, Nanjing Tech University, Nanjing 210009, China

bisphenols; endocrine disrupting chemicals; estrogenic activity; endocrine disrupting effects; quantitative structure-activity relationships

国家自然科学基金(No.201507038)；江苏省自然科学基金(No.BK20151100)；中央级公益性科研院所基本科研业务专项(2015)

杨先海(1985-)，男，工学博士,助理研究员，研究方向为环境内分泌干扰效应的计算模拟研究，E-mail: xhyang@nies.org；

*通讯作者(Corresponding author)， E-mail: ljn@nies.org

10.7524/AJE.1673-5897.20160331015

2016-03-31 录用日期：2016-05-05

1673-5897(2016)4-069-10

X171.5

A

简介：刘济宁(1977-)，男，绿色化学工学博士，副研究员，主要研究方向为化学品环境毒理学研究，发表学术论文70余篇。

杨先海, 刘会会, 杨倩, 等. 双酚A类似物雌激素干扰效应的定量结构-活性关系模型[J]. 生态毒理学报，2016, 11(4): 69-78

Yang X H, Liu H H, Yang Q, et al. Predicting estrogenic activity of bisphenols using quantitative structure-activity relationships [J]. Asian Journal of Ecotoxicology, 2016, 11(4): 69-78 (in Chinese)