JCIM｜基于相对结合自由能建模的主动学习指导先导化合物优化

2023年1月4日，来自美国卡内基梅隆大学的研究者在Journal of Chemical Information and Modeling上发表文章“Active Learning Guided Drug Design Lead Optimization Based on Relative Binding Free Energy Modeling”。

作者开发了一个主动学习模型，能够从巨大的分子库中筛选出结合亲和力改善的化合物，其仅通过少量次数的热力学积分计算相对结合自由能，效率高，而且有效。作者将模型应用于SARS-CoV-2木瓜蛋白酶样蛋白酶的抑制剂的筛选，找到了133种结合亲和力改善的化合物。

1 摘要

使用计算方法进行有效蛋白抑制剂的识别通常需要预测配体结合自由能(a ligand binding free energy, BFE)。基于分子动力学模拟(MD)的热力学积分(Thermodynamics integration, TI)方法是一种能够获得精确BFE的计算方法，但其计算成本高且耗时。在这项工作中，作者开发了一种高效的自动化工作流程，用于识别数千种同源配体中BFE最低的化合物，这本应需要数百次热力学积分计算。自动机器学习(Automated machine learning, AutoML)由主动学习(active learning)实现，可以无偏、高效地搜索性能最佳的分子数据集。

作者应用此工作流程来选择SARS-CoV-2木瓜蛋白酶样蛋白酶(SARS-CoV-2 papain-like protease, SARS-CoV-2 PLpro)的抑制剂，并找到133种具有改善的结合亲和力的化合物，其中包括16种具有改进超过100倍结合亲和力的化合物。作者获得的苗头化合物比率(hit rate)超过了传统药物化学家的方法。这证明了AL和AutoML与自由能模拟的结合提供了与朴素暴力方法相比至少20倍的加速。

2 引言

药物设计的苗头化合物到先导化合物(hit-to-lead)和先导化合物优化(lead optimization)阶段旨在通过改变苗头化合物的化学结构来发现先导化合物。典型的先导化合物优化过程包括首先化学合成多种化合物，然后测试其生物活性。这昂贵且耗时。基于结构的超大分子库虚拟筛选，其目的是尽量减少实验室合成和测试所选化合物的数量，已成为计算药物设计中的一种成功策略。尽管通过将配体与靶蛋白对接实现了高苗头化合物比率，但此类方法的两个主要局限性仍然存在：对接方法预测配体结合亲和力(ligand binding affnity)的能力有限，以及筛选数十亿化合物组成的库在技术上存在困难。

与对接方法不同，全原子分子动力学模拟方法(包括热力学积分(thermodynamics integration, TI)方法)可以高精度地预测配体结合亲和力，也称为结合自由能(binding free energy, BFE)。相对结合自由能(relative BFE, RBFE, ΔΔG)定义为新配体和原先导化合物之间的BFE差异。虽然GPU加速MD模拟的算法得到了改进，但为大量化合物计算多个RBFE仍然非常耗时且技术上难以处理。

因此，为了克服这个问题，作者开发了一种基于TI计算的RBFE的机器学习方法(即主动学习引导的先导化合物优化过程的自动化方法)。用于TI计算的化合物是用自动ML算法选择的，该算法旨在实现两个目标：(1) 选择良好的结合物用于TI计算，(2) 使用TI计算的RBFE改进ML模型对整个分子筛选库的RBFE的预测。为了实现这一双重目标，作者将TI RBFE计算与自动机器学习(AutoML)循环相结合，从而消除了模型选择偏差，并有效地利用了每次AL迭代的信息增益。这种AutoML–TI RBFE计算工作流程允许通过最少的TI RBFE运算识别紧密结合的配体。

图1 SARS-CoV-2 PLpro及其抑制剂的结构

(A) SARS-CoV-2 PLpro与GRL0617(PDB ID: 7JIR)复合物的结构。催化三联体的抑制剂和残基相应地显示为粉红色和绿色。(B) GRL0617，普通骨架(N-[(1R)-1-arylethyl]arenecarboxamide)以绿色显示。设计SARS-CoV-2抗病毒药物最吸引人的药物靶标之一是SARS-CoV-2木瓜蛋白酶样蛋白酶(SARS-CoV-2 PLpro)，这是一种负责处理病毒多聚蛋白并抑制宿主免疫功能的酶。PLpro有315个残基，由两个不同的结构域组成：一个小的N-末端泛素样结构域和一个“拇指-手掌-手指”催化结构域(见图1)。

背景知识：据悉，SARS-CoV-2木瓜蛋白酶样蛋白酶PLpro是加工病毒多蛋白来产生功能性复制酶复合物、并使病毒能够传播所需的必需冠状病毒酶。PLpro还涉及切割宿主蛋白上的蛋白质翻译后修饰，从而作为针对宿主抗病毒免疫反应的逃避机制。感染后，SARS-CoV-2-PLpro有助于从干扰素反应因子3(IRF3)切割SG15，并减弱I型干扰素反应。重要的是，用GRL-0617抑制SARS-CoV-2-PLpro会削弱病毒诱导的细胞致病作用，促进抗病毒干扰素途径并减少感染细胞中的病毒复制。这些结果突出了一种双重治疗策略，即靶向SARS-CoV-2-PLpro可以抑制SARS-CoV-2感染并促进抗病毒免疫。(来源于网络)

3 方法

3.1 数据集筛选和分子准备

图2 SMARTS过滤器用于小分子数据集挖掘。

(A) 步骤1，SMARTS模式的可视化，用于N-[(1R)-1-arylethyl]arenecarboxamide衍生物的非立体特异性搜索(non-stereospecific search )。(B) 步骤3，SMARTS模式的可视化，用于N-[(1R)-1-arylethyl]arenecarboxamide衍生物的立体特异性搜索(stereospecific search )。芳香族原子表示为虚线圆圈，脂肪族原子表示为实心圆圈。深灰色代表碳原子，红色代表氧原子，蓝色代表氮原子，浅灰色代表任何原子。

作者采用了三家化学品供应商(包括WuXi，Mcule和Enamine)的目录合并为一个数据集，总计约13亿个分子。数据的处理共计6步(图2)：

1.基于SMARTS模式(图2A)的分子虚拟筛选，使用Bemis–Murcko骨架，允许碳取代氮。

2.使用OpenEye QUACPAC和OEFlipper工具包枚举质子化、互变异构和立体异构状态，允许至多8个手性中心：仅保留中性分子用于进一步计算。

3.使用SMARTS模式(图2B)过滤：仅保留手性中心构型与参考配体(R-立体异构体)相同的分子。

4.使用OpenEye Omega工具包为过滤得到的9998个分子生成3D构象。

5.对生成构象分子，进行分子对接和姿态过滤。

作者之后介绍了分子动力学模拟的原理、方法和实验设置，具体可以参见原文。

3.2 机器学习模型开发和特征工程

特征工程与分子表示

作者构建了5种不同的分子表示：

1.RDKit分子指纹(path fingerprints with path length 7，2048维)。

2.Morgan指纹(Extended-Connectivity Fingerprints with radius 3，2048维).

3.3D分子指纹E3FP.

4.蛋白-配体扩展连接指纹(protein−ligand extended connectivity (PLEC) fingerprints)。

5.E3FP和PLEC指纹的拼接。

然后作者采用3种降维方法，包括PCA、MDS和TSNE，4个降维后的维度2、10、100、200，最后组合得到42种分子表示(见原论文附录表S1)。

自动化机器学习(automated machine learning, Auto ML)

作者将上述的42种分子指纹作为模型输入(视作超参数)，并采用Scikit-learn实现的7类机器学习算法，包括：Random Forest (RF)、Multi-Layer Perceptron (MLP)、Linear Regression (LR)、KNeighbors Regression (KNR)、SupportVector Regression (SVR)、Gaussian Process Regression (GPR)和GaussianProcess Regression with Tanimoto Kernel (GPT).

模型选择和超参数搜索在每个主动学习周期从头开始，并以5折交叉验证的平均绝对误差(MAE)作为模型选择标准。

主动学习(Active Learning, AL)

图4 在一个主动学习周期中的自动计算工作流的一般方案。工作流程包括2个主要模块(AutoML和TI RBFE)以及4个主要步骤。步骤1中计算ΔΔG的分子被描绘为彩色六边形。化学空间显示为2D t-SNE图。以周期0为例，介绍整个计算流程如下：Start of Workflow，采用Butina聚类方法将8175个分子聚类为45个簇，并选择质心分子进入步骤1，利用TI RBFE模块进行ΔΔG计算。然后得到累积的带标签数据(Cumulative Label Data)，接着进入步骤2，利用AutoML模块，进行模型选择和超参数(不同分子表示)搜索。之后，进入步骤3，根据最优的模型给所有分子打上标签。进入步骤4，选择样本，进行下一轮的利用TI RBFE模块进行ΔΔG计算。至此构成一个主动学习周期，8轮后，在End of Workflow处停止：选择带标记数据中的冠军作为优化后的分子。

主动学习过程是迭代进行，单个主动学习周期的过程见图4，整体流程如下：

1.从选择的分子的种子集开始，进行了TI RBFE计算以训练初始ML模型。

2.然后使用当前ML模型为下一轮TI RBFE的计算选择分子。

3.我们为步骤2中选择的分子利用TI计算了RBFE，并使用更新的具有RBFE标签的数据集重新训练ML模型。重复循环直到收敛。

整个主动学习由8个主动学习周期组成，具体见原论文附录表S2。

4 结果

4.1 主动学习8周期的过程分析

图5 (顶部)展示了主动学习8个周期的热力学积分计算的RBFE；(中部)10折交叉验证衡量ML模型的MAE；(底部)追溯绝对误差(retrospective absolute error)。

图5显示了所有AL循环中获得的TI RBFE。AL循环0采用多样选择(diverse selection)方式进行选择分子的初始化，以尽可能广泛地对所选分子库的化学空间进行采样。接着，为这组初始分子计算TI RBFE，并将其提供给AutoML模块，用于初始ML模型训练。对于接下来的5个AL周期(AL周期1–5)，作者采用平衡选择(balanced selection )方式进行选择分子。这5个AL循环的目标是获得关于所选分子库的化学空间的信息，而不是选择具有最低RBFE的分子。

随着AL循环的进行，ML模型的性能提高。10折交叉验证的平均绝对误差(MAE，图5，中间)达到1 kcal/mol。为了验证模型的收敛性，作者使用随机选择(random selection)的分子进行了第6次AL循环。分子的随机选择也有助于克服AL被困在化学空间的局部最小值中的可能问题。最后在周期7，作者采用了贪婪选择(exploitative selection)方式选择分子，并计算RBFE，得到实验结果。

图7 RBFE结果和主动学习周期内的ML模型演变。每个子图(由主动学习AL编号和相应的选择类型标记)显示分子的两个2D t-SNE图：(顶部)在各自的AL循环中选择用于热力学积分的RBFE计算的分子由热力学积分计算的ΔΔG着色，其余以灰色显示；(底部)分子由ML预测的ΔΔG着色。ΔΔG值的色条显示在底部。

图7显示了ML模型对分子的化学空间的感知的演变过程，以及为热力学积分计算RBFE选择的分子的分布。值得注意的是，在主动学习工作流(AL周期0)开始时，ML模型没有区分化学空间中富含低ΔΔG的有利结合物的特定区域(选择的分子广泛随机分布在化学空间的各个区域(顶部))。在随后的AL循环1–5中，采用平衡(balanced)选择，模型探索了多个区域(图7，AL 1–5，TI行)。随着周期迭代进行，ML模型的感知正在发生显著变化。ΔG分子密集分布的化学空间区域开始被识别。可以看到化学空间的各个区域在主动学习的第5个循环开始，着色逐渐稳定，ML模型收敛(图5亦能够证明)。在AL周期6(图7，AL 6，TI行)中，采用随机选择的分子，采样分子分布在化学空间中，大多数分子的ΔΔG为正。值得注意的是，模型的误差(图5)没有增加，这支持模型已经收敛的观察。因此，如上所述，在周期7中选择最低的30个分子即可结束主动学习过程，得到研究结果。

4.2 具有改进的预测结合亲和力的配体分析

图8 结合亲和力提高的配体。(A) ΔΔG为负的配体的公共骨架。“Ar”对应于含有六元芳香环的任何芳香体系。关于参考配体的骨架的化学修饰被圈出。(B) 与A中公共骨架对应的参考配体类似物。(C) 具有最高预测结合亲和力的配体。

N-[(1R)-1-arylethyl]arenecarboxamide的萘环中的两种常见修饰存在于预测结合亲和力提高的分子中(图8A中的骨架S1和S2)。第1种修饰(S1)是在萘环的4位用氟取代氢。第2种修饰(S2)包括将β-萘碳取代为氮，并在芳香环的7位添加甲氧基。为了评估这些修饰的相对重要性，作者计算了GRL0617→ M1的RBFE和GRL0617→ M2的RBFE(图8B)，发现结合亲和力分别提高−0.84和−0.99 kcal/mol。

具有改进的预测结合亲和力的配体的第3个常见结构特征是存在融合的5,6-和6,6-双环芳香族体系取代参考配体的苯环(图8C)。在ΔΔG为负的配体中，有35个(约26%)的分子具有类似的芳香体系。其中9个分子的预测结合亲和力提高了100倍以上(ΔΔG<2.73 kcal/mol；见附录图S2)。对于预测ΔΔG最高的配体(配体1–3，见图8C，此处的最高似乎应指ΔΔG数值上最小)，计算出的ΔΔG分别为−4.06、−4.05和−3.72 kcal/mol，对应于2.85、2.90和5.07 nM的解离常数。

图9 (A)参考配体、(B)配体1、(C)配体3的代表性结合姿势。配体和蛋白质残基的碳原子分别以绿色和灰色显示。氮、氧和氟原子分别以蓝色、红色和青色显示。

参考配体苯环取代基显示出与蛋白质的特异性相互作用：氨基与Gln269的酰胺基和Tyr268的羟基形成氢键，甲基与Tyr264、Tyr273和Leu162的侧链形成疏水性相互作用(图9A)。配体1和3与参考配体的代表性结合姿势如图9B、C所示。连接头(linker)的酰胺基与Gln269的主链氨基、Tyr264的羟基和Asp164的羧基形成氢键。参考配体和配体3的萘环以及配体1和2的异喹啉环与Pro248和Tyr268的侧链形成疏水相互作用。参考配体和配体2、3的苯环以及配体1的吡啶环与Gln269和Asp164侧链的脂族区域形成疏水相互作用。

论文中还进行了其他详细的实验及结果分析。

5 结论

先导化合物优化仍然是现代计算化学的一个重大计算挑战。计算密集的任务，例如相对结合自由能模拟的分子动力学，通常受到计算资源的可用性以及以高通量方式执行计算的困难的严重限制。例如，战胜COVID-19的登月计划利用全球的Folding@home(一个分布式计算软件)计算方案运行了5000多次自由能模拟。这项大规模的工作使用了数亿计算机小时，使得在寻找SARS-CoV-2主要蛋白酶的效力提高了100倍。这种资源很难获得。在这里，作者能够仅对配体的子集进行相对结合自由能计算，而不是通过将此类计算与主动学习方法(包括自动机器学习模型选择)耦合，来对先导化合物的所有可用类似物进行相对结合自由能计算。

使用Auto ML程序对分子进行选择，作者通过仅对253个配体进行TI RBFE计算，鉴定了133种潜在的预测其具有改善的结合亲和力的SARS-CoV-2 PLpro抑制剂。值得注意的是，由ML选择的70%配体的结合亲和力提高，而随机选择的配体只有11%提高了结合亲和力。本文开发的方法是通过利用现代计算方法加快药物设计中先导化合物优化的重要一步。

参考文献

Gusev F, Gutkin E, Kurnikova M G, et al. Active learning guided drug design lead optimization based on relative binding free energy modeling[J]. Journal of Chemical Information and Modeling, 2022.

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