抗乳腺癌候選藥物的優化模型

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

抗乳腺癌候選藥物的優化模型

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

（淮北師范大學 數學科學學院，安徽 淮北 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