摘要:本文以D3R CG2的Free Energy Set 2数据集的14个IC50跨度为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激动剂同系物:这些化合物具有公共的母核、结构相似,但IC50跨度高达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相关性系数R2=0.80 ± 0.08; RMSD = 1.31 kcal/mol。作者还同时用了Schrodinger/Prime MM-GBSA方法进行了ΔG预测,结果表明Pearson相关性系数R2=0.72 ± 0.11; RMSD = 9.31 kcal/mol。作者单独对其中的14个羧酸进行FEP+分析,结果表明RMSD仅降低了0.1。

FLARE算例 | 炼金术法FEP预测FXR激动剂的结合自由能-墨灵格的博客

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激动剂的结构、IC50(nM)及其预测的ΔG (kcal/mol)

Structure ID IC50(nM) ΔG(Exp) ΔG(Pred)
FXR_10 5640 -7.2 -7.1
FXR_12 58 -9.9 -8.5
FXR_74 655 -8.4 -9.2
FXR_76 41800 -6.0 -6.3
FXR_77 250 -9.0 -9.0
FXR_78 28.3 -10.3 -10.5
FXR_79 4150 -7.3 -7.4
FXR_81 2690 -7.6 -8.0
FXR_82 180 -9.2 -8.8
FXR_83 330 -8.8 -8.9
FXR_84 4540 -7.3 -7.5
FXR_85 297 -8.9 -8.5
FXR_88 540 -8.6 -8.9
FXR_89 735 -8.4 -8.1

FEP计算

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

图的生成

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

FLARE算例 | 炼金术法FEP预测FXR激动剂的结合自由能-墨灵格的博客

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

由于涉及到的14个分子IC50已知,本计算实验的主要目的时考察预测值与实验值的关系,因此生成图的时候采用了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,具有更好的精密度。

FLARE算例 | 炼金术法FEP预测FXR激动剂的结合自由能-墨灵格的博客

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

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

更重要的是,几个不同的生物等排体得到很好的排序,尤其是卤素取代苯环的衍生物的排序:1)苯环被噻吩取代对活性不利;2)苯环的氟取代对活性不利;苯环的氯、溴取代对活性有利。这是药物化学研究阶段非常重要的信息,而Flare/FEP可以准确在实验前的告诉你。FEP似乎对低活性化合物识别的特别好:预测为活性最差的5个化合物,确实也对应为IC50最高的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、分子量、打分函数相关性分析的统计学结果

Iterms SlogP MW Binding Affinity Score
Pearson Correlation (r) 0.717 0.644 -0.701
Spearman Rank Correlation(ρ) 0.713 0.626 -0.7
Root Mean Squared Deviation (RMSD) 1.07 589.438 54.99


与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实验值与预测值的决策系数R2

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所示。

FLARE算例 | 炼金术法FEP预测FXR激动剂的结合自由能-墨灵格的博客

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计算服务。