引言
在药物发现领域，预测化学修饰如何影响配体与靶蛋白的结合能力是理性药物设计的核心。相对结合自由能估算方法，例如自由能微扰，已成为评估配体结合亲和力变化的宝贵工具。这些基于物理原理的计算方法能够帮助药物化学家优先选择待合成的化合物，从而优化先导化合物优化流程。在参数设置得当的情况下，FEP方法的预测误差通常能控制在1 kcal/mol以内，这为工业界和学术界的基于模拟的化合物优先排序提供了可能。
尽管计算能力和力场质量不断提升，FEP的广泛采用仍受限于一些实际挑战：方法通常需要专家知识进行系统设置、耗费大量计算资源并进行细致的质量控制。诸如薛定谔公司的FEP+等商业平台提供了成熟、用户友好的解决方案，但在可访问性和定制化方面存在限制。包括OpenFE和pmx在内的开源框架则提供了具有竞争力的透明替代方案。然而，在自动化与计算效率方面，尤其是对于大规模的扰动网络或显著的化学结构变化，挑战依然存在。
研究团队此前已介绍了QligFEP，这是一个在Q分子动力学引擎中使用球形边界条件的模块化工作流程。该方法将模拟限制在结合位点周围的聚焦区域，减小了系统规模，并支持基于CPU的执行，同时保留了计算出的配体-蛋白质相互作用力的准确性。本文介绍的QligFEP v2.1.0版本，通过专用的命令行工具并在需要时辅以图形界面，在工作流执行、灵活性和可用性方面进行了重大改进。更新后的平台具备自动化的配体映射与网络生成、简化的输入准备、可配置的距离限制算法，以及带有稳健错误处理能力的增强型分析工作流。这些进展减少了人工干预，支持脚本化、可重复的执行，并使QligFEP适用于高通量应用。
方法
基准数据与准备
基准数据来源于IndustryBenchmarks2024存储库，包含JACS和Merck基准集。输入结构从指定的代码提交版本中获取，除了pfkfb3系统。大多数系统使用了提供的“准备结构”，而ptp1b和thrombin系统则因使用存储库提供的结构时模拟崩溃次数异常增多以及与实验值的相关性较差而进行了独立的蛋白准备。独立的蛋白准备从相同的源PDB标识符开始，进行了基于pKa的质子化状态分配和局部能量最小化。
结构精修
高质量的输入结构对于FEP计算至关重要。通过微小的调整解决了空间位阻冲突并使结构符合力场惯例。结构操作，包括视觉检查和能量最小化，使用了薛定谔公司的Maestro软件。选定的结合位点残基使用基于OPLS力场的集成工具进行了局部最小化。对于某些系统，基于晶体学证据调整了两个配体的构象，以改善比对并提升与实验数据的相关性。
QligFEP输入准备
扰动网络是使用改编自IndustryBenchmarks2024的脚本生成的。Kartograf基于笛卡尔坐标重叠，在同源配体之间映射原子，将转化限制在相同电荷的配体对之间。过滤功能确保了合理的转化选择。LOMAP评分用于评估扰动，生成具有默认布局的网络。cdk8系统需要手动映射以避免自动映射算法提出的次优扰动。
从SDF配体结构和PDB蛋白文件生成输入文件需要最少的用户干预。qparams工具生成Q兼容的配体库和参数文件。命令行程序qprep_prot创建以配体几何中心为球心、半径为25 Å、包含结合位点的SBC系统。球体被水合，并保持与重原子至少3.0 Å的最小距离。在单个边创建期间，会移除与特定转化中配体发生冲突的水分子。默认设置包括10个重复以及通过RestraintSetter算法自动放置限制。
力场与模拟协议
模拟对蛋白质使用Amber ff14SB力场，对水使用TIP3P模型，对配体使用OpenFF Sage 2.2.1力场并搭配AM1-BCC电荷。
模拟使用25 Å的球形液滴并采用SCAAS表面约束模型。边界水分子受到限制；球体外的蛋白质和配体原子受到谐波约束。10 Å以外的静电相互作用使用局部反应场方法近似处理，参与炼金术转化的原子除外。
相对结合自由能的计算采用Q的双拓扑方法。等效重原子之间的距离限制用于保持相空间重叠，当距离超过0.1 Å时激活。扰动跨越101个S型分布的λ窗口，从λ=0.5开始，分别独立地向λ=0和λ=1进行。生产运行对每个λ窗口采样10 ps。每条分支（蛋白质、水）进行10次独立的重复运行，使用唯一的随机种子以确保可重复性——每条边总共20次重复并行运行。
系统规模平均而言，蛋白质分支包含6069 ± 92个原子，水分支包含6359 ± 9个原子。每个重复的总模拟时间为1141 ps。在CPU集群上的典型执行时间为1小时36分钟。
限制设置算法
QligFEP v2.1.0引入了RestraintSetter算法，它基于不同级别的原子和取代基等效性，为双拓扑FEP自动定义位置限制。算法首先使用Kartograf基于笛卡尔坐标的原子对应关系来识别和映射分子间对应的环状结构。接着，算法使用用户指定的标准分层比较分子结构，例如元素、杂化状态、芳香性或所有重原子等效。算法还通过三个严格级别考虑环状结构周围的化学环境，用于环取代基的比较：宽松型、较不严格型和严格型。这种多层次方法适应了多样的分子对和研究需求。
分析
qligfep_analyze命令行工具使用Bennett接受率法计算自由能。对于每个λ窗口，前100步被舍弃。相邻λ对的BAR估计值被求和以获得总自由能变化。ΔΔGbind通过热力学循环结合水分支和蛋白质分支的结果计算得出，不确定性从10个重复的标准误差传播而来。
性能评估使用等级相关系数Kendall‘s τ和平均无符号误差。指标包括通过自助法计算的95%置信区间。统计显著性检验使用Mann–Whitney U检验并经过Holm-Bonferroni校正。此外，还提供了一个图形用户界面，用于交互式分析获得的结果、扰动网络以及模拟边所使用的输入结构。
绝对结合自由能使用基于状态函数的校正算法进行循环闭合校正。该方法利用自由能作为状态函数的特性，即围绕封闭热力学循环的净自由能变化必须为零。SFC算法为网络中的每个配体分配一个绝对结合自由能，并全局优化这些值以最小化与计算出的成对ΔΔGbind值之间的差异。
对比分析
将QligFEP在16个靶标上的性能与已报道的pmx/GROMACS搭配Amber ff99SB-ILDN和Sage 2.0力场的结果，以及薛定谔FEP+搭配OPLS3/OPLS3e力场的结果进行比较。本文讨论的主要比较使用未经后处理的原始ΔΔGbind值，以直接评估方法和力场性能。
结果与讨论
QligFEP协议的计算效率
QligFEP的一个关键优势是通过SBC实现的高计算效率。通过将模拟聚焦于结合位点区域，该方法在保持预测准确性的同时，显著降低了与全系统方法相比的计算需求。使用tnks2基准系统进行的并行扩展分析表明，在最多16个CPU核心时具有良好的效率。蛋白质分支模拟实现了4.6倍加速，水分支模拟实现了4.4倍加速。每条FEP边的总墙上时间从使用2核心时的105.76小时减少到使用16核心时的23.46小时，使用32核心时进一步减少有限。8个核心提供了最佳的资源效率，具有3.1倍的加速和较高的并行效率。超过8个核心则收益递减。按照当前AWS竞价实例定价，计算一条FEP边的成本低于1美元，这使得大规模的FEP活动在经济上可行。研究推荐使用8个核心，以平衡速度与成本效率。
SBC方法模拟的原子数远少于采用周期性边界条件的全系统模拟。尽管系统规模较小，模拟仍能捕捉到必要的蛋白质-配体相互作用。水分支通常比蛋白质分支需要更少的时间。工作负载分布在三种HPC架构上表现出了一致的性能，证明了其可移植性。计算成本因靶标而异，所有靶标的平均计算成本为每条边29.79 ± 3.79小时的墙上时间。
与成熟方法相比的准确性
在确立了计算效率后，我们评估了其相对于成熟FEP方法的准确性。对所有16个靶标的绝对ΔΔG预测误差进行的统计比较显示，QligFEP与PMX-Sage 2.0相比没有显著更大的偏差。与FEP+ OPLS3e相比，QligFEP仅在shp2靶标上显示出显著更大的偏差，这表明SBC模拟没有引入系统性误差，并且与商业和开源协议相比仍具竞争力。对于pfkfb3靶标，尽管存在离群值，QligFEP的偏差与FEP+ OPLS3e相比未达到显著性水平，而PMX-Sage 2.0则达到了，这表明QligFEP的性能介于两种方法之间。
各系统间的排序性能和准确性比较揭示了所有方法都存在靶标依赖性的差异。一些靶标表现持续强劲，而另一些则对所有方法都构成挑战。重要的是，这些变异在所有三种方法中都是可比的。
QligFEP在6/16的靶标上表现出优于PMX-Sage 2.0的排序能力，特别是在shp2、ptp1b和syk靶标上有显著改善。相比之下，QligFEP在eg5、jnk1、cmet、cdk8和pfkfb3等靶标上表现出较低的性能。与FEP+ OPLS3e相比，这种商业方法在大多数靶标上显示出优势。QligFEP在bace、shp2、syk和thrombin等靶标上表现出相当或略好的性能。
QligFEP的MUE值在几个靶标上普遍高于其他方法。FEP+ OPLS3e在hif2a、pfkfb3和shp2等靶标上显示出明显的MUE优势，这可能反映了其专有力场的益处。总体而言，尽管在某些靶标上MUE值较高，但QligFEP保持了与其他方法相当的排序性能。这对于先导化合物优化应用尤为重要，因为准确地对同源化合物进行排序通常比最小化绝对预测误差更重要。结合其较低的计算成本，QligFEP成为了大规模FEP实验的一个实用选择。
处理复杂的化学转化
QligFEP双拓扑实现的一个关键能力是处理显著的结构变化。图5展示了一个复杂的syk转化，涉及配体CHEMBL3265020和CHEMBL3264999。尽管取代基在大小和柔韧性上截然不同，但计算获得了与实验值合理的一致。RestraintSetter算法正确识别出吲哚和苯基部分应保持无限制，以避免强加不自然的相互作用。这个例子展示了QligFEP处理对单拓扑方法具有挑战性的转化的能力。
自动化的RestraintSetter提供了灵活的配置：宽松型、较不严格型和严格型比较。宽松的重原子比较被证明在大多数网络中通常是合理的。更严格的配置可能有益于具有不同柔韧性和相互作用特征的R基团的支架跃迁。用户应根据其配体系列的结构相似性和自由度考虑限制的适当性，并进行调整。
输入结构质量的重要性
结构质量对FEP性能有至关重要的影响，这在该领域已广为人知。无论输入来源如何，配体比对和网络设计都需要仔细关注，因为不充分的比对会影响单个扰动。此外，当使用Kartograf作为原子映射创建的解决方案时，次优的比对可能会影响整体网络拓扑。当共享较大最大公共子结构的配体表现出次优的结构重叠时，研究应用了基于MCS的比对。在ptp1b系统中识别出4个这样的配体，并优化了它们的比对以优化网络拓扑。对于shp2数据集，在连接SHP099-1-7和配体5的边上观察到了比对问题。尽管化学等价，但由于它们在相对结合方向上的结构对齐不佳，氯取代的苯基缺乏位置限制，导致该边的预测值严重偏高。基于MCS的配体相似性聚类提供了一种系统方法来识别需要精修比对的化合物，从而为FEP计算实现最佳的结构对应。除了研究团队的观察之外，配体比对对于FEP活动成功的重要性也在OpenFE联盟的基准测试中得到了强调。
网络设计的考虑
除了单个扰动的质量，整体网络拓扑也显著影响性能。对于cdk8系统，自动网络生成了具有挑战性的扰动，经常涉及环到非环原子的转化，有时甚至同时扰动两个独立的化学部分。手动重新设计网络带来了实质性改善：Kendall‘s τ从0.28增加到0.40，而平均无符号误差从1.94 kcal/mol降低到1.58 kcal/mol。这表明，深思熟虑的网络设计可以挽救具有挑战性系统的性能。除了QligFEP提供的自动化步骤之外，研究工作流具有灵活性，并为需要用户手动管理的情况提供了足够的工具。未来，将RestraintSetter衍生的相似性值整合到网络生成算法中，可能是改善我们应用中自动网络设计的有效策略，从而减少手动干预的需求。
结论与未来展望
QligFEP提供了一个高效的、基于命令行工作流的开源框架，用于相对结合自由能计算。涵盖16个蛋白质靶标和639个配体转化的全面基准测试表明，QligFEP在保持与开源及商业替代方案相当准确性的同时，所需的计算资源更低。SBC模拟策略将系统规模减少至6214 ± 159个原子，并支持在标准CPU集群上进行高效、可并行化的边计算，在多个蛋白系统中保持一致的运行时间。总之，研究结果支持QligFEP作为依赖GPU的全系统方法的一个可访问替代方案。统计分析证实，这种聚焦模拟方法没有引入系统性误差，在保持预测可靠性的同时提供了卓越的效率。
工作流的可重复性、准确性和计算可访问性的结合，使QligFEP成为加速化合物优化的实用解决方案，特别对计算基础设施有限的研究团队具有重要价值。未来的发展将侧重于通过整合软核势的混合拓扑方法来增强对具有挑战性扰动的稳定性，进一步改进限制处理、网络优化和自动化。此外，支持GPU和云计算的实现目前正在开发中。这些增强将巩固QligFEP作为结构药物设计中一个透明、可扩展平台的角色。
