基于Matlab的抗乳腺癌候选药物优化建模——2021年中国研究生数学建模竞赛D题
一、背景介绍
乳腺癌是目前世界上最常见，致死率较高的癌症之一。乳腺癌的发展与雌激素受体密切相关，有研究发现，雌激素受体α亚型（Estrogen receptors alpha, ERα）在不超过10%的正常乳腺上皮细胞中表达，但大约在50%-80%的乳腺肿瘤细胞中表达；而对ERα基因缺失小鼠的实验结果表明，ERα确实在乳腺发育过程中扮演了十分重要的角色。目前，抗激素治疗常用于ERα表达的乳腺癌患者，其通过调节雌激素受体活性来控制体内雌激素水平。因此，ERα被认为是治疗乳腺癌的重要靶标，能够拮抗ERα活性的化合物可能是治疗乳腺癌的候选药物。比如，临床治疗乳腺癌的经典药物他莫昔芬和雷诺昔芬就是ERα拮抗剂。
目前，在药物研发中，为了节约时间和成本，通常采用建立化合物活性预测模型的方法来筛选潜在活性化合物。具体做法是：针对与疾病相关的某个靶标（此处为ERα），收集一系列作用于该靶标的化合物及其生物活性数据，然后以一系列分子结构描述符作为自变量，化合物的生物活性值作为因变量，构建化合物的定量结构-活性关系（Quantitative Structure-Activity Relationship, QSAR）模型，然后使用该模型预测具有更好生物活性的新化合物分子，或者指导已有活性化合物的结构优化。
一个化合物想要成为候选药物，除了需要具备良好的生物活性（此处指抗乳腺癌活性）外，还需要在人体内具备良好的药代动力学性质和安全性，合称为ADMET（Absorption吸收、Distribution分布、Metabolism代谢、Excretion排泄、Toxicity毒性）性质。其中，ADME主要指化合物的药代动力学性质，描述了化合物在生物体内的浓度随时间变化的规律，T主要指化合物可能在人体内产生的毒副作用。一个化合物的活性再好，如果其ADMET性质不佳，比如很难被人体吸收，或者体内代谢速度太快，或者具有某种毒性，那么其仍然难以成为药物，因而还需要进行ADMET性质优化。为了方便建模，本试题仅考虑化合物的5种ADMET性质，分别是：1）小肠上皮细胞渗透性（Caco-2），可度量化合物被人体吸收的能力；2）细胞色素P450酶（Cytochrome P450, CYP）3A4亚型（CYP3A4），这是人体内的主要代谢酶，可度量化合物的代谢稳定性；3）化合物心脏安全性评价（human Ether-a-go-go Related Gene, hERG），可度量化合物的心脏 毒性；4）人体口服生物利用度（Human Oral Bioavailability, HOB），可度量药物进入人体后被吸收进入人体血液循环的药量比例；5）微核试验（Micronucleus，MN），是检测化合物是否具有遗传毒性的一种方法。
二、数据集介绍及建模目标
本试题针对乳腺癌治疗靶标ERα，首先提供了1974个化合物对ERα的生物活性数据。这些数据包含在文件“ERα_activity.xlsx”的training表（训练集）中。training表包含3列，第一列提供了1974个化合物的结构式，用一维线性表达式SMILES（Simplified Molecular Input Line Entry System）表示；第二列是化合物对ERα的生物活性值（用IC50表示，为实验测定值，单位是nM，值越小代表生物活性越大，对抑制ERα活性越有效）；第三列是将第二列IC50值转化而得的pIC50（即IC50值的负对数，该值通常与生物活性具有正相关性，即pIC50值越大表明生物活性越高；实际QSAR建模中，一般采用pIC50来表示生物活性值）。该文件另有一个test表（测试集），里面提供有50个化合物的SMILES式。
其次，在文件“Molecular_Descriptor.xlsx”的training表（训练集）中，给出了上述1974个化合物的729个分子描述符信息（即自变量）。其中第一列也是化合物的SMILES式（编号顺序与上表一样），其后共有729列，每列代表化合物的一个分子描述符（即一个自变量）。化合物的分子描述符是一系列用于描述化合物的结构和性质特征的参数，包括物理化学性质（如分子量，LogP等），拓扑结构特征（如氢键供体数量，氢键受体数量等），等等。关于每个分子描述符的具体含义，请参见文件“分子描述符含义解释.xlsx”。同样地，该文件也有一个test表，里面给出了上述50个测试集化合物的729个分子描述符。
最后，在关注化合物生物活性的同时，还需要考虑其ADMET性质。因此，在文件“ADMET.xlsx”的training表（训练集）中，提供了上述1974个化合物的5种ADMET性质的数据。其中第一列也是表示化合物结构的SMILES式（编号顺序与前面一样），其后5列分别对应每个化合物的ADMET性质，采用二分类法提供相应的取值。Caco-2：‘1’代表该化合物的小肠上皮细胞渗透性较好，‘0’代表该化合物的小肠上皮细胞渗透性较差；CYP3A4：‘1’代表该化合物能够被CYP3A4代谢，‘0’代表该化合物不能被CYP3A4代谢；hERG：‘1’代表该化合物具有心脏 毒性，‘0’代表该化合物不具有心脏 毒性；HOB：‘1’代表该化合物的口服生物利用度较好，‘0’代表该化合物的口服生物利用度较差；MN：‘1’代表该化合物具有遗传毒性，‘0’代表该化合物不具有遗传毒性。同样地，该文件也有一个test表，里面提供有上述50个化合物的SMILES式（编号顺序同上）。
建模目标：根据提供的ERα拮抗剂信息（1974个化合物样本，每个样本都有729个分子描述符变量，1个生物活性数据，5个ADMET性质数据），构建化合物生物活性的定量预测模型和ADMET性质的分类预测模型，从而为同时优化ERα拮抗剂的生物活性和ADMET性质提供预测服务。
三、需解决问题
问题1. 根据文件“Molecular_Descriptor.xlsx”和“ERα_activity.xlsx”提供的数据，针对1974个化合物的729个分子描述符进行变量选择，根据变量对生物活性影响的重要性进行排序，并给出前20个对生物活性最具有显著影响的分子描述符（即变量），并请详细说明分子描述符筛选过程及其合理性。
方法：熵权法
问题2. 请结合问题1，选择不超过20个分子描述符变量，构建化合物对ERα生物活性的定量预测模型，请叙述建模过程。然后使用构建的预测模型，对文件“ERα_activity.xlsx”的test表中的50个化合物进行IC50值和对应的pIC50值预测，并将结果分别填入“ERα_activity.xlsx”的test表中的IC50_nM列及对应的pIC50列。
方法：svm神经网络
2.1.1支持向量机理论
支持向量机（Support Vector Machine, SVM），最早由Vapnik于1995年提出，是基于统计学理论的人工学习方法，广泛应用于数据的分类，预测和回归。
其根本思想如图2-1所示，是用一个映射过程将高度非线性关系样本映射到一个n维空间，然后通过风险结构优化原理去处理分类问题。空间中的回归函数是f(x, ω),可以通过下列方程表示：
f(x,ω)=∑_(j=1)^n▒〖ω_j*〗 g_j (x)+b
其中g_j (x)是为映射函数，ω_j是权系数，b是阈值。
根据结构风险最小化原理，回归优化约束可以表示为:
min 1/2 ‖ω‖2+c∑_(i=1)n▒〖(ε_i+ε_j)〗
sunject to {█(y_i-f(x_i,ω)≤ε+ε_i*@f(x_i,ω)-y_i≤ε+ε_i*@ε_i,ε_i^≫0,i=1………n)┤
其中ε_i^ ，ε为松弛变量，通过求解上式的约束优化问题得到α_i*，确定好核函数之后通过α_i计算出阈值b，最终得到其最终的决策函数为：
f(x)=∑_(i=1)(n_sv)▒〖(α_i-α_i)k(x,x_i)〗
subject to 0≤α_i^*≤c,0≤α_i≤c
k(x,x_i)为核函数，c为惩罚参数。
图2-1 支持向量机原理
2.2预测模型的建立
2.2.1数据的获取及预处理
2.2.3支持向量机模型的建立
I．训练集的确定：
T={(x_1,y_1 ),⋯,(x_l,y_l )}∈〖（X×Y）〗^l
式中：{x_i,y_i }(i=1) ,x_i∈Rn为第i个组的影响因素，y_i∈Rn为第i个具体期望 输出值，来自于训练集中的数据。
II．选择适当的核函数，构造并求解最优化问题，本文选择核函数为径向基函数：
K(x,x_i )=exp(-gamma‖x-x_i ‖^2 )
根据结构最优化原则，回归优化的目标表示为：
〖min〗^α 1/2 ∑(i=1)j▒∑_(j=1)l▒y_i y_j α_i α_j exp(-gamma‖x_i-x_j ‖^2 )-∑_(j=1)^l▒α_j
subject to ∑_(i=1)^l▒y_i α_i=0，0≤α_i≤c，i=1,2,…,l
经过计算得到最优解：α*=〖（α_1,…,α_l*）〗T，c为惩罚参数，g(-gamma)核函数参数。
IV．计算阈值：
b*=y_j-∑_(i=1)l▒〖y_i α_i^ 〗 〖K（x_i-x_j）exp〗(-gamma‖x_i-x_j ‖^2 )
由III中的最优解a代入（10），计算得到的阈值为b=-0.5281。
V．构造决策函数
将计算得到的a*，b*，g，选择的核函数代入决策函数得到：
f(x)=sgn(∑_(i=1)^l▒ω_i exp(g‖x_i-x‖^2 )+b^*)
其中x是待预测样本，决策函数用来对样本数据检验，并测试其泛化能力。
问题3. 请利用文件“Molecular_Descriptor.xlsx”提供的729个分子描述符，针对文件“ADMET.xlsx”中提供的1974个化合物的ADMET数据，分别构建化合物的Caco-2、CYP3A4、hERG、HOB、MN的分类预测模型，并简要叙述建模过程。然后使用所构建的5个分类预测模型，对文件“ADMET.xlsx”的test表中的50个化合物进行相应的预测，并将结果填入“ADMET.xlsx”的test表中对应的Caco-2、CYP3A4、hERG、HOB、MN列。
方法：bp神经网络
如图2-2所示，输入为P为R×1维输入矢量，网络层由权值矩阵W，阀值矢量b，求和单元和传递函数运算单元f组成，S个神经元的输出构成了s×1维的神经网络输出矢量a.
a=f(W×P+b)
其中，输入层网络权值矩阵W和阀值矢量b的具体形式如下：
对于一个三层网络结构来说，将上图集训推广得到三层网络结构的表达式为：
a=f(L×W_1×P×f(L×W_1×P+b_1 )+b_2)
学习特征是神经网络的基本特征，其可通过调整网络权值和阈值来实现神经网络的学习和训练。根据学习过程的组织和管理方式，学习方法分为两类：监督学习和无监督学习。
如图2-3所示，BP神经网络，即误差反向传播算法，是最经典的监督学习算法之一，网络往往是一定数量的训练样本进行训练。训练样本通常由输入向量和目标向量组成。在学习过程中对神经网络的实际输出与目标输出进行比较，并根据比较结果或误差，根据一定的规则或算法对网络权值和阈值进行调整，使网络的输出逐渐接近目标值。
2.2.3 BP神经网络的建立
图2-5建立BP神经网路在训练过程的进程图，建模过程中使用的模型参数如表2所示。采用的步幅为10，动量因子0.01，训练次数10000，目标误差0.01，其中输入的层节点数为20，隐藏层节点数为7，传递函数分别为：logsig和tansig函数。
表 2 建模参数
参数 show lr mc epochs goal
数值 10 0.01 0.9 10000 0.01
训练集和测试集
问题4. 寻找并阐述化合物的哪些分子描述符，以及这些分子描述符在什么取值或者处于什么取值范围时，能够使化合物对抑制ERα具有更好的生物活性，同时具有更好的ADMET性质（给定的五个ADMET性质中，至少三个性质较好）。
方法：算法寻优(BP-MIV)
附件：
附件一：ERα_activity.xlsx
附件二：Molecular_Descriptor.xlsx
附件三：分子描述符含义解释.xlsx
附件四：ADMET.xlsx
最后，有相关需求，欢迎联系我们。
