“最后一步难题”：基于物理方法和人工智能的工具在虚拟筛选中小分子结合预测中的效果评估

《Journal of Chemical Information and Modeling》：The Last Mile Problem: A Critical Assessment of Physics-Based and AI Tools for Small Molecule Binding Prediction in Virtual Screening

　　高分辨率图像下载
MS PowerPoint 幻灯片

基于对接的虚拟筛选（VS）在药物或探针发现的初始阶段对于寻找潜在活性分子至关重要。然而，该方法仍然容易产生较高的假阳性率，这往往导致筛选活动失败。基于分子动力学（MD）的炼金术自由能方法为提高 VS 的命中率提供了有希望的解决方案，但它们对资源需求极高。结合炼金术绝对结合自由能（ABFE）计算的现实世界和基准研究可以帮助优化其在 VS 流程中的使用。在这里，我们提出了一个大规模的基准测试，以评估 ABFE 计算在 VS 工作流程中的相对价值。我们使用了两个数据集：一个是来自 PDBbind 数据库的 632 个配体-蛋白质复合物的精选集，用于评估 ABFE 的定量准确性；另一个是来自有用诱饵数据库（DUD-E）的 315 个结合子和诱饵的集合，用于评估其在 VS 环境中的预测能力。除了炼金术 ABFE，我们还对比了计算成本较低的全局状态物理方法以及五种机器学习（ML）模型。研究发现，预测器的排名与其计算成本一致，其中炼金术 ABFE 在两个基准测试中表现均较好。全局状态方法在识别 DUD-E 数据集中的活性分子方面表现良好，但在 PDBbind 中与实验值的关联程度较低。大多数 ML 模型在 PDBbind 上表现不错，这可能是由于训练数据的重叠，但在 DUD-E 上表现不佳，除了 GNINA 和 Boltz-2，它们表现出与全局状态物理方法相当的一般化能力。总体而言，采用 Boltz-2 作为初级筛选器，随后使用炼金术 ABFE 的分阶段方法，可以稳健且经济高效地丰富基于对接的 VS 活性分子列表。

引言
多年来，基于对接的虚拟筛选（VS）已成为启动药物发现项目的主要方法。此外，对于缺乏大量传统化合物收藏的小型制药公司和学术实验室来说，VS 是不可或缺的。与实验筛选方法不同，VS 不需要实际获得化合物；只需要获取一小部分经过计算优先排序的候选分子进行实验验证。现在有许多公开的对接算法（3?5），在性能和结合位点准确性方面都显示出价值。最近，开源的基于 ML 的方法扩展了对接工具的数量，带来了实质性的改进（6?8）。尽管取得了这些进展，但在计算机模拟筛选中仍然容易产生较高的假阳性率，这可能导致筛选活动失败——这是许多虚拟筛选者可能遇到的问题。直到最近，公共领域中关于虚拟筛选效率的定量估计仅限于已发表的案例研究。根据综合评论，实验确认的命中率通常超过 20%（9,10）。然而，这些估计很可能会受到生存者偏差的影响，即失败的筛选活动几乎从未被公开发表。因此，公共领域中的成功率看似很高，但实际上是因为社区只能接触到那些对配体有预先了解的目标的非常成功的虚拟筛选结果。为了更客观地了解 VS 的效率，2021 年启动了“计算命中发现实验的关键评估”（CACHE）挑战（10）。每个挑战的参与者需要针对特定目标进行 VS，并提交 100 个虚拟候选分子的列表以供采购和实验验证。CACHE 挑战 #1 的结果显示，23 名参与者中只有 7 人识别出了任何实验确认的活性分子（11），而在 CACHE 挑战 #2 中，这一数字仅为 5 个（12）。是什么因素促成了获胜策略的成功？有趣的是，成功的因素并不在于对接算法的选择，因为大多数参与者（包括获胜者）依赖于广泛使用的软件，如 AutoDock-GPU（13）、AutoDock Vina（14）或 Glide（15）。相反，成功似乎源于 VS 工作流程的其他方面，例如在对接前对配体的选择性预筛选（16）、改进的评分函数（6），或广泛使用基于 MD 的绝对和相对结合自由能（ABFE/RBFE）计算（17）。基于 MD 的炼金术自由能方法通常被认为是提高 VS 命中率的最终解决方案，有可能在少量实验测试过的候选分子中识别出真正的活性分子（18,19）。

炼金术自由能计算提供了一个严格的、基于物理的框架来估计结合亲和力（19?22）。炼金术方法通过将结合在蛋白质上的配体缓慢转化为无价物的过程来评估 ABFE。通常，这种配体湮灭是通过称为 λ-窗口的离散的、化学上不可能的步骤进行的。可以通过对所有采样配置的系统能量差异进行统计分析来计算连续炼金术状态之间的自由能。然后对这些炼金术能量求和，得到最终的 ABFE 值。同一 λ-窗口内状态之间的高相似性确保了两种状态具有相等的配置空间，这是该方法适用的前提。然而，这些计算并不完全准确。误差来源包括分子力学力场（这只是真实量子相互作用的近似）、窗口分辨率不足以及对配置空间的采样不足。此外，这些计算在计算上非常昂贵，需要大量的时间和资源。对于典型的 VS 研究，实验室可能只能负担得起最多 100 个化合物的 ABFE 计算。然而，100 个对接候选分子的列表可能根本不包含任何活性分子，整个列表都需要通过实验来验证才能获得可靠的结果。为了充分利用 ABFE 计算，必须将其应用于更大的分子集——最多可达 1,000 个分子——以确保实验测试的活性分子列表较少。即使计算资源继续变得更加强大和便宜，但在可预见的未来，对数百个配体-蛋白质复合物进行炼金术模拟的成本仍然可能很高。这样的资源密集型努力能否产生预期的结果？虽然支撑炼金术计算的理论非常严谨，但实际应用可能会受到不完美的力场或 λ-窗口数量不足的阻碍。此外，某些配体的错误对接姿势可能导致错误的 ABFE 预测。到目前为止，对于如何有效地应用资源密集型方法来获得用于实验验证的丰富活性分子列表，还没有明确的共识。此外，这最后一步中每个化合物的成本比所有先前步骤的成本高出几个数量级，几乎没有浪费的空间——这在物流领域通常被称为“最后一公里问题”。结合 ABFE 计算的现实世界和基准研究可以促进其在 VS 流程中的有效使用。

为了实现这一目标，我们提出了一项大规模的基准研究，据我们所知，这是首次评估 ABFE 计算在 VS 工作流程中的相对价值的研究。在这项研究中，我们首先通过执行亲和力预测来建立准确性基准，使用了来自 PDBbind 数据库的 632 个配体-蛋白质复合体（包括 206 种独特蛋白质和 606 种独特配体），这些复体的结合亲和力（Kd/Ki 值）和 3D 结构已经通过实验确定（23）。其次，我们评估了来自增强型有用诱饵数据库（DUD-E）的 135 个真实结合子和 180 个诱饵的亲和力（24）。应用了超几何概率分布函数（25）来估计每种基准技术在现实世界 VS 研究中可能实现的富集程度。我们总共评估了三种基于物理的方法：（i）炼金术 ABFE（22,26）、（ii）分子力学/泊松-玻尔兹曼表面积（MM/PBSA）（27）和（iii）MM/广义玻恩表面积（MM/GBSA）（27），以及七种机器学习（ML）模型：GNINA（6）、KDEEP（28）、OnionNet-2（29）、TopologyNet（30）、Yuel（31）、RF-score-VS（32）和 Boltz-2（33）。MM/PBSA 和 MM/GBSA 是全局状态方法，只需要对配体-蛋白质复合物和未结合的配体进行 MD 模拟。虽然全局状态方法每次亲和力计算的计算时间与炼金术 ABFE 技术相当（以小时计），但它们使用的计算节点数量少了一个数量级，从而实现了显著的并行化。因此，评估这些资源较少要求的替代方法是至关重要的。最后，我们还包括了基于 ML 的亲和力预测器作为更经济的选择。自从深度神经网络的出现（34）以来，ML 模型在亲和力预测方面显示出潜力。选择这五种 ML 模型是基于它们不同的架构、报告的性能以及源代码的可用性。

结果 数据集
本研究使用了两个数据集。其中一个数据集是 PDBbind 基准集，其中包含具有实验确定的亲和值和 X 射线结构的配体-蛋白质复合物。这是一个参考数据集，允许我们评估基准方法在定量预测已知结合体的配体-蛋白质结合自由能方面的准确性。第二个数据集是 DUD-E 基准集，用于评估方法使用对接姿势作为 MD 模拟和需要 3D 输入信息的 ML 方法区分结合体和非结合体的能力。这些基准方法在这一数据集上的表现将表明它们在 VS 环境中的实用性。参考亲和力数据集的主要要求是具有实验确定的 3D 结构和以 Kd/Ki 表达的亲和力数据。精制的 PDBbind 子集包含 5,316 个由 4,131 个独特配体和 1,409 个独特蛋白质组成的配体-蛋白质复合物（35）。由于这项研究的最终目标是评估 BFE 预测器在寻找活性分子方面的潜力，我们专注于那些可能作为药物发现项目先导药物的配体。为此，我们保留了符合 Lipinski 规则的配体（36），除了分子量限制在 250–600 Da 范围内。我们还排除了具有多个可电离基团的配体、含有磷原子的配体，以及具有潜在毒性或高反应性基团的配体（详见方法部分）——这些过滤器确保基准化合物具有与可能被认为是药物发现项目可行先导分子的分子相似的物理和化学特性。具有多个可电离基团或磷原子的化合物通常表现出较差的药代动力学特性，如低膜通透性，因此不适合作为先导优化的起点。此外，我们还排除了涉及金属蛋白的复合物，因为标准力场无法准确表示金属配位相互作用（37?39）。最终，保留了 632 个配体-蛋白质复合物的数据集用于所有基准计算（详见补充图 S1，其中包含数据集化合物的关键物理特性的概率密度函数（PDF），补充数据集 1 包含压缩的蛋白质结构 PDB 文件，补充数据集 2 包含 SMILES 格式的独特配体结构）。这些蛋白质结构属于 124 个蛋白质结构超家族（如 SUPFAM 数据库中定义的蛋白质激酶类或甲基转移酶类（见补充表 S1 中的蛋白质家族列表）（40））。每个蛋白质与 1 到 56 个配体相关联（详见补充数据集 3 中的每个目标的配体数量），九个目标与超过十个配体相关联。具有高配体数量的目标特别有助于确定是否可以更准确地预测单个目标的结合自由能，而不仅仅是整个数据集。最后，我们使用 Butina 聚类算法估计了化学多样性，Tanimoto 系数的纳入标准设置为 0.5（基于 ECFP4 指纹）。共获得了 379 个簇，簇的大小从 1（263 个单例）到 10（有五个簇包含超过 5 个化合物）不等。

结合体/非结合体数据集是从 DUD-E 数据库中选取的（24,41）。这样做的目的是在接近虚拟筛选活动条件的情况下测试基准方法。因此，需要统计上足够的已知结合体和诱饵数量 against 单一目标进行评估（这里的“诱饵”指的是具有与结合体相似物理特性的非结合体）。为了满足这一要求，同时保持炼金术 ABFE 计算所需的化合物总数在可承受范围内，目标的数量必须显著少于 PDBbind 数据集。之前，我们创建了一个包含 10 个结构多样性的 DUD-E 子集，广泛覆盖了蛋白质空间（42）。其中九个目标——ACE、ADRB1、FAK1、GRIK1、HMDH、MCR、PGH2、PRGR 和 TRYB1——也存在于参考 PDBbind 数据集中，并涵盖了最重要的药物类别，包括激酶、蛋白酶、GPCR、核受体和合成酶。该数据集共包含315个配体复合物——135个结合剂（每个目标15个）和180个诱饵（每个目标20个）。重要的是，这个DUD-E子集中大多数配体-蛋白质复合物的实验确定的3D结构以及定量亲和力数据都不可用，这是它的另一个特点，使其更适用于基于对接的虚拟筛选（VS）环境。本研究使用了AutoDock Vina方法来生成对接姿势。（14）DUD-E数据集的基于对接的配体-蛋白质复合物的压缩PDB文件集作为补充数据集4提供，而独特的配体结构则以SMILES格式作为补充数据集5提供。

PDBbind基准ABFE计算
我们首先通过使用理论严谨的abchemical ABFE计算为PDBbind数据集建立了一个参考基准。我们使用了GROMACS中实现的标准ABFE协议。（43,44）abchemical状态之间的自由能差异使用BAR方法估计。（45）abchemical ABFE协议中的关键参数是λ窗口的数量和每个λ窗口的大小。λ窗口的数量设置为25，这是近期出版物中常见的选择。（21,46?48）通过在从PDBbind数据集中随机抽取的77个具有均匀分布的实验结合亲和力的配体-蛋白质复合物上运行ABFE模拟来确定λ窗口的长度，窗口大小分别为1 ns、5 ns、10 ns和15 ns。由于abchemical计算的计算成本与λ窗口长度成正比，因此确定一个最佳的成本-准确性平衡长度非常重要。根据计算值和实验BFE值之间的皮尔逊相关系数（r）和斯皮尔曼等级相关系数（ρ）以及相应的散点图（图1A–D），可以预期使用较长的λ窗口的ABFE计算会更加准确。然而，考虑到这种计算所需的额外计算能力，10 ns和15 ns窗口的准确性提升相对有限。因此，整个研究中都使用了5 ns的λ窗口。

图1
图1. PDBbind数据集上abchemical（ΔGFEP）和实验（ΔGexp）BFE之间的相关性。分别在1 ns、5 ns、10 ns和15 ns的λ窗口长度下对一个试点子集（A–D）进行的模拟；以及带定向约束的设置（E），并对整个数据集（F）进行的模拟。

高分辨率图像
下载MS PowerPoint幻灯片

我们还研究了应用空间约束的效果，这些约束可以保持配体接近其初始结合姿势。通常建议使用这样的约束以提高热力学可逆性。因此，我们感兴趣的是在模拟环境中它们是否真的必要。为此，我们对用于窗口大小优化的相同77个复合物进行了带有定向约束的额外ABFE计算（图1E）。虽然应用约束后相关指标有所改善（皮尔逊r从0.73增加到0.77），但这种改进必须与VS分级的实际权衡因素相权衡。特别是，实施空间约束需要额外的λ窗口，这会将计算成本增加高达30%。此外，在高度多样化的配体和结合口袋之间自动化分配这些约束会在计算流程中引入显著的复杂性和潜在的失败点。因此，由于无约束设置仍然取得了相当高的相关系数r = 0.73，我们认为它提供了准确性、计算速度和自动化易用性的最佳平衡。这些实际考虑因素，加上先前研究表明通常可以在没有约束的情况下达到类似的准确性（20,49），促使我们决定对整个数据集使用无约束设置。

对于PDBbind数据集中的所有复合物，共运行了632次模拟。所有PDBbind复合物的计算BFE与实验BFE的散点图显示在图1F中（所有计算和实验BFE值作为补充数据集6提供）。整个数据集的皮尔逊r = 0.69和斯皮尔曼ρ = 0.68略低于在77个复合物试点子集上获得的值（r = 0.73和ρ = 0.70）（图1B）。我们试图了解这种下降的原因，以及BFE计算准确性如何随蛋白质或配体的性质而变化。这样的分析可能有助于预测abchemical ABFE计算在特定VS项目中的成功率。首先，我们检查了目标内相关性——定义为仅在结合同一蛋白质的不同配体集合之间计算的相关性——对于至少有10个配体的目标。我们发现，exp./calc.相关性在蛋白质之间有所不同，皮尔逊r的范围从0.44到0.75，斯皮尔曼ρ的范围从0.19到0.74，表明计算准确性取决于蛋白质-配体复合物的性质（补充图S2，补充表S2）。几个可能影响ABFE计算准确性的因素（并且可以预先考虑）包括蛋白质或配体的大小（通过分子量测量）、配体浸入蛋白质的比例（近似为每个配体原子的附近蛋白质原子的中位数），以及结合口袋的组成（表示为口袋残基的平均clogP）。我们测试了具有上述因素高于和低于中位数值的复合物的ABFE计算准确性是否有所不同。然而，所有分割子集的ABFE值与实验数据的相关性水平相似，因此不能用作预测ABFE预测准确性的有用标准（补充图S3）。

此外，我们还研究了abchemical方法的准确性是否取决于配体的亲和力。为此，我们计算了前50%结合剂（RMSE = 2.30 kcal/mol）和后50%结合剂（RMSE = 2.16 kcal/mol）的均方根误差（RMSE）。这两个子集之间的差异不显著，表明较强的结合剂并没有比较弱的结合剂模拟得更准确。

除了理论严谨的abchemical计算外，我们还探索了两种流行的基于实验的模拟方法的值。MM/PBSA（27）和MM/GBSA（27）都是末端状态技术；也就是说，MD模拟和随后的分析仅涉及结合的以及完全解耦的配体和蛋白质。MM/PBSA和MM/GBSA都使用了5 ns的MD轨迹。末端状态计算在三种模式下进行：（i）不进行熵校正，（ii）通过正常模式分析（NMA）进行熵校正（50,51），以及（iii）通过相互作用熵（IE）进行熵校正（52），后者比NMA计算成本更低。两种末端状态方法的有无熵校正的结果总结在图2A–F中（所有计算值作为支持信息数据集6的一部分提供）。总体而言，没有一种方法表现出令人满意的性能，并且与早期的证据一致，（53）经过熵校正的方法甚至比原始的MM/PBSA和MM/GBSA表现得更差。

图2
图2. 实验BFE（ΔGexp）与末端状态方法之间的相关性：（A）MM/GBSA（ΔGMM/GBSA），（B）MM/PBSA（ΔGMM/GBSA），（C）带有IE校正的MM/GBSA（ΔGMM/GBSA–IE），（D）带有IE校正的MM/PBSA（ΔGMM/PBSA–IE），（E）带有正常模式校正的MM/GBSA（ΔGMM/GBSA–nmode），以及（F）带有正常模式校正的MM/PBSA（ΔGMM/PBSA–nmode）。

高分辨率图像
下载MS PowerPoint幻灯片

目标内相关性，除了少数例外，与整个PDBbind数据集中观察到的结果一致（补充表S2）。MM/GBSA数据中只有两个目标（BRD4和WDR5）的皮尔逊r大于0.5，而MM/PBSA数据识别出三个这样的目标：BRD4、MK14和WDR5。

此外，我们还研究了末端状态方法的准确性是否取决于配体-蛋白质复合物的特定属性，包括蛋白质和配体的分子量、配体浸入蛋白质的比例以及结合口袋的亲脂性。如补充图S4和S5所示，这两种末端状态方法的准确性似乎不受这些因素的影响。

七个ML模型被选出来，以代表自重新关注将机器学习应用于化学以来开发的各种方法。例如，KDEEP（28）使用基于3D结构的3D卷积神经网络（CNN）进行训练；化学结构被描述为药效团特征。OnionNet-2（29）是一个图CNN；输入图节点是配体原子和整个蛋白质残基。TopologyNet（30）是一个1D CNN模型，其中配体和蛋白质由拓扑指纹表示。Yuel（31）是CNN（用于蛋白质）和图CNN（用于配体）的组合，它们被插入到一个多层感知器（MLP）中；输入配体表示为SMILES，蛋白质表示为序列。RF-Score-VS（32）是一个随机森林模型，它依赖于分子间原子-原子接触的发生次数来预测结合亲和力。上述大多数方法都是在PDBbind数据库的子集上训练和测试的。（23）最后，GNINA（6）和Boltz-2（33）与其他方法有所不同：GNINA的亲和力模型（结合CNN和基于物理的函数）嵌入在对接过程中，因此最终选择的配体姿势旨在最小化预测的BFE值，而Boltz-2是一个共折叠模型，它预测结合亲和力以及蛋白质-配体复合物结构。将这七种ML方法应用于632个PDBbind集以预测BFE值的结果总结在图3中（所有计算值作为补充数据集6的一部分提供）。两种方法，KDEEP和OnionNet-2，显示出与其原始出版物中报告的准确性相当。相比之下，另外两种方法TopologyNet和Yuel的性能明显较低，皮尔逊相关性略高于0.5，而原始来源报告的值约为0.75。Boltz-2的性能介于两者之间，皮尔逊r和斯皮尔曼ρ均为0.58。最后，GNINA和RF-Score-VS在与实验数据的相关性方面表现最弱。GNINA的皮尔逊r = 0.34（斯皮尔曼ρ = 0.52），而RF-Score-VS的皮尔逊r = 0.30（斯皮尔曼ρ = 0.35）。目标内相关性与总体趋势一致（详情见补充表S2）。

图3
图3. 实验BFE（ΔGexp）与ML模型之间的相关性：（A）KDEEP（ΔGKDEEP），（B）OnionNet-2（ΔGOnionNet–2），（C）Yuel（ΔGYuel），（D）TopologyNet（ΔGTopologyNet），（E）GNINA（ΔGGnina），（F）Boltz-2，以及（G）RF-ScoreVS。

在高分辨率图像
下载MS PowerPoint幻灯片

DUD-E基准性能指标
在建立了一组具有已知3D结构和BFE值的配体-蛋白质对的准确性基准后，我们评估了基准方法在区分DUD-E数据集中的结合剂和诱饵方面的表现。可以使用标准分类计数来计算灵敏度、特异性和ROC-AUC（54,55）等定量指标。然而，我们的基准数据集规模较小（每个目标15个活性分子和20个诱饵）以及活性配体的比例人为地偏高，使得这些标准指标不太能代表现实世界的虚拟筛选环境。在实际场景中，实验室通常会从更大的对接生成的命中列表中选择100个配体进行实验测试。根据CACHE挑战1和2中报告的命中率（12,56），初次对接筛选的真阳性率约为1%。因此，从100个化合物的选择中识别出真正的活性分子通常需要重新评估至少1000个对接命中，其中大约包含10个真正的活性分子。为了使我们的性能指标在这个实际背景下具有可解释性，我们出于两个特定目的使用了超几何分布。首先，我们用它来建立一个定量统计估计，以比较观察到的真阳性率与随机化合物选择之间的差异。它提供了在从包含恰好Npos个活性分子的Ntot个化合物总集中随机选择n个化合物时成功识别k个真阳性的确切概率：??(??)=(??????????)(??????????????????????)(??????????)P(k)=(Nposk)(Ntot?Nposn?k)(Ntotn)(1)在我们的基准设置中，Ntot = 35且每个目标Npos = 15，我们选择了n = 7个配体来评估成功率。提供显著预测价值的方法显示出的真阳性率在这个确切分布下几乎不可能通过随机机会观察到。其次，我们使用超几何分布将这些基准结果数学地投影到更现实的虚拟筛选场景Ntot = 1000、Npos = 10和n = 100上。通过将我们特定的基准设置与这个更大范围的假设分布之间的累积概率进行匹配，我们将有限数据集上的原始性能转化为一个有意义的估计，表明每种方法在实际的大规模筛选活动中的效果如何——在这些活动中，活性化合物极为罕见（详见方法部分）。总体而言，结合剂和诱饵的BFE值概率密度函数（PDFs）在图4中显示出明显不同的分布，尽管存在一定的重叠（所有基础数据都包含在补充数据集7中）。在实际的虚拟筛选（VS）环境中，由于诱饵的数量远多于活性化合物的数量，这种重叠可能会导致比我们的基准研究中观察到的更高的假阳性率。在这个基准测试中，59个预测为阳性的样本实际上是真实的阳性样本，对应的True Positive Rate（TPR）为94%。根据方程1计算，在随机选择的配体列表中观察到这种TPR的概率接近于零（P(k) ≈ 10^-20）。将这些统计数据外推到实际环境中（Ntot = 1,000，Npos = 10，n = 100），并映射相应的超几何分布，表明所有10个活性化合物都能被识别出来。此外，即使只有17个阳性样本被识别出来，也能捕捉到全部10个活性化合物。这一趋势在所有九个目标中都是一致的（见补充表S3），其中五个目标的TPR为100%，其余四个目标的TPR为86%。

图4 显示了活性化合物（红色）和诱饵（黑色）的计算BFE值的概率密度函数：(A) alchemical ABFE；(B) MM/GBSA；(C) MM/PBSA。

使用MM/GBSA和MM/PBSA计算的活性化合物和诱饵的结合自由能（BFE）值的PDFs在图4中显示出明显的分离（所有基础数据都包含在补充数据集7中）。MM/GBSA和MM/PBSA的总体TPR分别为70%和83%。根据方程1计算，在随机选择中观察到这些TPR值的概率极低（P(k) < 10^-6）。将MM/GBSA的性能外推到实际环境中（Ntot = 1,000，Npos = 10，n = 100），表明可以识别出8个活性化合物。对于MM/PBSA，在100个配体的命中列表中可以识别出所有10个活性化合物。值得注意的是，即使在只有38个配体的命中列表中，MM/PBSA也能识别出全部10个活性化合物。然而，这种总体高性能在九个目标中并不一致（支持信息表S3）。对于MM/GBSA，四个目标的TPR为57%，与随机选择的概率相当。对于HMG-CoA还原酶（HMDH）这一个目标，MM/GBSA的表现明显低于随机选择器，TPR仅为14%。MM/GBSA的总体高性能主要归功于它在这四个目标上的出色表现：FAK1（TPR = 100%），MCR（TPR = 86%），PRGR（TPR = 100%）和TRYB1（TPR = 100%）。相比之下，MM/PBSA在九个目标上的表现更为一致（支持信息表S3）。它仅在HMDH这一个目标上的表现明显低于随机选择器，TPR为29%，与随机选择器的结果相当。对于其余目标，MM/PBSA的TPR在71%到100%之间。

MM/PBSA相对于MM/GBSA的优越性能与理论预期一致。MM/PBSA通过数值方法解决了Poisson-Boltzmann方程，提供了更严格的电静力环境的物理表示。相比之下，MM/GBSA依赖于广义Born近似，这种方法计算速度更快，但简化了极性溶剂化自由能的处理，因此通常会导致较不准确的结合亲和力预测。

由于135个活性配体-蛋白复合物中有55个的X射线结构是实验确定的，我们进行了回顾性分析，研究结构数据的来源如何影响基于物理的方法的性能。在这55个配体中，43个的对接姿势与其相应的X射线结构非常相似。对于剩下的12个配体，对接姿势的偏差较大，与实验坐标的均方根偏差超过3 ?。使用这些对接姿势作为分子动力学的起始几何结构，而不是X射线结构，MM/GBSA的整体性能有所提高，而alchemical方法的性能略有下降；MM/PBSA的性能则保持不变。具体来说，对于MM/GBSA，在两个目标（ADRB1和GRIK1）上的性能从随机选择的基线（NTP = 4）提高到了合理的富集水平（NTP = 5）。相反，alchemical方法在PGH2上的命中数减少了一个（NTP为5而不是6）。在从不同对接姿势开始的12个alchemical计算中，有五个的计算显示出显著差异（大于2 kcal/mol），与从X射线结构开始的计算相比。这五种计算中有四种预测的亲和力较弱，从而增加了假阴性的风险。这种行为是可以预期的，因为天然的实验姿势应该对应于最低的结合自由能。尽管如此，我们的结果表明，即使这些被低估的亲和力仍然足以进行正确的分类，因此对总体真阳性率没有显著影响（详见支持信息数据集7中的BFE值）。

视觉上，五种机器学习（ML）模型计算的活性化合物和诱饵的结合自由能（BFE）值的概率密度函数（PDFs）几乎完全重叠（所有基础数据都包含在补充数据集7中）。这五种早期模型的总体TPR在20%到38%之间。在随机选择中观察到TopologyNet的20%的TPR的概率为11%，这是Ntot = 307，Npos = 135，n = 63的超几何分布所能达到的最高概率。观察到KDEEP的24%的TPR的概率也很高（约为2%）。在这些早期方法中，只有GNINA的总体TPR为38%，与随机选择器的概率（P(k) = 9.0 × 10^-4）显著不同。相比之下，两个较新开发的模型Boltz-2和RF-Score-VS分别取得了81%和79%的总体TPR，它们通过随机选择发生的概率几乎为零（P(k) = 2.1 × 10^-11和1.6 × 10^-10）。

相应地，在实际的环境中（Ntot = 1,000，Npos = 10，n = 100），四种早期的ML方法（KDEEP、Yuel、OnionNet-2和TopologyNet）在100个配体的命中列表中只能识别出1到3个活性化合物，其富集程度与随机选择相当。它们在九个目标上的表现一致较差（补充表S4），每种方法生成的命中列表中只有三个比随机选择更加富集的活性化合物。GNINA在实际的VS环境中识别出了10个中的5个活性化合物，在其他早期ML方法中表现突出。然而，由于它在九个目标中有五个目标上的表现较差，将其作为初步筛选工具是有风险的。相比之下，将较新的Boltz-2和RF-ScoreVS模型的性能外推到实际环境中表明，这两种方法都能成功识别出所有10个活性化合物。此外，与早期的ML方法不同，Boltz-2和RF-ScoreVS在九个目标上的表现都相当一致，仅在两个目标上表现较差——Boltz-2在HMDH上表现较差，RF-Score-VS在GRIK1和HMDH上表现较差。

这项工作填补了基于分子动力学（MD）方法的虚拟筛选关键评估中的一个空白。此外，我们还在一个标准化的基准协议中比较了计算密集型、基于物理的方法和机器学习技术。虽然之前已经讨论过在VS背景下应用alchemical和end-state方法，但迄今为止还没有大规模的评估研究。先前的大规模研究主要关注相对结合自由能（RBFE）计算，以支持化学序列中的先导优化；在这些研究中使用了商业软件FEP+。相比之下，ABFE评估目前还仅限于较小规模的研究。只有一项研究在类似VS的环境中专门测试了alchemical ABFE方法的性能，使用了少量的活性化合物和诱饵针对三个目标。然而，从包含等量活性化合物和诱饵的数据集中得出的ROC-AUC值并不能准确反映实际VS活动中的性能。有趣的是，这项研究还发现实验BFE值与计算出的ABFE值之间没有相关性，可能是由于使用了较短的λ窗口（1-3 ns）。

我们设计了这项研究来回答几个关键问题：alchemical ABFE计算与涉及多种蛋白质目标和配体的复合物的实验BFE之间的相关性如何？alchemical ABFE在区分结合剂和非结合剂方面的有效性如何？在实际的VS环境中，当非结合剂的数量远多于结合剂时，其性能如何？它与更经济实惠的BFE预测方法（包括基于物理的方法和ML模型）相比如何？为了确保可访问性和可重复性，我们使用了公开可用的软件（GROMACS）和标准的ABFE协议，而不是大多数先前大规模研究中使用的商业软件FEP+。尽管ABFE计算仍然是一个活跃的研究领域，但我们特意选择了易于非专家采用的标准化实现。具体来说，我们确定了一组在广泛的蛋白质目标上表现良好的参数。虽然采用可变长度的λ窗口可能会节省一些计算资源，但我们优先选择了更为保守的协议，以实现自动化、无需人工 curate 的执行。

选择基准数据集是任何回顾性虚拟筛选评估中的关键因素。特别是，我们用于类似VS基准的DUD-E数据集存在已知的局限性和偏差，例如诱饵有时可以根据简单的物理性质而非真实的结合相互作用来区分，以及给定目标下活性化合物的化学多样性有限。虽然开发了新的数据集（如MUV、DEKOIS 2.0和LIT-PCBA）来解决这些问题，但每个数据集都有自己的权衡。例如，MUV和DEKOIS 2.0覆盖的蛋白质目标范围较窄，而基于实验高通量筛选数据的LIT-PCBA可能会引入干扰和实验上的假阳性或假阴性。尽管存在这些不完善之处，DUD-E仍然是一个非常实用的选择，并且在该领域被视为事实上的黄金标准。使用DUD-E确保我们的性能指标可以直接与其他现有的虚拟筛选方法文献进行比较和定位。

我们的研究结果直观地根据计算成本对BFE预测器进行了排名。alchemical ABFE方法在定量预测和区分对接的结合剂和非结合剂方面的表现都更为出色。其中一个最意外的发现是end-state MD方法和ML模型之间的角色逆转。虽然end-state方法在PDBbind数据集上没有显示出计算BFE值与实验BFE值之间的相关性，但它们在DUD-E数据集中识别活性化合物的表现相对较好。相比之下，ML模型在PDBbind基准测试中的表现优于end-state方法，可能是因为这部分数据与它们的训练数据有重叠。有趣的是，两种ML模型Yuel和TopologyNet在我们的PDBbind基准测试中的表现不如在它们训练用的更广泛的PDBbind数据集中的表现，这可能是由于包含了噪声较大的IC50数据作为结合亲和力的代理导致的过拟合。GNINA是一个例外，它不像大多数ML模型那样表现，而是在定量PDBbind数据上的表现不如其他ML模型，但在其他方面表现得更好，表明它具有更高的泛化能力。虽然大多数ML模型在DUD-E数据集上的表现较差，生成的命中列表并不比随机选择更好，但GNINA的表现与MM/GBSA相当。总体而言，经过基准测试的ML模型在超出其训练范围的化合物上表现一致较差。像许多其他模型一样，这些模型是在PDBbind数据集上训练的，该数据集包含了几千个实验验证具有结合能力的蛋白质-配体复合物。由于在训练过程中它们从未接触到真正的非结合剂，因此它们在DUD-E数据集中无法区分结合剂和诱饵也就不足为奇了。在之前的研究中，领域适用性问题大多被忽视了，只有最近的一项研究（73）明确调查了它们对分布外（OOD）化合物的性能。然而，即使是那项研究也仅限于PDBbind数据集，并没有完全解决泛化到真实非结合剂的挑战。除了整体性能指标外，我们还分析了基准方法在DUD-E基准的九个目标上是否表现出目标特定的偏差（见补充表S2、S3）。分析显示不同方法之间的相关性非常弱，这表明它们的失败点主要是方法特定的。例如，基于最终状态的物理方法和Boltz-2在目标HMDH上表现不佳，真正的阳性结果（NTP）的数量处于或低于随机水平（≤4），而一些较旧的机器学习模型，如Yuel和OnionNet-2，对这个目标实现了真正的富集（见补充表S3）。相比之下，目标GRIK1对大多数较旧的机器学习方法来说极具挑战性，它们未能超过随机基线性能，但Boltz-2和一种炼金术方法成功得分。此外，目标PGH2对几乎所有机器学习方法都是一个重大挑战，只有GNINA和RF-ScoreVS实现了有意义的富集。炼金术结合自由能计算是唯一在整个目标面板中保持持续高表现的方法（见补充表S2）。最后，我们还有机会进行一个小规模的分析，研究在炼金术模拟中使用对接提出的起始配体几何结构对VS结果的影响。虽然使用对接是不可避免的，但它对整体真实阳性率的负面影响似乎很小。在DUD-E基准中的55个有活性配体子集中，大约80%的对接姿势与天然几何结构非常吻合——这一高成功率得到了广泛的社区对接基准的证实（74,75）。在12个不同的姿势中，有四个导致了预测亲和力的减弱，代表了潜在的假阴性，但实际上只有一个真正实现了这种风险。因此，基于这一有限的评估，人们可能会预期在典型的筛选活动中，由于错误的对接姿势，有2%到5%的有活性化合物（如果有的话）会被误判为假阴性。

除了预测准确性之外，任何VS分诊方法的实际效用在很大程度上取决于其计算成本和计算基础设施的可用性。直接比较这些成本往往具有挑战性，因为计算时间只是总体费用的一个维度，硬件要求和并行化能力同样至关重要。在我们的研究中，评估的方法通常可以分为三个不同的计算吞吐量层次：资源最密集的是基于物理的方法，每个配体的处理时间以小时计。完整的炼金术计算在我们的本地Hellbender集群上的专用高性能GPU上大约需要2小时，由于GPU的可用性通常有限，同时运行多个模拟是困难的。MM/PBSA和MM/GBSA每个BFE计算也需要大约2小时，但它们消耗的计算资源大约是炼金术方法的50分之一，允许在可用节点上进行广泛的并行执行。第二层包括对接/共折叠机器学习模型，如GNINA和Boltz-2，它们可以在几分钟内处理复合物。最后，最快的层次是机器学习模型，每个复合物只需要一次前向传递，并且只需要一小部分秒时间。由于这些标准机器学习计算所需的硬件资源非常少，它们可以大规模并行化，以实现极高的筛选吞吐量。

鉴于本研究的结果，可以自信地推荐使用炼金术ABFE方法来优先处理虚拟筛选中的化合物，因为它在所有目标上表现出一致的高性能，真实阳性率在86%到100%之间。然而，其显著的计算成本可能使其难以处理大型化合物库。在评估的其他方法中，AI模型Boltz-2作为一个高效的选择脱颖而出，其表现与基于最终状态的物理方法相当，但不同目标之间的一致性非常好，并且计算成本只有后者的一小部分。如果必须使用基于最终状态的方法，则MM/PBSA比MM/GBSA更受青睐。MM/GBSA在九个目标中有四个表现不佳，突显了实践者可能遇到的高变异性，并引发了对其可靠性的担忧。为了应对这种变异性和计算成本，我们建议采用多阶段的工作流程。实践者可以使用像Boltz-2这样的高效AI方法来预筛选大型化合物库（大约10^5个配体），然后对排名前几百个的命中进行严格的炼金术ABFE计算，以确保高精度地识别有活性化合物。我们希望这项研究能够鼓励在未来将基于MD的方法实际应用于具有挑战性的治疗目标的候选物发现。随着计算能力的提高和负担得起的基于云的解决方案的日益普及，这些方法预计将变得更容易获取。虽然机器学习模型尚未准备好广泛部署，但我们相信学习策略的进步，包括在PDBbind之外的数据集上进行训练，以及将配体-蛋白质相互作用整合到神经网络中的新方法，将推动未来的发展。此外，本研究生成的数据，包括MD轨迹，也可能有助于改进基于物理和机器学习的方法。这项工作还可以激发专门针对MD衍生数据训练的学习模型的开发，以预测结合亲和力。

**方法**
**PDBbind基准集**
使用一系列筛选标准评估了4,284个独特的配体，这些配体对应于PDBbind（23）中的精炼子集，包含5,316个蛋白质-配体复合物。数据于2022年10月26日从PDBbind v2020下载。筛选标准包括：1）改进的Lipinski规则（36）（分子量在250到600 Da之间，cLogP在1到5之间，氢键供体≤5，氢键受体≤10，可旋转键<14）；2）结构筛选（无磷原子，不超过两个可电离基团，至少有一个环，无长脂肪链，无大型融合环系统，无潜在有毒或高反应性基团）。所有筛选都使用RDKit库在Python中实现（76）。此外，只保留了活性数据以Kd/Ki表示并且关系为“=”的配体。总体而言，筛选后将数据集减少到1,286个蛋白质-配体复合物。相应的蛋白质结构从RCSB PDB下载（77），并从蛋白质结构中分离出配体。自动蛋白质准备（见下文）在自动化处理步骤中未能处理147个PDB文件（如Open Babel配体格式转换、蛋白质链分割、不完整的3D结构等），最终得到1,139个蛋白质-配体复合物。此外，由于标准力场的覆盖不足，移除了285个金属蛋白。一些系统因MD运行时的错误（如范围检查违规或内存相关问题）而失败。该集合后来进一步减少到632个复合物（或606个独特配体）。

**DUD-E基准集**
配体结构以SMILES格式从DUD-E网站下载（https://dude.docking.org/）。对于诱饵集，每个目标随机选择了200个化合物，并以“灵活”模式进行了Glide对接计算。然后根据其排名的Glide得分选择了每个目标的top 20个化合物。活性集的配体选择过程与诱饵相同；基于对接得分选择了135个化合物（每个目标15个）。135个活性化合物中有55个的X射线结构可用，并被用作MD模拟的起始坐标。

**数据准备**
**配体结构**
从PDB文件中提取结构，分配键顺序，并使用RDKit库的Python代码添加氢原子。得到的结构被保存为SD文件。
**蛋白质准备**
使用Schr?dinger Suite v2022–4（78）的工作流程处理蛋白质结构。该工作流程添加了氢原子，重新分配了键顺序，并生成了缺失的残基。蛋白质的质子化状态在pH 7.0下使用PROPKA进行分配。

**分子动力学数据**
使用Bash和Python脚本来准备MD输入文件。使用AMBER14SB力场（79）为蛋白质生成拓扑和参数文件，使用GAFF力场（80）为配体生成参数文件。将配体和蛋白质结构组合成单一的配体-蛋白质复合物，格式为GRO。使用bash脚本以批处理模式进行能量最小化和平衡MD运行。通过监控均方根偏差（RMSD）来评估MD模拟期间配体-蛋白质复合物的完整性。接下来，对通过完整性检查的复合物进行了MM/GBSA、MM/PBSA和FEP计算。最终，成功完成了632个复合物的MM/GBSA、MM/PBSA和FEP计算。

**对接**
我们使用AutoDock Vina 1.2.0（14）进行分子对接研究。目标蛋白质使用AutoDock Tools版本1.5.6准备为PDBQT文件。网格框在x、y和z维度上设置为40 × 40 × 40点，网格间距为0.375 ?。网格框的中心位于相应的X射线复合物中的配体处，以生成网格定义文件（GDF文件）。最初的SMILES格式配体结构使用Open Babel转换为PDB格式的三维几何结构。使用MGLTools的prepare_ligand4.py脚本向得到的结构中添加氢原子，并进一步保存为AutoDock对接所需的PDBQT文件。使用了默认的8个 Exhaustiveness参数值。

**分子动力学**
使用GROMACS v2022.2软件进行分子动力学（MD）模拟（79）。应用AMBER14SB力场（81）进行蛋白质能量计算，而配体参数化使用GAFF力场（80）。使用高斯16（82）在HF/6-31G*水平上计算ESP电荷，得到限制性静电势（RESP）电荷参数。使用ANTECHAMBER模块的AMBER18程序通过Python脚本acpype.py（83）准备GROMACS格式的配体拓扑。每个配体-蛋白质复合物放置在一个立方模拟盒中，复合物与盒子边界之间的最小距离为10 ?。模拟在显式溶剂中进行，使用TIP3P水模型，并用0.15 M NaCl中和。系统使用最速下降算法进行能量最小化，以解决空间冲突和不利的原子相互作用，当最大力达到2.4 kcal/mol/nm时终止。随后在298 K和1 atm下进行了100 ps的NVT和NPT系综平衡运行，应用239 kcal/mol/?2的重原子约束以稳定温度和压力。生产阶段包括5 ns的无约束MD模拟，使用速度重标度恒温器（84）。时间步长为2 fs，非键合相互作用的截止值为12 ?。使用Particle Mesh Ewald（PME）方法计算长程静电相互作用（85）。

**MM/GBSA和MM/PBSA计算**
MM/GBSA和MM/PBSA都使用以下方程进行BFE计算，参数设置如参考文献（86）所述：
Δ????????=Δ??????+Δ??????????+Δ????????‐?????????????Δ??
ΔGbind=ΔEMM+ΔGpolar+ΔGnon‐polar?TΔS
其中ΔEMM是由分子力学（MM）力场计算的气相系统能量，ΔGpolar是由Poisson–Boltzmann或Generalized Born模型计算的极性溶剂化自由能，ΔGnon–polar是由加权溶剂可访问表面积（SASA）估计的非极性溶剂化能量，TΔS是熵贡献。使用gmx_MMPBSA v1.6.2工具（87,88）进行了两种替代方法，即正常模式分析（50,51）和IE（52）来计算熵贡献。BFE计算使用了MD轨迹的最后1 ns（1,000个快照）。由于计算成本高，正常模式分析仅限于20帧。IE计算涉及所有1,000帧。

**炼金术自由能计算**
我们的工作流程的核心部分改编自GROMACS教程（89），进一步自动化以适合高通量使用，并针对密苏里大学的Hellbender集群进行了定制。在模拟之前，为两种耦合路径（配体在水中的耦合和配体与溶剂中的蛋白质耦合）生成了一组25个λ窗口。λ-窗口的范围从 λ = 0.0 到 1.0，其中 λ = 0.8 之前的间隔为 0.1，而 λ = 0.85 到 1.0 之间的间隔为 0.05。在前 12 个窗口中，逐渐激活范德华（vdW）势能，同时关闭库仑相互作用。在剩余的窗口中，逐步引入库仑相互作用。使用 Bennett 接受比率（BAR）方法（45）来计算状态之间的 ΔG。为了防止模拟过程中粒子重叠，定义了软核参数设置：sc–alpha = 0.5，sc–power = 1，和 sc–sigma = 0.3（90）。ΔH/λ 的值以 100 fs 的间隔记录。完整的 BFE 热力学循环包含 50 个 λ-窗口，每个窗口的分子动力学（MD）模拟时间为 5 ns。取向约束是使用 Boresch–Karplus 算法实现的（91）。通过使用 OpenBabel 库去除配体并添加氢原子来准备 Machine Learning BFE Predictors KDEEPProteins；处理后的结构保存为 PDB 格式。配体的制备方法是生成它们在预定义箱子内的三维坐标并添加氢原子；处理后的配体结构保存为 SD 格式。准备好的输入文件被上传到 OpenPlayMolecule.org，在那里于 2024 年 11 月 16 日在线运行 KDEEP 以计算 BFE 值。

对于 OnionNet-2，蛋白质和配体的制备过程与 KDEEP 相同。OnionNet-2 环境的搭建遵循 OnionNet GitHub 仓库中提供的安装说明。OnionNet-2 代码从 github.com/zchwang/OnionNet-2 下载，并使用准备好的输入文件于 2024 年 11 月 18 日运行以预测 BFE 值。

对于 TopologyNet，蛋白质和配体的制备过程也与上述两种方法相同。TopologyNet 环境的搭建遵循 TopologyNet GitHub 仓库中提供的安装说明。环境配置完成后，使用准备好的输入文件运行 TopologyNet 代码（2024 年 11 月 20 日）以预测 BFE 值。YuelProteins 以纯文本格式的一字母序列形式准备，配体以 SMILES 字符串形式准备，然后将它们保存在 Excel 文件中。Yuel 环境的搭建遵循相应的安装说明。代码和预训练模型从 bitbucket.org/dokhlab/yuel/src/master 下载，并于 2024 年 11 月 22 日运行以计算 BFE 值。

GNINAProteins 是通过使用 OpenBabel 库去除目标链来准备的；处理后的结构保存为 PDB 格式。未修改的配体结构用于定义对接箱坐标；配体没有进行进一步的修改。GNINA 环境是使用 github.com/gnina/gnina 上提供的 Google Colab 笔记本搭建的。GNINA 搭配在 2025 年 1 月 3 日运行，使用其默认评分功能来预测 BFE 值。

为了评估 BFE 预测器在虚拟筛选（VS）中的性能并将其结果外推到更大的数据集，我们使用了 Python 的 scipy.stats.hypergeom 模块中的超几何分布。进行了两项特定的分析。首先，使用超几何分布的概率质量函数（PMF）计算实验数据集中恰好观察到 kobs 个真实活性物的概率（P(kobs)）。PMF 的计算公式为 hypergeom.pmf(k,N,K,n)，其中 N 是总样本数量，K 是真实活性物的总数，n 是样本大小，kobs 是观察到的真实活性物数量。这项计算确定了观察到的真实活性物数量是否可以通过随机机会合理地产生，从而对富集度进行了定量评估。其次，为了推断更大数据集中预期的真实活性物数量（k2），将实验设置中观察到的累积概率 P(k≥k2) 与放大后的场景中的相应累积概率进行匹配。给定数据集的累积概率 P(k≥k2) 使用 1–hypergeom.cdf(k–1,N,k,n) 计算，其中 hypergeom.cdf 表示超几何分布的累积分布函数。对于放大的数据集（实际世界的 VS 场景），参数为 k2,N2,K2,n2，我们迭代求解 k2，使得累积概率 P(k≥k2) 与实验设置中的 P(k≥kobs) 紧密匹配。确定了满足这一条件的最小 k2 值。