应用于蛋白-小分子动态对接的等变几何扩散模型-DynamicBind 评测

DynamicBind是一个用于蛋白和小分子复合体系构象采样及其打分的等变几何扩散网络。

DynamicBind利用等变几何扩散网络平滑蛋白-小分子复合结构的势能面，以便促进蛋白-小分复合物构象的转化，从而实现高效的采样。DynamicBind 可以从未结合的蛋白质结构中准确地生成配体结合的特异性构象，而无需带配体蛋白结构或大量MD采样，俗称柔性对接。

与柔性对接不同，DynamicBind对口袋构象的处理可以是侧脸的旋转，loop的摆动，二级结构的变化，口袋构象变化更大。

一、模型简介

DynamicBind是上海交通大学郑双佳研究员和星药科技Lu wei等人于2024年2月5日发表在nature communication上的工作，题目为：DynamicBind: predicting ligand-specific protein-ligand complex structure with a deep equivariant generative model。

关于蛋白和小分子的结合力/活性预测，一直以来，有很多的文章不断地发表，是药物设计深度学习的重点研究方向之一。但是这些方法使用的蛋白表示有序列或三级结构，分子的表示一般是smiles。但是在蛋白和小分子的结合力时，往往只有一个结果，而没有给出蛋白和小分子的结合位置，以及小分子的结合pose。

当然，也有其他的深度学习模型，像docking一样，基于输入的蛋白三级结构和小分子的smiles，生成小分子的结合pose，例如：diffdock, tankbind等；但是这些模型往往又是仅给出结合位置和分子的pose，没有给出打分。

另一方面，这些模型使用的蛋白结构式静态的，在对接过程中蛋白也是刚性的。不能处理诱导性口袋，或者小分子导致的口袋的侧链旋转。(注：tankbind就是星药科技开发的模型)

DynamicBind很好的将这两个问题结合起来，一起解决了，实现了蛋白质和小分子的动态对接，以及结合构象采样。

DynamicBind 可以从未结合的蛋白质结构中准确地生成配体结合的特异性构象（注意：是小分子结合时的蛋白特异性构象，即，不同的小分子有不同的构象），而无需带配体蛋白结构或大量MD采样，而且还可以进行隐蔽口袋的对接。

DynamicBind甚至可以实现小分子引导的DFG-in到DFG-out的转变。从这一个角度来看，是第一次AIDD在功能上超越了传统CADD的柔性对接的能力范畴，实现了更高程度的整合（端到端）。

如作者展示的下图，在动态对接过程中，甚至可以模拟DFG-out到DFG-in的变化。（gif 无法播放。。。）

二、模型结构

DynamicBind预测蛋白-小分子的复合结构，同时预测蛋白可能出现的构象变化。

输入的蛋白结构为pdb结构，可以是AlphaFold的预测结果。输入的小分子是smiles或者是sdf格式。

DynamicBind将rdkit生成的小分子构象随机的放置在蛋白周围，然后使用逐渐变小的时间步长，模型逐渐平移和旋转小分子，同时调整其内扭转角。 在仅改变配体构象的最初五个步骤之后，DynamicBind模型随后在剩余步骤中同时平移和旋转蛋白质残基，同时修改氨基酸侧链扭转角。 

在每一步的迭代中，小分子和蛋白的特征（侧链扭转角，ca原子坐标，小分子的键，原子特征等）都会输入到SE3等变模块，然后经过读出机制更新。如下图：

与基于扩散模型训练中采用的添加噪音方法不同，传统方法通过使用不同幅度的高斯噪声扰乱真实状态来生成诱饵，DynamicBind采用类似变形的变换来产生蛋白质诱饵。 

DynamicBind的类变形转化生成的诱饵通常满足基本的化学约束，并且可以使模型能够专注于学习生物物理相关的状态变化事件。在此过程中，天然构象逐渐转变为 AlphaFold 预测的构象。而传统的高斯噪声扰乱只能让模型学习到符合化学规则的状态（添加噪音前的状态）。

在采样效率方面，当蛋白质与两种不同的配体结合时，模型可以预测不同的holo构象，模型生成的是特异性结合构象。 DynamicBind模型可以在 20 步内预测结合状态下的蛋白质构象，而需要数百万步的全原子 MD 模拟才能找到相同的蛋白结合状态。

三、模型性能

3.1 小分子pose的准确性和蛋白口袋构象

训练、验证和测试集，使用的是pdbbind数据集，按照时间进行划分。使用的测试集是2019年以来的300个结构，其中包含了53个多肽。同时，作者使用策展的Major Drug Target 数据集中2020 以后的结果，进一步扩展测试集，里面包含了激酶，核酸受体，GPCR等家族的蛋白。

如下图，DynamicBind预测的小分子pose在不同的RMSD阈值下，比DiffDock, Glid，TankBind等其他AI模型和传统基于力场的glide等对接工具，均展现出更高的达成比例。

33%的pdbbind测试集和39%的Major Drug Target测试集中小分子的RMSD小于2A，如下图。

同样，上图c中，深色展示了在RMSD小于2A且clash score <0.35条件下，测试集的达成比例。浅色的是RMSD小于5A且clash score <0.3条件下，测试集的达成比例。

显然，DynamicBind具有更好的结果，意味着，DynamicBind生成的小分子pose不仅仅更接近于晶体结构（RMSD更小），而且，与蛋白口袋的冲突更少。也意味着，DynamicBind能够处理蛋白口袋的构象变化，生成更接近holo状态下的构象。

如下图，即使在初始口袋与晶体结构的RMSD很大的情况下，DynamicBind预测的口袋RMSD要显著小于AlphaFold 预测的结果。

DynamicBind预测的contact-LDDT (cLDDT) 的打分与小分子的RMSD之间相关关系较好。在最优构象方面，随着采样的数量增加，成功率增加，但是相同成功率条件下DynamicBind所需的采样构象相对更少。

这些结果都说明DynamicBind有能力预测口袋相应的构象变化，生成的口袋构象更接近晶体结构。

在SI信息中，即使在与训练集中小分子相似度小于30%，60%，90%的测试集中，DynamicBind 预测的小分子在不同RMSD阈值下的达成比例，也优于其他AI模型，同时优于包括glide等基于力场的传统docking方法。（即，DynamicBind有能力预测没有见过的小分子。）

但是在没有出现在训练集中的蛋白（序列相似度小于30%，60%，90%），不同小分子RMSD阈值下的达成率，DynamicBind 的结果不如传统的基于力场的glide等docking工具结果，如下图。

意味着，DynamicBind无法应用于没有出现过新的靶点，或者不如传统的基于力场的方法好。当然，其他AI方法的结果类似，例如Diffdock。

说明，DynamicBind不能预测没有见过的蛋白。

3.2 口袋的构象变化

传统的对接方法，需要先确定好口袋的构象，要确定口袋正确的构象。当时，这一步很难，因为不同的小分子，往往会诱导口袋形成不同的构象，因此往往需要重复按照小分子去研究口袋的对接grid。另一种方案是使用ensemble docking方法。

但是，DynamicBind可以自动的根据小分子调整口袋的构象，DynamicBind可以输入AlphaFold预测的结果，而不需要先确定口袋的正确的构象，也意味着，对于不同的分子，只需要输入一个口袋构象即可。这一点，确实是传统CADD工具没有实现的。

作为例子，作者利用激酶的DFG-in和DFG-out进行举例，6UBW和7U3S，都是c-MET。

这两个例子的AlphaFold预测的口袋与小分子结合后的口袋RMSD差距都很大，分别为9.44 Å和6.02 Å。

且这两个体系均未出现在模型的训练集中。但是DynamicBind预测的小分子RMSD分别仅为0.49 Å和0.51 Å，说明预测的小分子pose基本正确。而DynamicBind预测的口袋RMSD分别仅为1.97 Å和1.19 Å说明，口袋构象预测非常准了。

并且在结构可视化分析中，两个小分子分别诱导的DFG-in和DFG-out也能被DynamicBind预测出来。

关于其他测试集中DFG-in和DFG-out体系的79个PDB结构，DynamicBind预测的口袋RMSD要明显小于AlphaFold预测的口袋RMSD。且口袋的RMSD比AlphaFold预测要大部分都有提升0～2Å（这种提升就不明显），有的甚至达到了8Å（这就属于显著提升）。

在测试集中，作者展示更多不同口袋构象变化的例子，从侧链到loop，到二结构的变化，如下图，白色为AlphaFold预测结构作为DynamicBind的输入，绿色为DynamicBind预测的结合构象。

晶体结构中的小分子为黄色，DynamicBind预测的小分子pose为蓝色。

图中，负的Δpocket RMSD意味着：相较于AlphaFold预测结构，DynamicBind预测口袋构象更接近晶体结构。

负的Δclash意味着，在预测的复合结构中，DynamicBind预测结果发生碰撞的更少。

这些结果说明，DynamicBind可以预测复合体中口袋构象的各种层次变化，即包括了侧链，也包括loop，也包括二级结构，涉及ps到毫秒不同时间尺度的结构变化。

3.3 隐蔽性口袋动态对接

DynamicBind对隐蔽性口袋的动态对接，作者使用的体系是SET domain-containing protein 2 (SETD2) 体系，其隐蔽性口袋的诱发小分子为EZM0414。

DynamicBind在训练时，有使用到SETD2同源家族，但是没有EZM0414的隐蔽性口袋，因为在2019年前，EZM0414的晶体结构还没有解出。

如下图C，白色是AlphaFold预测结构，直接与EZM0414发生冲突。图e为DynamicBind预测的结果，图f将预测结果与晶体结构进行叠合。

从结果来看，DynamicBind确实实现了隐蔽性口袋的动态对接，不仅仅是成功的放置了小分子，而且找到了较为合适的口袋构象（EZM0414小分子的RMSD为1.4Å，对应的口袋的RMSD 为2.16 Å）。

3.4 虚拟筛选能力

DynamicBind也增加了预测小分子和口袋的结合能力的模块。作者使用含有2616化合物-蛋白对的测试集。结果显示，DynamicBind的ROC有0.68，要好于其他AI模型的结果。

以，murD为例，预测的复合物体系构象，能比AlphaFold预测结果的直接叠合更合理。如下图。

四、DynamicBind评测

4.1 环境安装

复制项目代码

git clone https://github.com/luwei0917/DynamicBind.git

创建conda环境

conda env create -f environment.yml

安装Relax环境，用于结构弛豫。

conda create --name Relax python=3.8
source activate Relax
conda install -c conda-forge openmm pdbfixer libstdcxx-ng openmmforcefields openff-toolkit ambertools=22 compilers biopython

下载作者预训练好的checkpoint DynamicBind: Predicting ligand-specific protein-ligand complex structure with a deep equivariant generative model

解压在./目录即可。

4.2 案例测试-ABL1

这里测试的ABL1的uniprot ID 为 P00520 。ABL1有多个结构，抑制剂小分子，其中既有底物位点的，也有变构位点的。我们选择的含有底物抑制剂和变构抑制剂的结构是5MO4，如下图所示。

但是作为输入，使用的是ABL1的另一个晶体结构，1OPK。我们尝试将5MO4的小分子，底物小分子和变构小分子，分别动态对接到1OPK中。然后比较，对接位置和pose，以及口袋的构象变化。接下来中，紫色均为1OPK，蓝色为5MO4。

在变构口袋构象方面，1OPK和5MO4差别不大，RMSD仅有0.234 Å，如下图：

在底物位点口袋构象方面，RMSD差别还是很明显的，RMSD为0.554，特别是DFG位置的两个氨基酸。特别注意的是，由于5MO4在DFG位置，以及QYG位置有氨基酸的确实未解出，所以，将1OPK的口袋更狭窄。将5MO4的小分子对接到1POK中，DynamicBind需要调整DFG位置和QYG位置氨基酸的构象。在结果中，值得特别关注。

需要对接的底物小分子（Nilotinib）smiles为：Cc1ccc(cc1Nc2nccc(n2)c3cccnc3)C(=O)Nc4cc(cc(c4)n5cc(nc5)C)C(F)(F)F，对应2D结构如下图：

需要对接的变构小分子（asciminib）smiles为：c1cc(ccc1NC(=O)c2cc(c(nc2)N3CC[C@H](C3)O)c4ccn[nH]4)OC(F)(F)Cl，对应2D结构如下图：

将变构小分子和底物小分子均写入ligand.csv中，如下图：

将1OPK，删除水，删除原有的配体小分子，然后保存为1OPK.pdb。将1OPK.pdb和ligand.csv放入到新建的./ABL1文件夹中。然后在./ABL1文件夹中新建一个test文件夹，用于保存对接结果。

然后开始使用DynamicBind进行动态对接。执行如下命令：

python run_single_protein_inference.py \
  ABL1/1OPK.pdb  ABL1/ligand.csv \
  --savings_per_complex 40 \
  --samples_per_complex 10 \
  --inference_steps 20 \
  --header test \
  --device 0 \
  --python ~/anaconda3/envs/DynamicBind/bin/python \
  --relax_python ~/anaconda3/envs/Relax/bin/python \
  --results ABL1 -p

程序会先下载一个checkpoint，esm2_t33_650M_UR50D.pt，放在./esm_models/checkpoints/路径中，2.4G。第一次运行需要的下载，时间会有些长。

（注：作者一直在更新这个github项目，每次运行的结果不太一致，上次我是2月10日运行的，当时貌似无需下载）

运行输出：

输出结果保存在./ABL1/test/文件夹中，文件夹目录为：

.
├── affinity_prediction.csv
├── complete_affinity_prediction.csv
├── index0_idx_0
└── index1_idx_1

由于我们是同时对接两个分子，在ligand.csv文件中包含两个smiles，所以./ABL1/test/文件夹中有两个子文件夹（index0_idx_0， index1_idx_1），分别为两个分子的对接结果。index0_idx_0为底物小分子，index1_idx_1为变构小分子。

affinity_prediction.csv 内容如下，包含了每一个分子的 affinity 对接打分。可以看到，底物口袋小分子的结合力预测打分为7.31，而变构小分子的结合力打分为7.71。

name,affinity
idx_0,7.31790423465729
idx_1,7.717542942210824

complete_affinity_prediction.csv是每个小分子的top10的pose打分及其affinity打分，内容如下：

name,rank,lddt,affinity
idx_0,1,0.6664215,7.3399496
idx_0,2,0.6738703,7.2115693
idx_0,3,0.66106814,7.306659
idx_0,4,0.6627517,7.3648543
idx_0,5,0.6653233,7.27874
idx_0,6,0.6571611,7.3133626
idx_0,7,0.6605489,7.357269
idx_0,8,0.6572863,7.3163543
idx_0,9,0.65351903,7.3229322
idx_0,10,0.6513432,7.3701806
idx_1,1,0.6083162,7.9174366
idx_1,2,0.6065392,7.737226
idx_1,3,0.6057902,7.6436768
idx_1,4,0.59344906,7.744586
idx_1,5,0.59857726,7.9449034
idx_1,6,0.5886727,7.6790857
idx_1,7,0.5797746,7.6272044
idx_1,8,0.5735723,7.406704
idx_1,9,0.5832283,7.7300243
idx_1,10,0.5790154,7.726334

在index0_idx_0和index1_idx_1文件夹中，包含了Top10的蛋白构象，以及对应的小分子构象，例如：top1的小分子构象文件和蛋白构象文件分别为：rank1_ligand_lddt0.61_affinity7.92_relaxed.sdfrank1_complex_lddt0.61_affinity7.92_relaxed.pdb。其中，0.61是构象打分，7.92是小分子-口袋结合力打分。

我这里想测试的是DynamicBind预测的小分子pose，以及口袋构象。检查DynamicBind是否能真的动态对接小分子，根据对接的小分子，给出相应更为合理的口袋构象。

4.1.1 底物口袋对接结果

首先，查看底物小分子的结果，这里我们展示DynamicBind认为的最有结果，即rank1。

DynamicBind对接的底物小分子位置和pose与5MO4晶体的非常类似，没有出现分子pose颠倒的情况，只是在三氟苯的一侧，偏差有些大（紫色为DynamicBind预测pose）。

在口袋的构象方面，hinge loop 的构象，DynamicBind预测的结果（紫色）与5MO4（蓝色）和1OPK（绿色）类似。如下图。

python movie_generation.py \
  ABL1/test/index0_idx_0/ 1+2 \
  --device 0 \
  --python ~/anaconda3/envs/DynamicBind/bin/python \
  --relax_python ~/anaconda3/envs/Relax/bin/python

python run_single_protein_inference.py \
  TEM1/1JWP.pdb  TEM1/ligand.csv \
  --savings_per_complex 40 \
  --samples_per_complex 30 \
  --inference_steps 20 \
  --header test \
  --device 0 \
  --python ~/anaconda3/envs/DynamicBind/bin/python \
  --relax_python ~/anaconda3/envs/Relax/bin/python \
  --results TEM1 -p