董奕鑫，張歡歡，王昌會，陳昊，李孝誠

抗乳腺癌候選藥物的優(yōu)化模型

董奕鑫，張歡歡，王昌會，陳昊，李孝誠

（淮北師范大學(xué) 數(shù)學(xué)科學(xué)學(xué)院，安徽 淮北 235000）

根據(jù)華為杯中國研究生數(shù)學(xué)建模競賽D題所提供的ERα拮抗劑信息，綜合運用灰色關(guān)聯(lián)度分析、BP神經(jīng)網(wǎng)絡(luò)、決策樹、回歸模型等方法和理論，借助MATLAB，SPSS，GeoGebra等軟件，構(gòu)建了化合物生物活性的定量預(yù)測模型和ADMET性質(zhì)分類預(yù)測模型，并在此基礎(chǔ)上建立抗乳腺癌候選藥物的優(yōu)化模型，篩選出適合用于抗乳腺癌候選藥物的化合物．經(jīng)檢驗發(fā)現(xiàn)，模型均具有良好的性能，可將其應(yīng)用于虛擬藥物篩選流程，為計算機輔助藥物設(shè)計與藥物發(fā)現(xiàn)提供參考．

抗乳腺癌候選藥品；灰色關(guān)聯(lián)度分析；BP神經(jīng)網(wǎng)絡(luò)；決策樹分類預(yù)測模型

本文研究的問題引自2021年華為杯中國研究生數(shù)學(xué)建模競賽D題[1]．乳腺癌是目前世界上最常見、致死率較高的癌癥之一．乳腺癌的發(fā)展與雌激素受體密切相關(guān)，有關(guān)研究發(fā)現(xiàn)，雌激素受體亞型（ER）在乳腺發(fā)育過程中扮演了十分重要的角色[2]．因此，ER被認(rèn)為是治療乳腺癌的重要靶標(biāo)，能夠拮抗ERα活性的化合物可能是治療乳腺癌的候選藥物[3]．一個化合物想要成為候選藥物，不僅需要具備良好的生物活性（此處指抗乳腺癌活性），還需要在人體內(nèi)具備良好的藥代動力學(xué)性質(zhì)和安全性，合稱為ADMET（Absorption吸收、Distribution分布、Metabolism代謝、Excretion排泄、Toxicity毒性）性質(zhì)[4-5]．

目前，對于活性化合物的篩選，通常采用建立化合物活性預(yù)測模型的方法．即針對與疾病相關(guān)的某個靶標(biāo)（此處為ER），收集一系列作用于該靶標(biāo)的化合物及其生物活性數(shù)據(jù)，以一系列分子結(jié)構(gòu)描述符作為自變量，化合物的生物活性值作為因變量，構(gòu)建化合物的定量結(jié)構(gòu)-活性關(guān)系（Quantitative Structure-Activity Relationship，QSAR）模型[6]，進而使用該模型預(yù)測具有更好生物活性的新化合物分子，或者指導(dǎo)已有活性化合物的結(jié)構(gòu)優(yōu)化．

然而，當(dāng)前并沒有準(zhǔn)確的方法能夠篩選出同時具備良好的生物活性和ADMET性質(zhì)的化合物．為解決此問題，本文建立了化合物生物活性的定量預(yù)測模型和ADMET性質(zhì)分類預(yù)測模型，并根據(jù)這2個模型建立化合物的最優(yōu)預(yù)測模型，篩選合適的化合物作為治療乳腺癌癥候選藥物．具體篩選過程分為4步：第1步，構(gòu)建相關(guān)性分析模型，分別計算各分子描述符與生物活性的關(guān)系，找出相關(guān)性最大的若干個變量；第2步，選擇影響生物活性最顯著的前20個分子描述符作為自變量，構(gòu)建化合物對ER生物活性的定量預(yù)測模型；第3步，根據(jù)化合物的ADMET數(shù)據(jù)，分別構(gòu)建化合物性質(zhì)的分類預(yù)測模型；第4步，建立最優(yōu)化預(yù)測模型，尋找最優(yōu)化合物的分子描述符，以及這些分子描述符的取值或取值范圍，以此來確定抗乳腺癌候選藥物．

本文根據(jù)“華為杯”數(shù)學(xué)建模競賽D題所提供的ER拮抗劑信息（1 974個化合物樣本；每個化合物都有729個分子描述符變量，1生物活性的值，5個ADMET性質(zhì)（表現(xiàn)Caco-2，CYP3A4，hERG，HOB，MN對應(yīng)值），以分子描述符變作為自變量，生物活性數(shù)據(jù)作為因變量，同時根據(jù)其ADMET性質(zhì)，構(gòu)建相關(guān)模型．

1 基于灰色關(guān)聯(lián)度分析的自變量篩選

為找出對生物活性最具有顯著影響的自變量，需要構(gòu)建相關(guān)性分析模型來分別計算自變量與生物活性的關(guān)系．灰色關(guān)聯(lián)分析方法是衡量因素間關(guān)聯(lián)程度的一種方法，用于尋求系統(tǒng)中各子系統(tǒng)（或因素）之間的數(shù)值關(guān)系，能夠為一個系統(tǒng)發(fā)展變化態(tài)勢提供量化的度量，非常適合動態(tài)歷程分析[7]．因此，本文選擇灰色關(guān)聯(lián)度分析法計算反映生物活性與不同自變量之間貼近程度的關(guān)聯(lián)度，通過比較各關(guān)聯(lián)度的大小來判斷自變量對生物活性的影響程度[8]．

關(guān)聯(lián)系數(shù)是因變量列與自變量列在各個化合物對象的關(guān)聯(lián)程度值，關(guān)聯(lián)程度值不止一個，將各個化合物對象的關(guān)聯(lián)系數(shù)集中為一個值，即求其平均值[10]，具體公式為

將729個自變量對因變量的關(guān)聯(lián)度按照從大到小的順序排列起來．

根據(jù)灰色關(guān)聯(lián)度的模型，利用MATLAB軟件編程，求得與生物活性灰色關(guān)聯(lián)度最高的20個自變量（見圖1）.

圖1 前20個因素與生物活性的灰色關(guān)聯(lián)度

2 基于BP神經(jīng)網(wǎng)絡(luò)的生物活性預(yù)測

利用灰色關(guān)聯(lián)分析模型初步篩選出與生物活性關(guān)聯(lián)度較高的20個自變量后（后續(xù)研究均基于此），需要構(gòu)建化合物對ER生物活性的定量預(yù)測模型．由于BP神經(jīng)網(wǎng)絡(luò)能夠很好地對非線性模型進行預(yù)測，因此選擇使用BP神經(jīng)網(wǎng)絡(luò)構(gòu)建生物預(yù)測模型[11]．

BP網(wǎng)絡(luò)的基本結(jié)構(gòu)包括1個輸入層，1個輸出層，1個或多個隱含層（或稱為隱層）[12]．基于BP神經(jīng)網(wǎng)絡(luò)的基本結(jié)構(gòu)以及算法流程，基于已有的數(shù)據(jù)，把包含1 974個化合物的20個自變量值以及pIC50值的數(shù)據(jù)矩陣作為訓(xùn)練集，訓(xùn)練和構(gòu)建BP神經(jīng)網(wǎng)絡(luò)[13]．

該網(wǎng)絡(luò)的結(jié)構(gòu)是一個多層前饋神經(jīng)網(wǎng)絡(luò)，輸入維數(shù)為20，輸出維數(shù)為1，具體結(jié)構(gòu)參數(shù)設(shè)置如下：

（1）層數(shù)選擇．Kolmogorov定理指出[14-15]，只要不限制隱含層節(jié)點數(shù)，單隱含層的BP神經(jīng)網(wǎng)絡(luò)就可以實現(xiàn)任意的非線性映射．因此，單隱含層的三層BP神經(jīng)網(wǎng)絡(luò)可以滿足本模型的要求．

（2）各層節(jié)點數(shù)．根據(jù)BP神經(jīng)網(wǎng)絡(luò)節(jié)點選擇的要求，確定輸入層、輸出層、隱含層的節(jié)點數(shù)分別為20，1，15．

（3）函數(shù)選擇．基于模型構(gòu)建的需要，選取logsig函數(shù)、train函數(shù)以及S型函數(shù)分別作為模型的傳遞函數(shù)、訓(xùn)練函數(shù)以及神經(jīng)元轉(zhuǎn)換函數(shù)．

（4）初始參數(shù)設(shè)定．基于網(wǎng)絡(luò)需要，設(shè)置迭代次數(shù)為2 000次，學(xué)習(xí)速率的初始值為0.01，訓(xùn)練結(jié)束的目標(biāo)精度為0.1．

依據(jù)BP神經(jīng)網(wǎng)絡(luò)的基本結(jié)構(gòu)，構(gòu)建了生物活性預(yù)測評價模型的BP神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)（見圖2） ．

圖2 BP神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)

BP神經(jīng)網(wǎng)絡(luò)訓(xùn)練參數(shù)見圖3．利用MATLAB對訓(xùn)練過程中的誤差進行分析，訓(xùn)練過程擬合度分析見圖4．

圖3 BP神經(jīng)網(wǎng)絡(luò)訓(xùn)練參數(shù)

圖4 訓(xùn)練過程擬合度分析

由圖3可以看出，BP神經(jīng)網(wǎng)絡(luò)訓(xùn)練模型共進行21次迭代（驗證集誤差不再降低，為防止過度擬合，利用early stop原則，模型訓(xùn)練停止），模型梯度值為0.024 685，且模型在迭代中趨向優(yōu)化.

由圖4可以看出，BP神經(jīng)網(wǎng)絡(luò)訓(xùn)練過程中的擬合系數(shù)值均在0.8左右，除個別異常點之外，訓(xùn)練樣本點幾乎均在擬合直線上或者附近呈對稱分布，表明BP神經(jīng)網(wǎng)絡(luò)模型訓(xùn)練的效果理想，可用于預(yù)測其他化合物pIC50的值．

為方便藥物篩選，需要將求出的pIC50轉(zhuǎn)化為IC50_nM值．借助GeoGebra軟件[16]對訓(xùn)練集中1 974個化合物結(jié)構(gòu)式的IC50_nM列及對應(yīng)的pIC50列進行回歸擬合分析，擬合效果見圖5．

圖5 IC50_nM與pIC50擬合分析

所得對數(shù)函數(shù)的擬合模型為

將BP神經(jīng)網(wǎng)絡(luò)得到的pIC50預(yù)測值代入式（3），便可對其他化合物的IC50_nM進行預(yù)測．

利用建立的神經(jīng)網(wǎng)絡(luò)生物預(yù)測模型，對50種新化合物對應(yīng)的pIC50值進行預(yù)測，結(jié)果見表1．

表1 預(yù)測集的IC50值和pIC50值

3 決策樹分類預(yù)測

為找出具備ADMET性質(zhì)的化合物，基于729個自變量，針對1 974個化合物的ADMET數(shù)據(jù)，分別構(gòu)建５個化合物的分類預(yù)測模型．由于ADMET數(shù)據(jù)是由0和1組成，且屬于一個因變量多個自變量的分類預(yù)測問題，因此采用SPSS構(gòu)建深度學(xué)習(xí)的CHAID算法決策樹分類預(yù)測模型[18]．

由于樣本足夠大，可將研究數(shù)據(jù)分為訓(xùn)練數(shù)據(jù)集和驗證數(shù)據(jù)集，并使用前者構(gòu)建決策樹模型，后者決定樹的適合大小，以獲得最優(yōu)模型[19]．決策樹分類預(yù)測模型先對1 974個化合物進行訓(xùn)練并檢驗獲得分類情況，再導(dǎo)入50個預(yù)測化合物的729個自變量，根據(jù)所得的分類規(guī)則獲得預(yù)測結(jié)果．

利用SPSS進行訓(xùn)練數(shù)據(jù)分析，利用分割樣本驗證，其中訓(xùn)練樣本與檢驗樣本分別占比70%與30%．在決策樹框中，將待預(yù)測化合物性質(zhì)的量作為因變量，將影響ADMET的729個分子描述符量作為自變量．選擇CHAID生長法，利用SPSS生成決策樹，并輸出預(yù)測結(jié)果．

以HOB為例，利用ADMET性質(zhì)的決策樹分類預(yù)測模型研究HOB與自變量的關(guān)系．決策樹框中，將HOB作為因變量，729個分子描述符作為自變量．選擇CHAID生長法，利用SPSS生成決策樹，并輸出預(yù)測結(jié)果．基于決策樹模型HOB被分為12類（見圖6）．

圖6 HOB性質(zhì)的分類預(yù)測模型

根據(jù)SPSS生成的決策樹表（見表2）可知，化合物的HOB性質(zhì)主要依據(jù)BCUTc-1l，maxdO，maxHCsatu，VP-7，ETA_Beta_ns_d，minHBa這6個自變量進行分類（由于版面限制，拆分值只保留小數(shù)點后4位數(shù)）．

表2 HOB性質(zhì)分類樹

類似地，通過改變決策樹框中的因變量，可得出其余4個ADMET性質(zhì)的分類情況（見表3）．

表3 ADMET性質(zhì)的分類情況

基于所構(gòu)建的5個ADMET性質(zhì)的分類預(yù)測模型與分類規(guī)則，導(dǎo)入表1中對應(yīng)的50個化合物的729個分自變量，獲得其各對應(yīng)的5個性質(zhì)的預(yù)測結(jié)果（見表4）．

表4 50個化合物的預(yù)測結(jié)果

4 最優(yōu)化綜合模型

根據(jù)提供的ER拮抗劑信息，借助灰色關(guān)聯(lián)度分析、BP神經(jīng)網(wǎng)絡(luò)、回歸模型以及決策樹模型等構(gòu)建一個包含化合物生物活性的定量預(yù)測模型和ADMET性質(zhì)的分類預(yù)測模型．候選藥物的篩選模型需要將二者結(jié)合，形成最優(yōu)化綜合模型[20]，模型建立過程見圖7．

圖7 最優(yōu)化綜合模型的建立

優(yōu)化模型需要篩選出ADMET中有3個及3個以上最好性質(zhì)的化合物．以灰色關(guān)聯(lián)模型選出的20個與生物活性關(guān)聯(lián)度最大的分子描述符作為自變量，以pIC50值為目標(biāo)函數(shù)，借助BP神經(jīng)網(wǎng)絡(luò)模型建立新的關(guān)系模型，以求出最優(yōu)．

關(guān)于ADMET性質(zhì)，本文采用二分類法提供相應(yīng)的取值．為方便計算，對CYP3A4、MN中的0和1進行替換．替換后，５個ADMET性質(zhì)中，1均代表該化合物具有好的藥物性質(zhì)，0代表該化合物的具有不好藥物性質(zhì)．若求至少3個性質(zhì)較好的化合物，利用sum函數(shù)，對５個性質(zhì)的取值求和，和大于3即符合要求．經(jīng)篩選，共645組化合物滿足性質(zhì)要求，形成新的數(shù)據(jù)集．

將645組數(shù)據(jù)作為新的數(shù)據(jù)集訓(xùn)練網(wǎng)絡(luò)，其中600組作為新的神經(jīng)網(wǎng)絡(luò)訓(xùn)練集數(shù)據(jù)，其余45組作為預(yù)測集數(shù)據(jù)，獲得新的神經(jīng)網(wǎng)絡(luò)數(shù)據(jù)．網(wǎng)絡(luò)訓(xùn)練后的測試樣本的預(yù)測值（新的BP神經(jīng)網(wǎng)絡(luò)產(chǎn)生）與期望值（題目中給出的數(shù)據(jù)）非常接近（見圖8），經(jīng)過迭代，達(dá)到最佳性能0.001 259 4．

根據(jù)BP神經(jīng)網(wǎng)絡(luò)中訓(xùn)練過程的誤差對比（見圖9）可知，預(yù)測值與期望值誤差較小，表明BP神經(jīng)網(wǎng)絡(luò)模型訓(xùn)練的效果理想，可以基于此網(wǎng)絡(luò)篩選分子描述符．

圖8 均方誤差隨訓(xùn)練次數(shù)的變化

圖9 誤差對比

利用Fminsearch函數(shù)計算生物活性達(dá)到最大值時各自變量的取值范圍，結(jié)果見表5．

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

基于建立抗乳腺癌候選藥物的優(yōu)化模型，經(jīng)過分析，建議選用包含20個分子描述符的化合物作為抗乳腺癌候選藥物．

5 結(jié)語

本文針對華為杯中國研究生數(shù)學(xué)建模競賽D題中的抗乳腺癌候選藥物等問題，通過綜合運用灰色關(guān)聯(lián)度分析、BP神經(jīng)網(wǎng)絡(luò)、決策樹與回歸模型等方法和理論建立模型，此模型在藥物生物活性預(yù)測及ADMET分類預(yù)測上均取得良好性能，較好地解決了候選藥物的優(yōu)化問題．在構(gòu)建模型的過程中，通過對數(shù)據(jù)擴增以及預(yù)測算法模型的迭代優(yōu)化等方式進一步增強了預(yù)測工具的通用性．可將該模型進行推廣，應(yīng)用于虛擬藥物篩選流程，為計算機輔助藥物設(shè)計與藥物發(fā)現(xiàn)提供新思路，具有較好的借鑒意義．

[1] 中國學(xué)位與研究生教育學(xué)會．華為杯中國研究生數(shù)學(xué)建模競賽D題[EB/OL]．（2021-09-01）[2021-10-16]．https://cp

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

[2] 路珩，張一奇．雄激素受體在雌激素受體陽性乳腺癌患者中的表達(dá)及其臨床意義[J]．中國現(xiàn)代醫(yī)學(xué)雜志，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] 張翠鋒，謝海棠，潘國宇．大分子藥物的吸收、分布、代謝、排泄和毒性特征及藥代模型的應(yīng)用[J]．藥學(xué)學(xué)報，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] 劉雅紅，賀利民，梁智斌，等．用于預(yù)測化合物活性的兩級擬合QSAR模型的構(gòu)建方法：中國：102930113B[P]．（2015-02-03）[2021-10-16]．https://wenku.baidu.com/view/4c597ce9ce1755270722192e453610661fd95ac3?fr=xueshu_top．

[7] 鄧聚龍．灰理論基礎(chǔ)[M]．武漢：華中科技大學(xué)出版社，2002．

[8] 虞曉芬，傅玳．多指標(biāo)綜合評價方法綜述[J]．統(tǒng)計與決策，2004（11）：119-121．

[9] 羅黨，劉思峰．灰色關(guān)聯(lián)決策方法研究[J]．中國管理科學(xué)，2005（1）：102-107．

[10] 韓中庚．?dāng)?shù)學(xué)建模方法及其應(yīng)用[M]．北京：高等教育出版社，2005．

[11] 謝良旭，李峰，謝建平，等．基于融合神經(jīng)網(wǎng)絡(luò)模型的藥物分子性質(zhì)預(yù)測[J]．計算機科學(xué)，2021，48（9）：251-256．

[12] 潘斌．?dāng)?shù)學(xué)建模教程[M]．北京：化學(xué)工業(yè)出版社，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的統(tǒng)計教學(xué)課例賞析：兩個變量的線性相關(guān)（第2課時）[J]．?dāng)?shù)學(xué)通報，2016，55（12）：20-23．

[17] 黃忠裕．初等數(shù)學(xué)模型[M]．北京：科學(xué)出版社，2013．

[18] 李琳，陳德釗，束志恒，等．基于預(yù)處理的決策樹在化學(xué)數(shù)據(jù)挖掘中的應(yīng)用[J]．分析化學(xué)，2005（8）：1091-1094．

[19] 郭曉龍，蔣艷，邱路．決策樹分類模型預(yù)測蛋白質(zhì)相互作用的應(yīng)用研究[J]．生物醫(yī)學(xué)工程學(xué)雜志，2013，30（5）：952-956．

[20] 顧耀文，張博文，鄭思，等．基于圖注意力網(wǎng)絡(luò)的藥物ADMET分類預(yù)測模型構(gòu)建方法[J]．?dāng)?shù)據(jù)分析與知識發(fā)現(xiàn)，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

安徽省自然科學(xué)研究項目（1908085MF186）；安徽省高校自然科學(xué)研究重點項目（KJ2019A0589）；安徽省質(zhì)量工程項目

（2020jyxm1670，2020jxtd）

董奕鑫（1998-），女，山東臨沂人，在讀碩士研究生，從事數(shù)學(xué)教學(xué)評價、數(shù)學(xué)建模研究．E-mail：dongyixin1998@163.com