1. 引言
近年来乳腺癌已成为致死率较高的癌症之一,乳腺癌是乳腺上皮细胞在多种致癌因子的作用下,发生增殖失控的现象。乳腺癌的发展与雌激素受体密切相关,有研究发现,雌激素受体α亚型在不超过10%的正常乳腺上皮细胞中表达,但大约在50%~80%的乳腺肿瘤细胞中表达;而对ERα基因缺失小鼠的实验结果表明,ERα确实在乳腺发育过程中扮演了十分重要的角色。因此,ERα被认为是治疗乳腺癌的重要靶标,能够拮抗ERα活性的化合物可能是治疗乳腺癌的候选药物。目前,在药物研发中,为了节约时间和成本,通常采用建立化合物活性预测模型的方法来筛选潜在活性化合物。根据提供的ERα拮抗剂信息(1974个化合物样本,每个样本都有729个分子描述符变量,1个生物活性数据,5个ADMET性质数据),构建化合物生物活性的定量预测模型和ADMET性质的分类预测模型,从而为同时优化ERα拮抗剂的生物活性和ADMET性质提供预测服务。建立数学模型研究下列问题:
1、针对1974个化合物的729个分子描述符进行变量选择,分析变量对生物活性影响的重要性进行排序,给出前20个对生物活性最具有显著影响的分子描述符(即变量),然后说明分子描述符筛选过程及其合理性。
2、结合问题一,选用不超过20个分子描述符变量,构建化合物对ERα生物活性的定量预测模型,对50个化合物进行IC50值和对应的pIC50值预测,并将结果分别填入对应的表。
3、利用ADMET数据,分别构建化合物的Caco-2、CYP3A4、hERG、HOB、MN的分类预测模型。通过5个分类预测模型,对50个化合物进行相应的预测,并将结果填入对应的表。
4、寻找并阐述化合物的分子描述符,以及这些分子描述符在什么取值或者处于什么取值范围时,能够使化合物对抑制ERα具有更好的生物活性,同时具有更好的ADMET性质,在给定的五个性质中,至少选择三个性质较好。
2. 问题分析
研究发现,乳腺癌的发展与雌激素受体密切相关,而对ERα基因缺失小鼠的实验结果表明,ERα在乳腺发育过程中扮演了十分重要的角色;因此能够拮抗ERα活性的化合物可能是治疗乳腺癌的候选药物。所以通过预测化合物的生物活性来筛选药物是高效且快速的方法,为了解决这个问题,我们进行了具体分析。
2.1. 问题一的分析
问题一要求我们建立模型针对分子描述符进行变量筛选,将变量对生物活性影响的重要性进行排序,并按照降序排列筛选出前20个对生物活性最具有影响力的分子描述符。
2.2. 问题二的分析
问题二要求我们选择不超过20个分子描述符变量,构建化合物对ERα生物活性的定量预测模型,然后使用构建的预测模型,对文件“ERα_activity.xlsx”的test表中的50个化合物进行IC50值和对应的pIC50值进行预测。首先我们可以选用第一题所得到的前20个对生物活性最具有显著影响的分子描述符表示为自变量,生物活性表示为因变量;其次构建多元非线性回归方程解出各因变量的系数,最后得到预测模型。
2.3. 问题三的分析
问题三要求分别构建化合物的Caco-2、CYP3A4、hERG、HOB、MN的分类预测模型,分类是指定性输出预测,传统的多元回归分析预测值不能落入0-1区间,因此不能针对解决二分变量(0,1编码)的预测问题,因此,本题将应用全连接的单层神经网络分别构建化合物的Caco-2、CYP3A4、hERG、HOB、MN的分类预测模型。首先我们利用题一的模型找出分别影响Caco-2、CYP3A4、hERG、HOB、MN的10个特征,然后通过找出来的10个特征分别使用全连接的单层神经网络构建分类预测模型。
2.4. 问题四的分析
问题四要求寻找化合物的某些分子描述符并得到其取值范围,使得化合物对抑制ERα具有更好的生物活性,同时具有更好的ADMET性质,基于问题一和问题三筛选出的分子描述符,我们采用博弈论的方法将问题四看成是卖方与买方的对策问题,进而用博弈论中的零和二人混合对策来进行数学描述,在极大化化合物对抑制ERα具有更好的生物活性的情况下,同时具有更好的ADMET性质。
3. 模型假设
1、所有化合物的生物活性值、分子描述符、ADMET性质的数据测量正确;
2、没有其他的分子描述符会对化合物的生物活性值和ADMET性质有影响;
3、建立模型时假设不考虑人为因素对数据的影响;
4、在优化主要变量时认为所提出的预测模型结果准确。
4. 符号说明
由于文章符号较多,为了方便阅读,将所有符号整理如下
5. 问题一的模型建立与求解
5.1. 模型的建立
问题1模型建立过程如图1。
5.2. 模型的求解
5.2.1. 数据预处理
数据预处理是从数据中检测,纠正或者损坏不准确,不适用模型的记录的过程,为了提高要进行观察和分析的数据的质量,我们在进行建模前去掉全为0的分子描述符的列 [1] 。
利用Python实现随机森林算法,综合考虑到算法速度和算法准确率,设定K = 500,M = 100。运行程序得到504个分子描述符的贡献度排名,将贡献度由大到小排序。考虑到下一步的高相关性滤波操作会对进一步变量进行降维,故在这一步中先筛选出排名处于前60的分子描述符,如下表1所示。可以发现排名为60的分子描述符对pIC50值的归一化贡献度只有0.284%,且前60个分子描述符贡献度有83.78%,说明选取贡献度排名前60名的分子描述符已经可以有足够的信息来预测pIC50值。

Table 1. Name and contribution size of the top 60 molecular descriptors with contribution ranking
表1. 贡献度排名前60的分子描述符名称及贡献度大小
5.2.2. 高相关性变量去耦合
本题中各变量间为非线性关系,选择距离相关系数作为衡量变量间相关性的指标 [1] 。对于两个随机向量,
,
,记
为观察到的随机样本,x,y间的距离相关系数(dCor)可以定义为:
(5.1)
其中:
(5.2)
(5.3)
(5.4)
(5.5)
同理可计算:
(5.6)
(5.7)
根据题目要求,去耦合后的分子描述符之间应该具有中相关以下的相关性,故将距离相关系数(dCor)的阈值设为0.6。基于此,对表1中的分子描述符进行相关性检验。首先计算两两分子描述符间的距离相关系数,然后判断两者的相关系数是否大于所设阈值。若大于,则说明两分子描述符相关性高,对两分子描述符的贡献值排名靠后的那一分子描述符采取删除操作。若小于,则说明两分子描述符间并不强相关,则将两分子描述符都暂时予以保留。之后重复上述操作,直至遍历结束。使用python编写去耦合程序,程序流程图如图2所示。

Figure 2. Decoupling program flowchart
图2. 去耦合程序流程图
根据此去耦合程序,对表1中分子描述符进行高相关性滤波处理,前60个贡献值高的分子描述符经滤波处理后剩余20个分子描述符,如下表2所示。降维后的这20个分子描述符是对生物活性最具有显著影响的分子描述符,供下文建立模型时使用 [2] 。

Table 2. Residual molecular descriptors after filtering processing
表2. 滤波处理后剩余分子描述符
5.3. 合理性评价
从分子描述符降维过程中采用的算法及处理流程来看:随机森林得到的分子描述符贡献度排名,保留排名靠前的分子描述符从而保证了所选取变量与因变量之间的相关性,而设计的基于距离相关性高相关度变量滤波算法则保证了降维后分子描述符之间的独立性,可根据所提取的分子描述符之间的距离相关系数计算结果,对分子描述符之间的相关程度进行可视化,如图3所示。可见选取的前20个分子描述符之间相关性低,独立性较好。

Figure 3. Correlation distribution map of the top 20 main molecular descriptors
图3. 前20个主要分子描述符的相关性分布图
6. 问题二的模型建立与求解
6.1. 数学模型的建立
问题二的流程图如图4所示,生物活性通常受到多个因素的影响,根据问题1我们选择前20个对生物活性最具有显著影响的分子描述符作为变量,建立用于预测生物活性的多元线性回归模型,见式
(6.1)
式中y为生物活性值;x为20个生物活性影响分子描述符;
为待求系数;
为随机误差,
(6.2)
将式(6.1)写成误差方程的形式,如下
(6.3)
式中,
为改正数,由于pIC50由多变量控制,因此,可将式(6.2)写为矩阵形式,见式(6.4):
(6.4)
式中,
,
,
,
由最小二乘准则可得,待定系数解算公式为,
,将
带入式(6.1)即可求得预测pIC50值。
由于影响因子之间单位不一致,容易出现量纲不同的情况,因此,本文在解算时,使用参数中心化,以去除量纲的影响。
将上述的待定系数解算公式进行参数中心化变更,见式(6.5):
(6.5)
式中,
,
但在使用MATLAB软件求解过程中,我们发现经过7次对异常点进行处理后的时序残差图仍存在较多的异常点,如下图5所示,所以下面我们考虑构建多元非线性回归模型。

(a)(b)
Figure 5. Time series residual comparison chart, (a) Residual diagram of the first running time sequence; (b) Residual diagram of the seventh running time sequence
图5. 时序残差对比图,(a) 第一次运行时序残差图;(b) 第七次运行时序残差图
与多元线性回归相对应的是多元非线性回归,然而多元非线性回归能够对因变量的影响产生更加重要的意义。如本题中预测pIC50中,共有20个变量对因变量产生影响,应用多元线性回归的实用意义很大,在实际应用中变量之间的非线性变化,以及交互项对因变量的的预测产生重要的影响,非线性变化的加入也能够让因变量得到更加准确的效果。多元非线性回归的公式如下所示:
(6.6)
这里的指非线性函数,包括幂次、指数、曲线、对数函数等。这里还应该注意的是,由于不同的指标具有不同的单位,不同的单位加入到公式运算中,虽然得到了如上公式中的系数,但这时的系数值大小不能代表该指标的重要性权重。所以本例中应用了Python软件进行标准化的操作,得到了标准化指标。
基于第二问主要变量的确定,记pIC50为y。以下将构建pIC50与20个主要变量的回归模型。自变量符号如下表3所示:
首先我们选用第一题所得到的前20个对生物活性最具有显著影响的分子描述符作为20个变量,其次我们使用python构建多元非线性回归得到预测模型。最后我们将结果分别填入相应的test表中。
6.2. 数学模型的求解
我们用python软件的相关性分析功能,对各自变量之间进行相关性检验,发现变量之间没有很高的相关性,所以我们对一些数值较大的变量进行取自然对数处理,以此来降低其数值大小。最后我们将train表中的数据作为模型训练样本,test表中的数据作为预测样本并与实际值进行比较。
根据python软件运行结果,得到新的多元非线性回归方程模型:
(6.7)
多元非线性回归预测值与实际值对比图如下图6所示,我们可以得到pIC50预测值与实际值之间的差值在0.2以内,有着较高的准确度,因此可以用来进行pIC50值的预测。

Figure 6. Comparison between predicted and actual values of multiple nonlinear regression
图6. 多元非线性回归预测值与实际值对比图
7. 问题三的模型建立与求解
7.1. 数学模型的建立
问题三的求解过程如图7,通过利用题一的模型我们可以分别得到前10个对化合物的Caco-2、CYP3A4、hERG、HOB、MN最具有显著影响的分子描述符,如下表4所示,供下文分别建立Caco-2、CYP3A4、hERG、HOB、MN损失预测模型时使用。
通过对分子描述符之间的相关程度进行可视化,如图8所示。可见分别选取的前10个分子描述符之间相关性低,独立性较好。


(a) [Caco-2](b) [CYP3A4]

(c) [hERG] (d) [HOB]
(e) [MN]
Figure 8. Correlation distribution map of the top 10 main molecular descriptors
图8. 前10个主要分子描述符的相关性分布图
从数据类型上来看,分子描述符之间可能具有高度非线性关系。全连接单层神经网络是人工智能网络的一个典型算法,其本身具有很强的非线性映射能力,非常适合解决一些非线性问题 [3] 。另外其网络拓扑结构简单,而且具有较高的计算精度。因此在此问中,我们将样本数据集数据进行标准化处理之后,采用train数据用于全连接单层神经网络模型训练,使用test数据对模型精度进行验证以评价模型合理性,本题我们选择只有一层的全连接神经网络,即二分类logistic回归模型。
对于维度为m + 1特征为x样本的二分类问题,有负类记为0,正类记为1,即对于类别y,有
我们期望找到一个
,使得
(7.1)
其中,θ为待优化的参数,使得在对未知类别的样本x0分类时,
为样本为正类的概率。即分类准则如下:
(7.2)
在线性回归中,我们常找一组参数
计算
(7.3)
设置阈值T,通过
与T的大小关系判断正负类。
而在Logistic回归中,我们引入Sigmoid函数
(7.4)
Logistic回归取hypothesis function为
(7.5)
即
等价于正类的概率,由Sigmoid函数图像可知,当
时,判定为正类,当
时,判定为负类。
与线性回归问题类似,Logistic同样需要定义代价函数使用梯度下降法优化参数
由于Sigmoid函数的使用,若使用与线性回归相同的二次损失函数,优化问题将变为非凸问题,即可能存在很多局部最优解。Logistic回归采用以下损失函数:
(7.6)
对于样本数目为n的训练集,定义目标函数为:
(7.7)
优化目标为:找到令
最小的θ。
7.2. 模型求解和分析
各个化合物的分类预测模型的交叉熵损失图如下图9所示,从图中可以观察到预测值与实际值的趋势大体一致,误差较小,所以模型较为合理,可以用来预测ADMET性质的值。

(a) [Caco-2] (b) [CYP3A4]
(c) [hERG](d) [HOB]
(e) [MN]
Figure 9. Cross entropy loss graph of top 10 molecular descriptors of five ADMET properties
图9. 五种ADMET性质的前10个分子描述符的交叉熵损失图
8. 模型评价
8.1. 模型的优点
1) 充分考虑了各分子描述符与生物活性值之间的非线性关系,使用了随机森林回归、距离相关系数等适用于处理非线性特征的方法。所获得的主要变量意义明确,符合实际。
2) 针对主要变量与ADMET性质之间复杂的关系,选用了全连接的单层神经网络机器学习算法,同时数据集与测试集分开,所得到的结果最大相对误差很小,所建立的预测模型精度、鲁棒性表现优秀。
3) 使用博弈论模型,符合分子描述符实际筛选情况,能更好的找到既具有较高生物活性和较好ADMET性质的分子描述符。
4) 算法速度快,响应性好。
8.2. 模型的缺点
1) 模型过于理想化,应该考虑细分因素深度探究,升华模型。
2) 本文构建的神经网络预测模型训练样本数较少,日后可获得更多数据对其进行修正。
3) 多元非线性回归模型有更好的函数形式,由于变量过多,我们没有进行多次尝试得到更优的一个函数。