抗乳腺癌候选药物的优化模型

董奕鑫，张欢欢，王昌会，陈昊，李孝诚

抗乳腺癌候选药物的优化模型

董奕鑫，张欢欢，王昌会，陈昊，李孝诚

（淮北师范大学 数学科学学院，安徽 淮北 235000）

根据华为杯中国研究生数学建模竞赛D题所提供的ERα拮抗剂信息，综合运用灰色关联度分析、BP神经网络、决策树、回归模型等方法和理论，借助MATLAB，SPSS，GeoGebra等软件，构建了化合物生物活性的定量预测模型和ADMET性质分类预测模型，并在此基础上建立抗乳腺癌候选药物的优化模型，筛选出适合用于抗乳腺癌候选药物的化合物．经检验发现，模型均具有良好的性能，可将其应用于虚拟药物筛选流程，为计算机辅助药物设计与药物发现提供参考．

抗乳腺癌候选药品；灰色关联度分析；BP神经网络；决策树分类预测模型

本文研究的问题引自2021年华为杯中国研究生数学建模竞赛D题[1]．乳腺癌是目前世界上最常见、致死率较高的癌症之一．乳腺癌的发展与雌激素受体密切相关，有关研究发现，雌激素受体亚型（ER）在乳腺发育过程中扮演了十分重要的角色[2]．因此，ER被认为是治疗乳腺癌的重要靶标，能够拮抗ERα活性的化合物可能是治疗乳腺癌的候选药物[3]．一个化合物想要成为候选药物，不仅需要具备良好的生物活性（此处指抗乳腺癌活性），还需要在人体内具备良好的药代动力学性质和安全性，合称为ADMET（Absorption吸收、Distribution分布、Metabolism代谢、Excretion排泄、Toxicity毒性）性质[4-5]．

目前，对于活性化合物的筛选，通常采用建立化合物活性预测模型的方法．即针对与疾病相关的某个靶标（此处为ER），收集一系列作用于该靶标的化合物及其生物活性数据，以一系列分子结构描述符作为自变量，化合物的生物活性值作为因变量，构建化合物的定量结构-活性关系（Quantitative Structure-Activity Relationship，QSAR）模型[6]，进而使用该模型预测具有更好生物活性的新化合物分子，或者指导已有活性化合物的结构优化．

然而，当前并没有准确的方法能够筛选出同时具备良好的生物活性和ADMET性质的化合物．为解决此问题，本文建立了化合物生物活性的定量预测模型和ADMET性质分类预测模型，并根据这2个模型建立化合物的最优预测模型，筛选合适的化合物作为治疗乳腺癌症候选药物．具体筛选过程分为4步：第1步，构建相关性分析模型，分别计算各分子描述符与生物活性的关系，找出相关性最大的若干个变量；第2步，选择影响生物活性最显著的前20个分子描述符作为自变量，构建化合物对ER生物活性的定量预测模型；第3步，根据化合物的ADMET数据，分别构建化合物性质的分类预测模型；第4步，建立最优化预测模型，寻找最优化合物的分子描述符，以及这些分子描述符的取值或取值范围，以此来确定抗乳腺癌候选药物．

本文根据“华为杯”数学建模竞赛D题所提供的ER拮抗剂信息（1 974个化合物样本；每个化合物都有729个分子描述符变量，1生物活性的值，5个ADMET性质（表现Caco-2，CYP3A4，hERG，HOB，MN对应值），以分子描述符变作为自变量，生物活性数据作为因变量，同时根据其ADMET性质，构建相关模型．

1 基于灰色关联度分析的自变量筛选

为找出对生物活性最具有显著影响的自变量，需要构建相关性分析模型来分别计算自变量与生物活性的关系．灰色关联分析方法是衡量因素间关联程度的一种方法，用于寻求系统中各子系统（或因素）之间的数值关系，能够为一个系统发展变化态势提供量化的度量，非常适合动态历程分析[7]．因此，本文选择灰色关联度分析法计算反映生物活性与不同自变量之间贴近程度的关联度，通过比较各关联度的大小来判断自变量对生物活性的影响程度[8]．

关联系数是因变量列与自变量列在各个化合物对象的关联程度值，关联程度值不止一个，将各个化合物对象的关联系数集中为一个值，即求其平均值[10]，具体公式为

将729个自变量对因变量的关联度按照从大到小的顺序排列起来．

根据灰色关联度的模型，利用MATLAB软件编程，求得与生物活性灰色关联度最高的20个自变量（见图1）.

图1 前20个因素与生物活性的灰色关联度

2 基于BP神经网络的生物活性预测

利用灰色关联分析模型初步筛选出与生物活性关联度较高的20个自变量后（后续研究均基于此），需要构建化合物对ER生物活性的定量预测模型．由于BP神经网络能够很好地对非线性模型进行预测，因此选择使用BP神经网络构建生物预测模型[11]．

BP网络的基本结构包括1个输入层，1个输出层，1个或多个隐含层（或称为隐层）[12]．基于BP神经网络的基本结构以及算法流程，基于已有的数据，把包含1 974个化合物的20个自变量值以及pIC50值的数据矩阵作为训练集，训练和构建BP神经网络[13]．

该网络的结构是一个多层前馈神经网络，输入维数为20，输出维数为1，具体结构参数设置如下：

（1）层数选择．Kolmogorov定理指出[14-15]，只要不限制隐含层节点数，单隐含层的BP神经网络就可以实现任意的非线性映射．因此，单隐含层的三层BP神经网络可以满足本模型的要求．

（2）各层节点数．根据BP神经网络节点选择的要求，确定输入层、输出层、隐含层的节点数分别为20，1，15．

（3）函数选择．基于模型构建的需要，选取logsig函数、train函数以及S型函数分别作为模型的传递函数、训练函数以及神经元转换函数．

（4）初始参数设定．基于网络需要，设置迭代次数为2 000次，学习速率的初始值为0.01，训练结束的目标精度为0.1．

依据BP神经网络的基本结构，构建了生物活性预测评价模型的BP神经网络结构（见图2） ．

图2 BP神经网络结构

BP神经网络训练参数见图3．利用MATLAB对训练过程中的误差进行分析，训练过程拟合度分析见图4．

图3 BP神经网络训练参数

图4 训练过程拟合度分析

由图3可以看出，BP神经网络训练模型共进行21次迭代（验证集误差不再降低，为防止过度拟合，利用early stop原则，模型训练停止），模型梯度值为0.024 685，且模型在迭代中趋向优化.

由图4可以看出，BP神经网络训练过程中的拟合系数值均在0.8左右，除个别异常点之外，训练样本点几乎均在拟合直线上或者附近呈对称分布，表明BP神经网络模型训练的效果理想，可用于预测其他化合物pIC50的值．

为方便药物筛选，需要将求出的pIC50转化为IC50_nM值．借助GeoGebra软件[16]对训练集中1 974个化合物结构式的IC50_nM列及对应的pIC50列进行回归拟合分析，拟合效果见图5．

图5 IC50_nM与pIC50拟合分析

所得对数函数的拟合模型为

将BP神经网络得到的pIC50预测值代入式（3），便可对其他化合物的IC50_nM进行预测．

利用建立的神经网络生物预测模型，对50种新化合物对应的pIC50值进行预测，结果见表1．

表1 预测集的IC50值和pIC50值

3 决策树分类预测

为找出具备ADMET性质的化合物，基于729个自变量，针对1 974个化合物的ADMET数据，分别构建５个化合物的分类预测模型．由于ADMET数据是由0和1组成，且属于一个因变量多个自变量的分类预测问题，因此采用SPSS构建深度学习的CHAID算法决策树分类预测模型[18]．

由于样本足够大，可将研究数据分为训练数据集和验证数据集，并使用前者构建决策树模型，后者决定树的适合大小，以获得最优模型[19]．决策树分类预测模型先对1 974个化合物进行训练并检验获得分类情况，再导入50个预测化合物的729个自变量，根据所得的分类规则获得预测结果．

利用SPSS进行训练数据分析，利用分割样本验证，其中训练样本与检验样本分别占比70%与30%．在决策树框中，将待预测化合物性质的量作为因变量，将影响ADMET的729个分子描述符量作为自变量．选择CHAID生长法，利用SPSS生成决策树，并输出预测结果．

以HOB为例，利用ADMET性质的决策树分类预测模型研究HOB与自变量的关系．决策树框中，将HOB作为因变量，729个分子描述符作为自变量．选择CHAID生长法，利用SPSS生成决策树，并输出预测结果．基于决策树模型HOB被分为12类（见图6）．

图6 HOB性质的分类预测模型

根据SPSS生成的决策树表（见表2）可知，化合物的HOB性质主要依据BCUTc-1l，maxdO，maxHCsatu，VP-7，ETA_Beta_ns_d，minHBa这6个自变量进行分类（由于版面限制，拆分值只保留小数点后4位数）．

表2 HOB性质分类树

类似地，通过改变决策树框中的因变量，可得出其余4个ADMET性质的分类情况（见表3）．

表3 ADMET性质的分类情况

基于所构建的5个ADMET性质的分类预测模型与分类规则，导入表1中对应的50个化合物的729个分自变量，获得其各对应的5个性质的预测结果（见表4）．

表4 50个化合物的预测结果

4 最优化综合模型

根据提供的ER拮抗剂信息，借助灰色关联度分析、BP神经网络、回归模型以及决策树模型等构建一个包含化合物生物活性的定量预测模型和ADMET性质的分类预测模型．候选药物的筛选模型需要将二者结合，形成最优化综合模型[20]，模型建立过程见图7．

图7 最优化综合模型的建立

优化模型需要筛选出ADMET中有3个及3个以上最好性质的化合物．以灰色关联模型选出的20个与生物活性关联度最大的分子描述符作为自变量，以pIC50值为目标函数，借助BP神经网络模型建立新的关系模型，以求出最优．

关于ADMET性质，本文采用二分类法提供相应的取值．为方便计算，对CYP3A4、MN中的0和1进行替换．替换后，５个ADMET性质中，1均代表该化合物具有好的药物性质，0代表该化合物的具有不好药物性质．若求至少3个性质较好的化合物，利用sum函数，对５个性质的取值求和，和大于3即符合要求．经筛选，共645组化合物满足性质要求，形成新的数据集．

将645组数据作为新的数据集训练网络，其中600组作为新的神经网络训练集数据，其余45组作为预测集数据，获得新的神经网络数据．网络训练后的测试样本的预测值（新的BP神经网络产生）与期望值（题目中给出的数据）非常接近（见图8），经过迭代，达到最佳性能0.001 259 4．

根据BP神经网络中训练过程的误差对比（见图9）可知，预测值与期望值误差较小，表明BP神经网络模型训练的效果理想，可以基于此网络筛选分子描述符．

图8 均方误差随训练次数的变化

图9 误差对比

利用Fminsearch函数计算生物活性达到最大值时各自变量的取值范围，结果见表5．

表5 分子描述符及其取值范围

基于建立抗乳腺癌候选药物的优化模型，经过分析，建议选用包含20个分子描述符的化合物作为抗乳腺癌候选药物．

5 结语

本文针对华为杯中国研究生数学建模竞赛D题中的抗乳腺癌候选药物等问题，通过综合运用灰色关联度分析、BP神经网络、决策树与回归模型等方法和理论建立模型，此模型在药物生物活性预测及ADMET分类预测上均取得良好性能，较好地解决了候选药物的优化问题．在构建模型的过程中，通过对数据扩增以及预测算法模型的迭代优化等方式进一步增强了预测工具的通用性．可将该模型进行推广，应用于虚拟药物筛选流程，为计算机辅助药物设计与药物发现提供新思路，具有较好的借鉴意义．

[1] 中国学位与研究生教育学会．华为杯中国研究生数学建模竞赛D题[EB/OL]．（2021-09-01）[2021-10-16]．https://cp

ipc.acge.org.cn//pw/detail/2c90800c7c2f10dc017c34baa9180cdd.

[2] 路珩，张一奇．雄激素受体在雌激素受体阳性乳腺癌患者中的表达及其临床意义[J]．中国现代医学杂志，2021，31（18）：55-59.

[3] Pizon M，Lux D，Pachmann U，etal．Influence of endocrine therapy on the ratio of androgen receptor（AR）to estrogen receptor（ER）positive circulating epithelial tumor cells（CETCs）in breast cancer[J]．Journal of translational medicine，2018，16（1）：356-364．

[4] 张翠锋，谢海棠，潘国宇．大分子药物的吸收、分布、代谢、排泄和毒性特征及药代模型的应用[J]．药学学报，2016，51（8）：1202-1208．

[5] Feinberg E N，Joshi E，Pande V S，et al．Improvement in ADMET Prediction with Multitask Deep Featurization[J]．Journal of medicinal chemistry，2020，63（16）：8835-8848．

[6] 刘雅红，贺利民，梁智斌，等．用于预测化合物活性的两级拟合QSAR模型的构建方法：中国：102930113B[P]．（2015-02-03）[2021-10-16]．https://wenku.baidu.com/view/4c597ce9ce1755270722192e453610661fd95ac3?fr=xueshu_top．

[7] 邓聚龙．灰理论基础[M]．武汉：华中科技大学出版社，2002．

[8] 虞晓芬，傅玳．多指标综合评价方法综述[J]．统计与决策，2004（11）：119-121．

[9] 罗党，刘思峰．灰色关联决策方法研究[J]．中国管理科学，2005（1）：102-107．

[10] 韩中庚．数学建模方法及其应用[M]．北京：高等教育出版社，2005．

[11] 谢良旭，李峰，谢建平，等．基于融合神经网络模型的药物分子性质预测[J]．计算机科学，2021，48（9）：251-256．

[12] 潘斌．数学建模教程[M]．北京：化学工业出版社，2017．

[13] Jiang Dejun，Lei Tailong，Wang Zhe，et al．ADMET evaluation in drug discovery 20 prediction of breast cancer resistance protein inhibition through machine learning[J]．Journal of Cheminformatics，2020，12（1）：603-617．

[14] Hecht-Nielsen R．Nearest matched filter classification of spatiotemporal patterns[J]．Applied optics，1987，26（10）：1892-1899．

[15] Hecht-Nielsen R．Counter propagation networks[J]．Applied optics，1987，26（23）：4979-4983．

[16] 吴纯良．基于GeoGebra的统计教学课例赏析：两个变量的线性相关（第2课时）[J]．数学通报，2016，55（12）：20-23．

[17] 黄忠裕．初等数学模型[M]．北京：科学出版社，2013．

[18] 李琳，陈德钊，束志恒，等．基于预处理的决策树在化学数据挖掘中的应用[J]．分析化学，2005（8）：1091-1094．

[19] 郭晓龙，蒋艳，邱路．决策树分类模型预测蛋白质相互作用的应用研究[J]．生物医学工程学杂志，2013，30（5）：952-956．

[20] 顾耀文，张博文，郑思，等．基于图注意力网络的药物ADMET分类预测模型构建方法[J]．数据分析与知识发现，2021，5（8）：76-85．

Optimized model of anti-breast cancer candidate drugs

DONG Yixin，ZHANG Huanhuan，WANG Changhui，CHEN Hao，LI Xiaocheng

（School of Mathematical Sciences，Huaibei Normal University，Huaibei 235000，China）

The research data were obtained from the information of ERantagonists provided by the D problem of Huawei Cup，a mathematical modeling competition for graduate students in China．The quantitative prediction model for the biological activity of compounds and the ADMET property classification prediction model were constructed by combining the methods and theories of gray correlation analysis，BP neural network，decision tree and regression model with the help of MATLAB，SPSS，GeoGebra．On basis of it，an optimized model of anti-breast cancer candidate drugs was established，and compounds suitable for anti-breast cancer candidate drugs were screened out．After testing，it is found that the models have good performance and can be applied to the virtual drug screening process to provide reference for computer-aided drug design and drug discovery.

anti-breast cancer candidate drug；gray correlation analysis；BP neural network；decision tree classification prediction model

O22

A

10.3969/j.issn.1007-9831.2022.06.006

1007-9831（2022）06-0030-08

2021-11-20

安徽省自然科学研究项目（1908085MF186）；安徽省高校自然科学研究重点项目（KJ2019A0589）；安徽省质量工程项目

（2020jyxm1670，2020jxtd）

董奕鑫（1998-），女，山东临沂人，在读硕士研究生，从事数学教学评价、数学建模研究．E-mail：dongyixin1998@163.com