用机器学习预测药物在靶点上的停留时间

传统上，药物-靶点复合物形成和解离的速度，不被认为是影响药物在体内作用或持续时间的主要因素。2006 年引入了药物-靶点停留时间的概念后，这种传统的观点受到了挑战。

本文为《Artificial Intelligence in Drug Design》一书第8章的内容，介绍了影响停留时间的因素，以及两种用机器学习预测药物-靶点停留时间的模型。

摘要

药物在靶点上的停留时间（即与特定蛋白靶点结合的持续时间），在某些蛋白家族中，对药效的影响比结合亲和力更重要。为了在药物发现中对停留时间进行有效的优化，需要开发能够预测该指标的机器学习模型。预测停留时间的主要挑战之一是数据的匮乏。本章概述了目前所有可用的配体动力学数据，提供了一个迄今为止最大的GPCR-配体动力学数据资料库（公开来源）。为了帮助解读动力学数据的特征（这些特征的纳入对于预测停留时间的计算模型是有益的），我们总结了影响停留时间的实验证据。最后，概述了用机器学习预测停留时间的两种不同的工作流程。第一个是根据配体特征训练的单靶点模型；第二个是根据分子动力学模拟产生的特征来训练的多靶点模型。

1 简介

药物在靶点上停留时间，即配体解离率的倒数，对于某些靶点来说，对药效的影响比平衡结合亲和力更大。已经有几个将停留时间与体内药效关联起来的研究发表。对作用于12个不同靶点的50种药物的分析显示，70%的停留时间长的药物比停留时间短的同类药物具有更高的药效。停留时间与药效相关的证据在GPCRs中最为丰富。

研究发现，毒蕈性乙酰胆碱M3受体激动剂的药效仅与其停留时间有关，而与它们的结合亲和力无关。同样，对于A2A腺苷受体激动剂，停留时间是唯一被发现与体内药效相关的因素。停留时间和药效之间的相关性也被观察到与拮抗剂有关。抗组胺拮抗剂的停留时间被发现与它们抑制细胞中H1组胺受体的能力相关。

因此，在许多不同的GPCR例子中，激动剂的功效和拮抗剂的抑制强度已被证明与停留时间更相关（相比于结合亲和力）。应该注意的是，延长停留时间不仅影响药效，还可能对药物剂量间隔产生影响。噻托品与伊普拉托品相比，停留时间长50倍，两者都是M3毒蕈碱受体的配体，这意味着噻托品的用药频率可以降低。

配体在靶点蛋白上的停留时间与在非靶点蛋白上的停留时间的差异，决定了产生非靶点副作用的概率。传统上，靶点选择性是以非靶点蛋白与靶点蛋白的结合亲和力值的比率来衡量的（平衡选择性）。正如Copeland等人所指出的，药物在血浆中的浓度不是恒定的，因此，药物对不同蛋白质的解离率决定了药物的时间选择性。换句话说，在给药时，最初的选择性是由靶点蛋白和非靶点蛋白之间的结合亲和力的差异决定的，但是，随着血浆中药物浓度的降低，在靶点蛋白和非靶点蛋白上的停留时间的差异也决定了药物的整体选择性。

人们对在药物发现中优化药物-靶点的停留时间表示关注。数学模型表明，只有当停留时间超过清除时间时，才会对药物的占有率负责，而对于市场上的许多小分子药物，情况并非如此。还有一个问题是蛋白质的周转率；如果一个靶点有很高的周转率，可能就没有必要延长停留时间。例如，一种停留时间很长的酪氨酸激酶抑制剂，停留时间为一周，但受体周转率很高，24小时后体内占用率不到50%。延长停留时间应在蛋白质周转率和配体清除时间的背景下考虑。然而，延长靶点占有率可能不是长停留时间导致更好药效的唯一机制。持续的拮抗剂结合与长的药物-靶点停留时间可能会阻止瞬时激动剂的结合，因为靶点已经被占据。

由于这些不同的发现，有人提出，在药物发现的苗头化合物-先导化合物和先导化合物优化阶段，考虑药物-靶点停留时间是至关重要的。然而，为了有效地做到这一点，我们需要计算工具来预测和合理化停留时间，这是一个比结合亲和力更难通过实验来确定的属性。机器学习（ML）多年来一直处于药物发现的前沿，部分原因是它可以对人类无法理解的复杂数据进行回归，在寻找计算预测药物-靶点停留时间的方法时，它是一个很好的候选者。为了做出准确的ML回归模型，我们需要选择正确的输入特征，然后将其输入ML模型，给出预测值输出。特征的选择对于ML回归模型能否成功做出准确的预测，并与实验结果趋于一致至关重要。

1.1药物-靶点停留时间长短的特征

为了延长停留时间而不影响平衡结合亲和力，需要实现过渡结合状态的稳定性。这些过渡状态的细节目前很难通过实验来评估，但可以通过计算模拟（分子动力学）来观察。随着最近X射线自由电子激光器的发展，用时间分辨晶体学鉴定配体结合过渡状态在实验上变得更加可行。然而，在X射线自由电子激光器的可用性变得更加广泛之前，需要配体在结合状态下的信息（计算模拟和实验确定），来帮助揭示药物-靶点停留时间的分子决定因素。

在计算研究和实验研究中都有关于赋予药物-靶点长停留时间的重要特征的建议。对辉瑞公司2000个有停留时间值的化合物数据库的调查显示，延长的停留时间和配体的大小之间有关联。如图1所示，对所有可用的GPCR-配体动力学数据（500个化合物）的分析支持了这一点，它揭示了分子大小和停留时间之间的弱的正相关。在某些情况下，发现分子量只与停留时间相关，而与结合亲和力无关。在其他情况下，停留时间和配体分子量之间有很强的相关性。在这些情况下，应注意确保用于确定停留时间的方法优于这种简单的线性相关，而不是仅仅将配体大小作为确定停留时间的指标。

图1配体大小和药物-靶点停留时间的相关性。

散点图显示了500个GPCR配体的药物-靶点停留时间与配体分子量（MW）的相关性（上），以及停留时间与配体中的环数（no_rings）的相关性（下）。最小二乘法线性回归线显示为橙色虚线，这些相关性的强度显示为R2

除了配体的大小，与水的相互作用也是已知的决定停留时间的重要因素。被"掩埋"的亲水相互作用，即被水屏蔽的相互作用，已经在计算和实验中被证明可以延长停留时间。这些相互作用具有更高的能量屏障，这意味着它们更稳定，更少瞬时性。水也被证明是赋予GPCR受体（A2A受体）长停留时间的一个重要因素，在结合状态下溶解度可能降低的化合物，停留时间增加。在我们以前发表的使用转向分子动力学的文章中，观察到结合和未结合的配体之间的溶解度变化与停留时间密切相关，部分原因是溶解度能量变化较大的配体与蛋白质的亲水相互作用更可能被"掩埋"。

有人指出，结合位点的灵活性也会影响停留时间。配体可以通过稳定蛋白质内的相互作用来影响结合部位的稳定性，例如ZM-241,385可以稳定A2A受体中的蛋白质内盐桥，增加停留时间。通过定点诱变破坏该盐桥，可以使该配体的停留时间减少8倍，从84分钟减少到5分钟。前面提到的毒蕈碱亚家族拮抗剂噻托品，对M3的停留时间比M2受体亚型长10倍。分子动力学模拟显示，M2受体的第二个胞外环更加灵活。配体与受体的相互作用降低了ECL2的蛋白质灵活性，已被证明会增加停留时间。

匹配分子对（Matched molecular pair, MMP）分析最近被用来尝试理解结构-动力学关系。MMP分析一般在解释配体结合率方面比较成功，配体极性的增加导致结合率的降低。只有少数的MMP转化被确定为停留时间明显延长，而结合亲和力和结合速率基本上没有受到干扰。其中一个转化是去除噻托品的羟基（去羟基噻托品），它使M3毒蕈碱受体的停留时间减少了56倍。噻托品中的这个羟基已知会与Asn507形成氢键。由于周围的芳香族残基，这种相互作用可能是一种被"掩埋"的亲水作用，因此，羟基的存在延长了停留时间。

1.2 以前用于预测停留时间的ML方法

有许多已发表的方法试图用ML来预测药物-靶点的停留时间，但在数量上远远少于基于ML预测结合亲和力值的方法。这可能是由于两个原因：首先，与结合亲和力相比，严重缺乏停留时间的训练数据；其次，相对而言，最近才发现停留时间在药物发现中的作用。训练数据的缺乏可以体现在目前大多数的停留时间预测方法都是在少量（100个以下）的化合物上训练的。此外，这些方法只用于两个蛋白质靶点，即HIV-1蛋白酶和HSP90，进一步突出了训练数据的稀缺性问题。

表1显示了已经用于预测药物-靶点停留时间的不同ML方法的摘要。最早发表的方法之一是QSKR（定量结构动力学关系）模型，使用主要是水的VolSurf描述符来预测37种HIV-1蛋白酶抑制剂的停留时间。另一种方法是使用COMBINE分析，在PLS（偏最小二乘法）模型中使用特定蛋白质残基与配体之间的静电和范德瓦尔斯相互作用作为特征（具有不同的权重）。另一种基于ML的方法，使用来自配体与HSP90解离的随机加速分子动力学轨迹的蛋白质-配体相互作用指纹（fingerprints）。通过这样做，不仅开发了一个支持向量（SV）回归模型，可以预测潜在的HSP90配体的停留时间，还可以指出配体与特定蛋白质残基的相互作用对延长停留时间的重要性，有助于指导基于结构的HSP90配体的药物设计。所有这些方法的一般精确度都在1个对数单位左右。这些研究主要是在小的数据集上验证的，使其真正的预测性难以评估。

表1 已发表的预测药物-靶点停留时间的ML方法

部署这些ML模型的能力，在很大程度上取决于所研究的特定蛋白质系统，因为只有当有足够的蛋白质动力学结合数据来训练ML模型时，我们才能预测化合物在蛋白质靶点上的停留时间。这就限制了这些方法在研究良好的蛋白质系统中的应用，并使其不太适合于药物发现，因为药物发现通常涉及针对新型蛋白质靶点（第一类靶点）开发药物。

2 材料

本节详细说明了需要安装的Python库，以执行这里详述的ML方法。

1.要安装PyQSAR，需要创建一个Python 2.7环境。

conda create -name py2 python=2.7

2.对于配体特征的生成，将使用Mordred。这必须安装在Python 2.7环境中。

3. 在单独的Python 3.7环境中，安装以下软件包：matplotlib, RDKit, pandas, 和scikit-learn。

3 方法

本章介绍了两种可以预测药物-靶点停留时间的ML方法。第一种是只用配体的方法，换言之，只用配体的特征来训练模型。第二种方法结合了配体和其蛋白靶点的特征。执行第二种方法需要结构数据或高质量的同源模型。

3.1 动力学结合数据

无论是哪种方法，都需要训练数据来训练监督下的ML模型。这需要在实验中确定指定靶点的配体停留时间值。人们可以通过搜索活性类型从ChEMBL获得配体动力学数据，然而，这些数据比结合亲和力终点的数据要少得多，而且往往很不全面。最近出版的一个数据库，KOFI-DB（http://koffidb.org），包含了从表面等离子体共振得到的配体结合动力学参数。目前在这个数据库中大约有1000个koff值。到目前为止，最大的动力学数据集合是KIND（KINetic Dataset），它包含了3812个条目，这些条目是由欧盟-IMI联盟K4DD（http://k4dd.eu）为广泛的靶点（包括离子通道、激酶和GPCRs）整理出来的数据。尽管这似乎令人印象深刻，但它比同等的结合亲和力数据库要小得多。对于单一靶点，配体结合的动力学数据是相当稀少的，一般少于100条。动力学数据的匮乏是使用ML来预测停留时间的主要挑战之一。

在本章所示的例子中，将使用从文献中搜集的GPCR数据和人工得到的数据。这536个条目大约是KIND中GPCRs条目数量的两倍，可以从https://potterton48.github.io下载。温度校正的停留时间值的计算是为了尽量减少不同研究小组在进行动力学检测时在不同出版物中使用的温度差异的影响。所有的停留时间值都用Arrhenius方程（见公式1）校正到所有实验中使用的平均温度（294.15 K）。由于温度的变化很小（低于15K），频率因子（A）被认为是常数。

公式1

确保训练数据对测试化合物有相对的代表性是很重要的。出于这个原因，NK1受体的肽的条目被从GPCR训练集的例子中删除。

3.2 单一靶点的仅有配体特征的QSKR

第一种方法，仅在配体特征上训练的ML模型，是一个单一靶点的QSKR多线性回归模型。一个开源的Python库，PyQSAR，将被用来进行QSKR建模。

1.创建一个有三列的表格：化合物名称、化合物SMILES字符串和相关的实验确定的、温度校正的停留时间。这些数据将被用作QSKR模型的训练和测试数据。

2.调查停留时间值的分布；它是否遵循正态分布？如果不是，可能要对数据进行转换以达到这个目的。在GPCR的例子中，为此目的使用了log10转换。

3. 在macOS/Linux上用"conda activate py2"启动python 2.7环境，或在Windows上用"activated py2"，以便使用PyQSAR。然后启动一个Jupyter笔记本会话。

4. 导入以下库：pandas, numpy, Mordred, RDKit, multiprocessing, pyqsar, 和scikit-learn。

5. 使用pandas库，将带有SMILES字符串和相关居住时间值的表格加载到Python 2.7。

6.为了使用 Mordred 生成化合物的特征，首先必须根据 SMILES 字符串为每个化合物创建 RDKit 分子。

mols = [Chem.MolFromSmiles(mol) for mol in df['SMILES'].values.tolist()]

这个命令循环查看DataFrame中 "SMILES"列的每个值，并从中生成一个RDCit分子，将其存储在一个列表中。

7.使用 Mordred 生成大约1500个配体特征。这可能需要一些时间，但在一个CPU上每秒大约可以处理10个化合物。Mordred是免费提供的，但是诸如Dragon等程序可以用来生成更多的配体特征。

8.删除没有生成数值的特征。当一个描述符不适用于某个特定的配体时就会发生这种情况。如果只有有限的配体对某一描述符缺失数值，可以用例如中位数的数据来填补。

9.使用Scikit-learn对每一列的数据进行缩放。

10.现在将数据分成训练集和测试集。通常情况下，80:20的分割是为了最大限度地增加用于训练的数据量。有几种不同的方法可以用来分割数据，最简单的是随机分割。

11. 对特征进行聚类，以找到高度相关的特征，减少需要使用遗传算法搜索的总体特征的数量。

12.使用遗传算法进行特征选择，只从一个聚类中选择一个特征，以防止高度相关的特征同时被选中。需要提供靶点信息和实验确定的停留时间。特征选择函数中的"components"参数决定了QSKR模型中所包含的最终特征数量（见注5）。这个步骤可能需要几分钟到几小时，取决于配体的数量和计算能力。该命令将返回具有最佳预测能力的所选特征的名称。

13.将所选特征的测试和训练数据都保存为.csv/pickle文件。同时将实验确定的、经温度校正的停留时间值保存为.csv/pickle文件。

14.打开Jupyter Python 3笔记本，通过关闭Python2.7环境 "conda deactivate py2"，然后重新启动Jupyter笔记本。加载以下模块：scikit-learn、pandas、matplotlib。

15.将保存的特征数据和靶点数据（实验确定的、经过温度校正的停留时间）加载到Python中的两个独立的pandas'DataFrames中。

16. 使用在特征选择阶段选择的特征在训练数据上训练多线性回归模型。使用该训练模型来预测测试集的值。将应用的质量指标是均方根误差（RMSE）和R2。对于A1受体的动力学数据，取得了以下结果。RMSETrain = 0.30, RMSETest = 0.48, R2Train = 0.76, R2Test = 0.67。

17.使用Matplotlib或任何其他绘图软件/软件包绘制结果（A1受体QSKR模型的结果见图2）。通过调查图上的异常值，可以确定模型在某些情况下失败的原因，以及可以做什么来改进它。

图2 预测和实验确定的停留时间的相关图

QSKR模型对A1受体动力学配体数据的预测和实验确定（Expt.）的停留时间（RT）的相关图。蓝色的、开放的圆圈是训练数据，橙色的、不透明的三角形是测试数据。

3.3 根据分子动力学模拟获得的特征训练的多靶点QSKR模型

在没有足够的数据开发单靶点QSKR模型的情况下，可以使用多靶点模型来增加数据量。多靶点模型需要蛋白质的一些信息，要么通过配体所针对的蛋白质的明确表征（例如，通过输入蛋白质序列），要么通过配体与蛋白质的相互作用。这些模型朝着预测配体动力学速率的更普遍化的模型发展，而不是只适用于单一受体或单一靶点的单一配体系列的模型。

为了将介绍中提到的COMBINE工作流程扩展到多靶点模型，可以使用蛋白质家族编号方案来寻找等价残基，作为训练数据输入模型。例如，对于GPCRs，可以使用GPCRdb修改的Ballesteros和Weinstein编号方案来寻找配体能量值（VdWs和静电）到具体残基位置（见图3）。

图3 使用蛋白质家族的编号方案来寻找等同的蛋白质-配体相互作用。

该图展示了如何使用蛋白质家族的编号方案（GPCRdb对GPCRs的编号方案），来寻找等价的残基来分配相互作用能量。这个例子显示了A1（PDB登录号：5UEN）和A2A受体（PDB登录号：3PWH）第一螺旋中的两个等效残基，分别为蓝色和桃红色。小热图显示了这两个残基位置的相互作用能量、范德瓦尔斯（VdW）和静电（Elect）的例子，深色代表较强的相互作用。

使用蛋白质-配体相互作用来训练多靶点模型的问题之一是获得足够的数据来训练模型。PDBbind数据库通过整理PDB中所有的蛋白质-配体结构和相关的结合亲和力数据来解决这个问题。由于缺乏配体动力学数据，配体动力学数据和相关的PDB结构之间的重叠度非常小。因此，人们不得不使用蛋白质-配体结构的预测结构来增加数据。为了增加这些预测姿势（对接）的可靠性，人们可以在对接后使用短分子动力学（MD）模拟的集合。在MD数据上训练模型已被证明在预测对数P等值方面有良好的表现。下面简要介绍一种开发以MD模拟数据为基础训练的多目标QSKR模型的方法。

1.获得蛋白质-配体的起始结构。倾向于使用X射线或冷冻电镜结构作为起始结构。如果不行，应该用对接法来预测配体的结合位置。请记住，所获得的数据的质量（输出）在很大程度上取决于起始结构的质量（输入）。

2.使用尽可能自动化的设置进行高通量的分子动力学模拟。有一些工具可以在这方面提供帮助，如HTMD。目标是进行集体模拟以确保结果的可重复性。这些集体模拟是复制的，唯一的区别是分配给每个原子的起始速度。

3. 从模拟中获取属性，如氢键、RMSD、RMSF、蛋白质-配体相互作用指纹。VMD、Chimera或MDAnalysis可以帮助提取这些特征，努力获得适合于系统和要建模的特征。

4.应用几种不同的ML方法，使用本章中概述的分割策略。使用验证集或通过k-means验证来评估哪种方法（和超参数）是最好的。

4 注释

1.PyQSAR 是一个Python 2.7库。由于Python 2.7已经通过了它的"生命末期"，因此不再得到适当的支持，只有需要的工作才应该在Python 2.7中进行(因为它调用了PyQSAR)。Python 3应该作为其他一切的默认使用。为了尽可能容易地支持这些不同的Python环境，应该使用Anaconda。

2.所有发表的GPCR-配体动力学数据被收集到一个数据库中。从每个公布的数据中获得的主要数据是配体名称、SMILES字符串、动力学测定的温度、Kon、koff、KD和Ki。室温被假定为294.15 K。停留时间被计算为koff的倒数。

3.可以使用配体的其他特征，如扩展连接性指纹（ECFP）。可以测试多个描述符，看哪个描述符对有关任务的性能最好。

4.对于任何一种拆分，重要的是要确保目标值（本例中的停留时间）在训练集和测试集中都有类似的分布。如果是随机拆分，由于测试集中的化合物在结构和相关目标值上与训练集非常相似，很容易高估模型的能力。为了避免这种情况，可以进行将配体结构考虑在内的分割。基于时间的拆分是另一种选择，可以重现药物发现项目的情景。

5.大多数QSAR模型有少于十个特征。在所示的例子中，已经选择了四个，但特征的数量可以变化，以确定什么是最适合一个给定的系统。一般来说，目标是使用最小数量的特征，在训练数据中产生相当好的准确性，以减少过拟合的机会。理想情况下，如果数据允许，可以使用验证集来研究模型中应该包含多少个特征以获得最佳性能。

6.RMSE的优点是与目标数据的单位相同，所以更容易理解。R2显示了模型的表现是比随机预测好还是差。训练和测试RMSE之间的差异越大，越有可能发生过拟合。

7.对于GPCR QSKR模型，计算了以下特征，因为它们有一定程度的实验证据表明它们有助于停留时间：水和配体之间的相互作用能量，GPCR第二细胞外环的RMSD，结合点的RMSD，以及配体大小的测量。

参考资料

https://link.springer.80599.net/protocol/10.1007/978-1-0716-1787-8_8

----------- End -----------