摘要：本文以D3R CG2的Free Energy Set 2数据集的14个IC₅₀跨度为1000倍的FXR激动剂为例，分析了Flare/FEP预测结合亲和力的性能。结果表明，预测值与实验值具有极其显著的相关性：Kendall Tau=0.64, Pearson R=0.90； MUE=0.35；RMSD=0.35。预测值可以显而易见地将高活性化合物与低活性化合物区分为两类。此外，借助于Cloudam云计算资源，整个计算可以在1个小时左右完成。这说明，借助云计算取之不尽的资源，FEP可以帮助药化/计算团队快速、可靠地对结构相似的同系物进行排序，从而加速项目的进度。

肖高铿/2020-07-27

背景

在药物开发过程中，人们通常希望对活性分子的整体或部分进行替换以改善其特性（活性，选择性，ADME等）或迁移到新的化学空间以获得结构新颖性。在苗头化合物或先导化合物优化项目中，药物化学科研人员通过常规的生物等排体替换、或者类似于SPARK的软件进行基团替换，很容易就获得侧链或骨架替换的一系列化合物结构。问题是：应该优先合成哪个侧链/骨架？本文以活性差异1000倍的几个FXR激动剂同系物来模拟药化科研人员的这个问题，演示FLARE/FEP如何快速地帮你做出决策。

FXR激动剂算例介绍

非酒精性脂肪肝(NASH)是临床未满足治疗的疾病，FXR是NASH的一个热门靶标。在2016年发起的D3R Grand Challenge 2^[1]挑战赛第2阶段的Free Energy Set 2中包含了18个FXR激动剂同系物：这些化合物具有公共的母核、结构相似，但IC₅₀跨度高达1000倍，是非常好的方法学验证数据集。D3R CG2以ΔG实验值与预测值的相关性系数Kendall Tau与以及RMSE为主要指标评价预测结果。

Schindler等人^[2]用Schrodinger/FEP+预测了这18个化合物的相对结合自由能。从结构的净电荷看，其中14个化合物带一个负电荷，另外4个为中性化合物。在该研究中，作者将18个化合物docking到FXR_12与FXR复合物晶体结构中，为了满足FEP对电荷一致的要求，将羧基按中性处理。计算结果如图1所示，FEP+具有良好的排序能力：ΔG预测值与实验值的Kendall Tau =0.5 ± 0.12； Pearson相关性系数R²=0.80 ± 0.08; RMSD = 1.31 kcal/mol。作者还同时用了Schrodinger/Prime MM-GBSA方法进行了ΔG预测，结果表明Pearson相关性系数R²=0.72 ± 0.11; RMSD = 9.31 kcal/mol。作者单独对其中的14个羧酸进行FEP+分析，结果表明RMSD仅降低了0.1。

Figure 1. Schindler等人对18个FXR激动剂的ΔG预测结果

数据准备

用Flare下载FXR_12与FXR的复合物结构(PDB 5Q0W)，并在默认条件下（pH=7.0）准备蛋白结构，并将配体FXR_12从蛋白中提取分离。其它的化合物以FXR_12为模板，在结合位点里用Flare的结构编辑器直接修改、能量最小化后保存。

鉴于FEP方法不能同时处理不同电荷的分子，本文仅以D3R CG2的Free Energy Set 2数据集中的FXR激动剂的14个带负电荷的羧酸分子（见Table 1）为研究对象，考察FLARE/FEP的预测性能。

准备好的蛋白结构：5Q0W_protein.pdb

准备好的配体结构(带FEP预测结果)：fxr-agonist.sdf

准备好用于FEP计算的Flare项目文件(已经包含一个计算结果)：D3R-FXR-FEP.flr

Table 1. 14个FXR激动剂的结构、IC₅₀（nM）及其预测的ΔG (kcal/mol)

FEP计算

详细的操作流程请参见：云计算教程 | 用FLARE FEP计算相对结合自由能^[3]。这里仅对关键的参数进行描述，如果没有说明则为默认参数。

图的生成

FLARE的FEP计算支持两种类型的图：normal与Star图。Normal图模式精度高,也更耗计算资源；Star图适合于研究结构与活性关系研究，但精度有所损失。本文采用以FXR_12为中心的Star图模式（见Figure 2)。实际上，FXR_16比FXR_12更适合做为Star图的中心分子，因为FXR_16与磺酰基连接的末端芳环为苯环，更加容易转化为其它取代苯环或噻唑环，也就更容易计算成功。

Figure 2. 以FXR_12为中心的星型、单向网络

由于涉及到的14个分子IC₅₀已知，本计算实验的主要目的时考察预测值与实验值的关系，因此生成图的时候采用了Benchmark模式，而不是Product模式。

FEP计算参数

FEP计算时小分子力场采用AMBER/GAFF2，并分配AM1-BCC电荷；体系设置10Å的溶剂盒子；每个体系进行了4nM的分子动力学模拟；Figure 2中每对分子之间的link采用单向模式（quick mode）:每个方向上仅计算从小分子->大分子的转化，每个方向默认计算9个λ窗口。鉴于FXR_12->FXR_79有较大的结构变化，因此设置15个λ窗口;FXR_12->FXR_85也设置了15个λ窗口。

硬件资源：云计算加速FEP计算

如Figure 2所示，共有13个Link，总共有111(9*9+2*15)个λ窗口需要计算，因此至少可以使用111个GPU进行加速计算。在本实验中，采用深圳云端软件公司(Cloudam)的云E算力平台的计算资源进行计算，broker设置最大使用200张GPU，每节点1张Tesla T4卡+4个CPU核心+8GB的内存。

具体的云计算资源使用方法见：云计算教程 | 用FLARE FEP计算相对结合自由能^[3]

结果

关于预测性能

ΔG的实验值与预测值见Table 1, Kendall τ与Pearson R相关性系数分别为0.64（p-value=0.001）与0.90（p-value=0.0000145），这表明实验值与预测值之间具有极其显著的相关性；MUE为0.35；RMSD为0.50，比Schindler等人^[2]报道的FEP+预测结果RMSD=1.31要低0.8，具有更好的精密度。

Figure 2. 结合自由能ΔG的实验值与预测值比较

根据Figure 3的ΔG实验值与预测值散点图可以直观地将结合亲和力高、低的化合物聚类区分开：ΔG大于-8.0与小于-8.0的化合物被明显地区分。尤其是IC₅₀大于1000nM的化合物都被被归属于ΔG大的化合物，而IC₅₀最低的几个化合物都聚集在ΔG预测值最小的化合物里。这说明使用Flare/FEP可以帮助药化工作人员可靠地做出化合物排序决策。

更重要的是，几个不同的生物等排体得到很好的排序，尤其是卤素取代苯环的衍生物的排序：1）苯环被噻吩取代对活性不利；2）苯环的氟取代对活性不利；苯环的氯、溴取代对活性有利。这是药物化学研究阶段非常重要的信息，而Flare/FEP可以准确在实验前的告诉你。FEP似乎对低活性化合物识别的特别好：预测为活性最差的5个化合物，确实也对应为IC₅₀最高的5个化合物；在预测活性最高的5个化合物里，包含了实际活性最强的化合物FXR_78。

与不用FEP的方法比较：logP、分子量与打分函数

用Flare计算的SlogP与分子量通常做为空模型，分析SlogP、分子量与实验结合自由能(pIC50)的关系；同时，也用Ligandscout的Binding Affinity Score来预测结合亲和力，分析预测值与实验结合自由能(pIC50)的关系。好的FEP方法，应该要优于空模型与打分函数。结果表明，SlogP、分子量、结合亲和力打分值与pIC50具有强相关性，Pearson相关性系数分别为0.717、0.644与-0.701,远低于Flare FEP的0.9，这说明使用更贵的FEP计算是有价值的，更多的统计学指标见表2。

表2. 结合自由能（pIC50）与logP、分子量、打分函数相关性分析的统计学结果

与pyAutoFEP比较

在2021年的时候，Carvalho等人发表了自动相对结合自由能计算工作流pyAutoFEP^[4]，刚好与本文采用了同样的测试集，除了化合物FXR_10之外。在pyAutoFEP测试中，使用AMBER03/GAFF/AM1-BCC without REST2时性能表现最佳：Pearson Correlation r=0.71；RMSE=1.07 kcal/mol。这个预测精度与不用模型的SlogP差不多，这说明了：在这个例子里，使用pyAutoFEP没有比空模型SlogP更好，也就是用了pyAutoFEP没有比不用更好。

关于机时与费用

本实验在Cloudam上计算完成，Cloudam提供了多种GPU型号可供选择，本次计算实验采用价格最低的Tesla T4进行计算。峰值时Flare调用120多张卡，全部计算实验大约在1.15小时内完成，总共机时费大约600元。这说明，此类FEP计算可以在非常合理地时间内完成，为计算和药化团队提供快速、可靠的结果对化合物进行排序与决策。云计算近乎取之不尽的资源为高精度的FEP计算提供了算力保障，真正地加速了项目研发进度。同时，在本计算中，每个化合物的平均计算费用约为42元，略低于之前^[3]预计的费用。

相关性系数与RMSD计算

Kendall Tau相关性系数计算

1
2
3
4
5

import scipy.stats as stats
exp = [-7.2,-9.9,-8.4,-6.0,-9.0,-10.3,-7.3,-7.6,-9.2,-8.8,-7.3,-8.9,-8.6,-8.4]
pred = [-7.1,-8.5,-9.2,-6.3,-9.0,-10.5,-7.4,-8.0,-8.8,-8.9,-7.5,-8.5,-8.9,-8.1]
tau, p_value = stats.kendalltau(exp, pred)
print(tau,p_value)

1

0.6404494382022472 0.0016945408238685901

Pearson相关性系数R的计算

1
2
3
4
5

import scipy.stats as stats
exp = [-7.2,-9.9,-8.4,-6.0,-9.0,-10.3,-7.3,-7.6,-9.2,-8.8,-7.3,-8.9,-8.6,-8.4]
pred = [-7.1,-8.5,-9.2,-6.3,-9.0,-10.5,-7.4,-8.0,-8.8,-8.9,-7.5,-8.5,-8.9,-8.1]
r, p_value = stats.pearsonr(exp, pred)
print(r,p_value)

1

0.8960225522291245 1.4532438595517099e-05

spearman相关性系数R的计算

1
2
3
4
5

import scipy.stats as stats
exp = [-7.2,-9.9,-8.4,-6.0,-9.0,-10.3,-7.3,-7.6,-9.2,-8.8,-7.3,-8.9,-8.6,-8.4]
pred = [-7.1,-8.5,-9.2,-6.3,-9.0,-10.5,-7.4,-8.0,-8.8,-8.9,-7.5,-8.5,-8.9,-8.1]
rho, p_value = stats.spearmanr(exp, pred)
print(rho,p_value)

1

0.7836644591611479 0.000911883738030973

ΔG实验值与预测值的RMSD计算

1
2
3
4
5
6

from sklearn.metrics import mean_squared_error
exp = [-7.2,-9.9,-8.4,-6.0,-9.0,-10.3,-7.3,-7.6,-9.2,-8.8,-7.3,-8.9,-8.6,-8.4]
pred = [-7.1,-8.5,-9.2,-6.3,-9.0,-10.5,-7.4,-8.0,-8.8,-8.9,-7.5,-8.5,-8.9,-8.1]
mse = mean_squared_error(exp, pred)
rmse = mse**0.5
print(mse,rmse)

1

0.24714285714285705 0.4971346468944375

ΔG实验值与预测值的决策系数R²

1
2
3
4
5

from sklearn.metrics import r2_score
exp = [-7.2,-9.9,-8.4,-6.0,-9.0,-10.3,-7.3,-7.6,-9.2,-8.8,-7.3,-8.9,-8.6,-8.4]
pred = [-7.1,-8.5,-9.2,-6.3,-9.0,-10.5,-7.4,-8.0,-8.8,-8.9,-7.5,-8.5,-8.9,-8.1]
r2 = r2_score(exp, pred)
print(r2)

1

0.8026803535785574

FEP预测的MUE计算

FEP的MUE取自Flare/FEP，如下图4所示。

Figure 4. FEP计算完毕的网络图与活性的实验值与预测值比较

文献

1. D3R GC2. Drug Design Data Resource. https://drugdesigndata.org/about/grand-challenge-2
2. Schindler, C.; Rippmann, F.; Kuhn, D. Relative Binding Affinity Prediction of Farnesoid X Receptor in the D3R Grand Challenge 2 Using FEP+. J. Comput. Aided. Mol. Des. 2018, 32 (1), 265–272. https://doi.org/10.1007/s10822-017-0064-z.
3. 肖高铿,李青松. 云计算教程 | 用FLARE FEP计算相对结合自由能. http://blog.molcalx.com.cn/2020/07/17/flare-fep.html

4. Carvalho Martins, L.; Cino, E. A.; Ferreira, R. S. PyAutoFEP: An Automated Free Energy Perturbation Workflow for GROMACS Integrating Enhanced Sampling Methods. J. Chem. Theory Comput. 2021, acs.jctc.1c00194. https://doi.org/10.1021/acs.jctc.1c00194.

联系我们

我们不仅提供软件，还提供云计算资源与FEP计算服务。