腾讯建设网站视频视频下载,wordpress目录分页,汕头有几个区几个县,金点子published at nature communication (2024.01.24) code link paper link
摘要
尽管在预测静态蛋白质结构方面取得了重大进展#xff0c;但蛋白质的内在动态性#xff0c;受到配体调节#xff0c;对于理解蛋白质功能和促进药物发现至关重要。 传统的对接方法#xff0c;常… published at nature communication (2024.01.24) code link paper link
摘要
尽管在预测静态蛋白质结构方面取得了重大进展但蛋白质的内在动态性受到配体调节对于理解蛋白质功能和促进药物发现至关重要。 传统的对接方法常用于研究蛋白质-配体相互作用通常将蛋白质视为刚性体。虽然分子动力学模拟可以提出适当的蛋白质构象但由于生物学相关平衡状态之间罕见的转换它们在计算上要求很高。现实中以5xco为例可以发现与配体结合的KRAS和单体KRAS有着很不一样的构象 本文提出了DynamicBind一种深度学习方法采用等变几何扩散网络来构建一个平滑的能量景观促进不同平衡状态之间的高效转换。DynamicBind准确地从未结合的蛋白质结构中恢复出特定于配体的构象而无需整体结构或广泛采样。 它在对接和虚拟筛选基准测试中展现出了最先进的性能实验揭示DynamicBind能够适应广泛的大型蛋白质构象变化并在未见过的蛋白质靶标中识别隐蔽口袋。
介绍
目前有一些预测结构方法已经很不错了比如机器学习和分子动力学或者蒙特卡洛的方法它们通过生成一个结构集合比如AF虽然它只使用蛋白质序列生成几个构象但它却已经有着非常准确的预测能力。 但问题是药物分子的治疗效果源自它们与目标蛋白质的某些构象的特异性结合从而通过改变这些蛋白质的构象景观来调节基本的生物活动也就是蛋白质本质上是动态的。 虽然大家也都知道动态性是很重要的但由于计算成本传统对接都是将蛋白质视为刚体或者最多是允许蛋白质部分灵活通常就是侧链扰动。这带来的结果就是当使用AlphaFold预测的无配体蛋白质(apoproteins)结构作为对接输入时将产生的配体姿态预测与配体结合的共晶全结构(holo-structures)对齐不佳。
本文提出DynamicBind一个为“动态对接”设计的几何深度生成模型。与将蛋白质视为大多数刚体的传统对接方法不同DynamicBind能有效地将蛋白质构象从其初始的AlphaFold预测调整到类全态(holo-like state)。该模型能够在预测期间处理广泛的大型构象变化例如在激酶蛋白中众所周知的 DFG-in 到 DFG-out 转换这是分子动力学(MD)模拟等其他方法难以克服的挑战。 模型通过学习一个漏斗形的能量景观来实现高效采样大型蛋白质构象变化的效率其中生物学相关状态之间的转换最小化了障碍。训练期间使用类似形态转换的创新手段生成假体(decoy generation)详见“DynamicBind架构”和“方法”部分。当前方法与Boltzmann生成器有相似之处因为它允许直接且高效地从学习的模型采样低能量状态。然而与通常仅限于其训练系统的传统Boltzmann生成器不同DynamicBind是一个可泛化的模型能够处理新的蛋白质和配体。
接下来文章分为六个部分首先概述DynamicBind模型然后将DynamicBind与当前对接方法进行baseline测试接着突出该方法在特定配体方式下采样大型蛋白质构象变化的能力然后具体演示它可以处理的构象变化范围通过案例研究展示其预测隐蔽口袋的能力最后展示其在使用抗生素数据集进行全蛋白质组虚拟筛选任务中的应用。
结果
结果总是一波大吹特吹总的来说就是配体姿态的预测rmsd较好碰撞较少pocket的rmsd较低设计的置信度有效效率高
DynamicBind架构
DynamicBind执行了“动态对接”这是一个在容纳大量蛋白质构象变化的同时进行蛋白质-配体复合物结构预测的过程。DynamicBind接受类似于天然未结合状态的结构在本研究中是AlphaFold预测的构象PDB格式以及小分子配体在几种广泛使用的格式如SMILES或SDF格式。在推理过程中模型使用RDKit生成的种子构象随机地将配体放置在蛋白质周围。然后在20次迭代的过程中更多细节见“模型架构”通过逐渐减小的时间步长模型逐渐平移和旋转配体同时调整其内部扭转角。在最初的五步中仅改变配体构象然后在剩余的步骤中模型同时平移和旋转蛋白质残基同时修改侧链Chi角。
如图1a所示每一步中蛋白质和配体的特征及坐标都被输入到一个SE(3)-等变的相互作用模块中。随后蛋白质和读出模块生成当前状态的预测平移、旋转和二面角更新。模型的更多细节在“蛋白质构象转换”中给出。 一个SE(3)-等变的模型能够保持其对旋转和平移的不变性或协变性这意味着如果输入数据在这种情况下是蛋白质和配体的三维结构经历了一个刚体变换例如整体旋转或平移模型的输出也会以相同的方式变换。 与在基于扩散的模型训练中采用的传统协议不同本文方法采用了一种类形态转换(morph-like transformation)来产生蛋白质变体decoys。 morph-like transformation: 平滑地从一个形状或结构转换到另一个形状或结构的过程用于生成在结构上连续且符合生物物理约束的蛋白质构象变体。使得模型在训练的时候可以更好地学习到结构的变化。 在这个过程中原始构象逐渐过渡到AlphaFold预测的构象。蛋白质的结构在许多方面受到高度限制残基通过肽键连接且键长受化学原理的制约。当使用高斯噪声生成变体时模型主要只学习恢复到最化学稳定的构象通常是添加噪声之前的构象。在当前任务中配体结合的整体构象是未知的最容易获得的蛋白质结构是AlphaFold预测的这通常与整体构象有显著不同。鉴于AlphaFold预测的结构通常已符合大部分化学约束预测仅用高斯噪声生成的变体上训练的模型如何准确预测生物相关的长时间尺度转换是一个挑战这是本文主要关注的问题。相比之下模型通过类形态转换生成的替代品通常满足基本的化学约束允许我们的模型专注于学习生物物理相关的状态变化事件。在无偏分子动力学模拟中如DFG“in”和“out”转换等亚稳态之间的转换由于全原子力场固有的现实而崎岖的能量景观而不频繁发生能量景观有过多的山峰和山谷想要跨过一个峰形成变化会比较困难。
相比之下本方法具有一个显著更加漏斗形的能量景观有效降低了生物学意义上的状态之间的自由能障碍。因此类似于其他Boltzmann生成器方法当前方法在抽样与配体结合相关的替代状态方面显示出显著提高的效率。为了阐明这些差异已经包含了一个示意图图1b。 Boltzmann: 利用深度学习技术特别是生成对抗网络GANs和变分自编码器VAEs来高效采样复杂系统能量景观中的低能态的工具。 fig1 全息状态(holo)以粉色表示初始apo和模型预测的构象以绿色表示。天然配体以青色表示预测的配体以橙色表示。模型接受蛋白质和配体的特征及其当前构象作为输入。输出结果包括预测的更新配体和每个蛋白质残基的全局平移和旋转配体的扭转角和蛋白质残基的chi角的旋转以及两个预测模块结合亲和力A和置信度得分D。在训练阶段模型旨在学习从apo样构象转换到全息构象的过程。在推理阶段模型迭代更新初始输入结构二十次。b示意图显示当蛋白质与两种不同的配体结合时我们的模型可以预测两种不同的全息构象。我们的模型可以在20步内预测出结合的蛋白质构象而寻找同一结合状态的所有原子MD模拟需要数百万步。 配体姿态预测更准确还优化了AF预测的构象
为了评估方法首先使用了PDBbind数据集并且按照以往的研究通过时间顺序的方式对训练集、验证集和测试集进行了分割。由于PDBbind测试集包含大约300个2019年的结构其中包括许多非小分子配体53个为多肽因此我们通过精选的主要药物靶标MDT测试集扩大了评估的范围。MDT集包括599个在2020年或之后存入的结构这些结构既有类药配体也有来自四大蛋白家族的蛋白质激酶、GPCRsG蛋白偶联受体、核受体和离子通道详细信息请参见“数据集构建”。这些蛋白家族代表了约70% FDA批准的小分子药物的靶标。
在测试中采用了一个更具挑战性和现实性的场景只使用AlphaFold预测的蛋白构象作为输入而不是使用全息结构作为输入文章还假设全息蛋白构象是不可获得的。之所以这么做是因为全息构象展示了与共结晶配体的强烈形状和电荷互补性这只会简化了配体姿态预测。相比之下apo构象或AlphaFold预测的构象可能与通过叠加晶体结构获得的移植配体发生冲突也就是将配体直接align到apo蛋白。
如图2a和2b所示DynamicBind预测的配体RMSD均方根偏差在各种阈值以下的案例数量超过了其他baseline。特别是它在PDBbind测试集上配体RMSD低于2Å5Å的比例分别为33%65%在MDT测试集上分别为39%68%。
评估模型时仅依靠配体的RMSD可能会偏向于深度学习模型如DiffDock、TankBind和DynamicBind因为它们对冲突的容忍度更高而可能对基于力场的方法如GNINA、GLIDE、VINA不利后者严格执行范德华力。显著的冲突可能会阻碍基于结构的药物设计中的相互作用分析掩盖关键的分子相互作用并使分子改进设计复杂化。因此本文使用配体RMSD和冲突得分按照Hekkelman等人定义来评估成功率。图2c显示了使用严格标准配体RMSD 2 Å冲突得分 0.35和更宽松的标准配体RMSD 5 Å冲突得分 0.5的成功率。在更严格的条件下DynamicBind的成功率0.33是最好的基线DiffDock0.19的1.7倍。此外DynamicBind已经展示了降低口袋RMSD相对于初始AlphaFold结构的能力即使在原始口袋RMSD较大的情况下见图2d。这一观察突出了当前方法能够处理大规模构象变化恢复全态结构而其他方法可能会遇到困难。鉴于我们的模型能够生成多样的构象我们开发了contact-LDDTcLDDT评分模块这一概念受到AlphaFold的LDDT分数的启发。该模块的目的是从预测输出中选择最合适的复合物结构。如图2e所示我们预测的cLDDT与实际配体RMSD有很好的相关性表明其在选择高质量复合物结构方面的有效性。以配体RMSD低于2Å为真正例的auROC得分为0.764。虽然我们的cLDDT评分函数有效但还有改进的潜力。完美的选择可以将我们的成功率从0.33提高到0.5如图2f所示。 a, bDynamicBind在预测PDBbind数据集a和主要药物靶标MDT数据集b中的配体姿态方面跨不同的RMSD阈值超越了其他方法。c深色和浅色阴影分别代表在严格配体RMSD 2 Å冲突得分 0.35和宽松配体RMSD 5 Å冲突得分 0.5标准下的成功率。d由DynamicBind预测的蛋白构象更接近天然态这可以通过绑定位点周围的口袋RMSD较低得到证明。eDynamicBind预测的contact-LDDTcLDDT得分与配体RMSD高度相关并且是预测真实配体RMSD低于2 Å的良好指标auROC 0.764。f随着生成样本数量的增加成功率提高。图c–f展示了结合PDBbind和MDT测试集的结果。关于单个数据集的结果以及在30%、60%和90%最大配体和蛋白序列相似性截止值下过滤的结果详细信息见补充图1-4。源数据提供为源数据文件。 即使没有理想的选择模型本方法在显著性能上也远远超过了DiffDock和顶级力场基方法GLIDE。由于Glide产生的样本数量变化较大通常是因为其过滤方案排除了不现实的构象Glide的最佳性能使用一条平线表示。这条线反映了由Glide最有效样本确定的成功率。DynamicBind的卓越性能源于其能够进行显著的蛋白质构象变化从而更好地匹配蛋白质和配体。
为了评估模型对新蛋白和配体的泛化能力分析了根据配体和蛋白序列的最大相似性划分的训练集的结果补充图3和4。这项分析揭示了DynamicBind在处理新配体时表现良好超过了其他方法但在处理新蛋白时效果较差此时它被具有预定义真实结合口袋的经典对接方法超越。其他深度学习方法在处理新蛋白时也显示出类似的下降这暗示了需要更大的训练集和改进的归纳偏差。此外考虑到在新蛋白上识别结合位点是一个活跃的研究领域深度学习方法在盲目全局对接中遇到的挑战可能是不同方法共有的。总的来说DynamicBind在新配体方面的熟练程度在药物发现中非常重要突出了其在识别对创建有效、特异性药物至关重要的蛋白质构象变化方面的潜力。
DynamicBind能够捕获配体特异性的蛋白质构象变化
传统的对接协议通常将蛋白质构象采样作为与对接过程分开的步骤进行。然而在许多情况下两个不同的配体可能适合于相互排斥的蛋白质构象。例如c-Met激酶可以采取两种不同的构象分别对应于活性状态和非活性状态通常被称为Asp–Phe–Gly (DFG)-in和DFG-out构象。DFG基序可以翻转出来随后阻塞或打开蛋白质的不同区域。在之前的对接模型中必须将蛋白质预设为正确的构象才有可能识别出配体的适当结合姿态。与此相反DynamicBind利用AlphaFold预测的蛋白质构象可以动态调整蛋白质构象找到能够容纳感兴趣配体的最优构象。作为一个代表性案例对于PDB 6UBW预测的配体RMSD为0.49Å口袋RMSD为1.97 Å而AlphaFold结构的口袋RMSD为9.44 Å。对于PDB 7V3S预测的配体RMSD为0.51 Å口袋RMSD为1.19 ÅAlphaFold 6.02 Å。这两个配体在训练集中之前都未曾见过。在定量分析中测试集79个PDB结构中只有七个蛋白质是同时采取了DFG-in和DFG-out两种构象这些通过激酶–配体相互作用指纹和结构KLIFS网络服务器进行标注。
图3f和g展示了这些蛋白质通过它们的UniProt ID标记从相同的初始结构开始随着I型抑制剂的结合逐渐向DFG-in构象移动并且在与II型抑制剂相互作用时倾向于DFG-out构象。此外图3h揭示了大多数预测的蛋白质结构与初始AlphaFold结构相比显示出较低的口袋RMSD。这些结果证明了DynamicBind能够捕获配体特异性的构象变化。这意味着即使特定构象与初始提供的蛋白质结构不同DynamicBind也能识别出能与蛋白质的其他可能构象很好结合的化合物。 AlphaFold预测的结构以白色显示晶体结构的蛋白质和配体分别以粉色和青色表示。模型的预测分别以绿色和橙色显示蛋白质和配体。天冬氨酸-苯丙氨酸-甘氨酸DFG残基的侧链以棍棒图显示。红色箭头突出显示了从AlphaFold结构到晶体结构的显著构象变化。输入构象是AlphaFold预测的构象。a 当配体84S (b) 与c-Met蛋白结合时蛋白质采用DFG-in构象。当配体5I9 (d) 与同一蛋白结合时蛋白质采用DFG-out构象。我们对两种配体的预测(c, e)与晶体结构非常吻合。配体RMSD分别为0.49 Å和0.51 Å。从初始AlphaFold到DFG-in和DFG-out的口袋RMSD改进分别为7.47 Å和4.83 Å。在测试集中七种蛋白通过其UniProt IDs识别包含DFG-in和DFG-out结晶全构象它们的初始AlphaFold和预测结构的口袋RMSD分别在(f. n 39)和(g. n 34)中显示用于DFG-in全构象和DFG-out全构象其中中心线标记中位数盒边缘指示上下四分位数须延伸至1.5倍四分位距个别点以点表示。h 所有79个PDB从AlphaFold到口袋RMSD改进的直方图。 DynamicBind覆盖多尺度蛋白质构象变化
DFG-in/out构象已被广泛研究通过使用集合对接ensemble docking技术其中同时利用两种构象的蛋白进行对接可以部分解决一些挑战。然而集合对接提高了计算成本并且可能不适用于较少被表征的构象。 在本节中提供了从皮秒级到毫秒级跨六种不同构象变化的综合分析每种变化都通过我们的PDBbind测试集中的一个案例来示例。在图4中晶体结构以粉色表示AlphaFold结构以白色表示我们的预测以绿色表示。天然配体以青色表示我们预测的配体以橙色表示。Δpocket RMSD衡量的是预测蛋白结构与AlphaFold结构之间的口袋RMSD差异基于与晶体结构的比较。负的Δpocket RMSD表明与晶体结构相比预测的结构与AlphaFold预测的结构更为接近。Δclash衡量的是预测的蛋白-配体对与AlphaFold结构中移植配体之间的冲突得分差异。负的Δclash表示预测复合物中的冲突更少。在图4a中天然配体与叠加的AlphaFold结构的一个侧链发生冲突在预测中这个侧链向天然构象移动从而解决了冲突。在图4b中AlphaFold结构中的一个酪氨酸挡住了一部分口袋在预测和天然结构中这部分口袋变得可访问。在图4c中一个灵活的环与配体相交它在预测中移开与天然结构一致。在图4d中靠近配体结合位点的α螺旋转变为环。在图4e中在热休克蛋白Hsp90α中观察到一个重要的二级结构运动从关闭状态转变为开放状态。在图4f中AKT1激酶的两个域聚合形成了一个之前不存在的口袋。综上所述当前模型能够预测与配体结合相关的多种类型的构象变化当配体结合口袋在AlphaFold预测的构象中不够宽敞或未形成时。
DynamicBind揭示了对药物发现至关重要的隐蔽口袋
蛋白质的动态特性经常会产生隐蔽口袋。这些在蛋白质动态过程中出现的隐蔽口袋可以揭示在静态结构中未发现的可药用位点从而将先前被认为“不可药用”的蛋白质转变为潜在的药物靶标。我们使用含SET域的蛋白2SETD2一种组蛋白甲基转移酶作为案例研究展示了DynamicBind在揭示这些隐蔽口袋方面的实用性。SETD2对于多发性骨髓瘤MM和弥漫大B细胞淋巴瘤DLBCL的治疗至关重要它有一个被一种高度选择性的化合物EZM0414靶向的隐蔽口袋该化合物目前正在进行I期临床试验。如图5a、b所示训练集中所有的SETD2同源体根据蛋白质Smith–Waterman相似性超过0.4定义都与S-腺苷甲硫氨酸SAM或Sinefungin类似物共结晶用线条表示。Sinefungin及其类似物通过占据SAM位点广泛抑制甲基转移酶使得对SETD2的选择性抑制变得具有挑战性。在2019年之前没有任何SETD2或其同源体与在EZM0414靶向的位点上结合的化合物一起结晶的结构。因此模型没有在任何结合到这个新发现位点的化合物的结构上进行训练。在图5c中AlphaFold结构及其表面以白色显示。隐蔽位点看起来被阻塞导致与移植的EZM0414发生大量冲突。图5d确认EZM0414为一个未见过的配体即使是最相似的Tanimoto配体也与EZM0414有显著偏差。图5e显示了模型预测的蛋白质-配体复合物结构以SETD2的AlphaFold预测结构和EZM0414的SMILES表示作为输入。图5f将预测与SETD2-EZM0414复合物的晶体结构PDB 7TY2叠加。结果的配体RMSD为1.4 Å口袋RMSD为2.16 Å。此外补充信息中包含了几个与训练集序列相似性低的Cryptosite数据集的案例。
DynamicBind在抗生素基准测试中实现了更好的筛选性能
在基于靶标的药物发现中筛选潜在药物候选物和反向筛选为特定化合物识别蛋白质靶标都至关重要。这些过程需要准确预测蛋白质和化合物之间的结合亲和力即在蛋白组水平上两者相互作用强度的度量。因此本模型中添加了一个亲和力预测模块使用从PDBbind数据集中获得的实验测量的结合亲和力数据进行训练。为了在现实世界的虚拟筛选场景中评估DynamicBind使用了一个最近发布的抗生素实验基准。这个数据集包括2616个蛋白质-化合物对这些在训练阶段中都未遇到过。它特别包含了与218种活性抗菌化合物配对的大肠杆菌必需蛋白质组中的12种蛋白质。图6a显示DynamicBind不仅超越了常见的对接方法如VINA和DOCK6.9还超过了最佳的基于机器学习的重新评分方法实现了接收器操作特征曲线下面积auROC的平均值为0.68。基线数据直接来源于基准论文。这种性能提升归因于DynamicBind的动态对接能力它优化了AlphaFold结构使其更接近自然状态从而实现了更精确的结合亲和力估计。如图6b所示预测的蛋白质murD结构更紧密地围绕配体形成了在初始AlphaFold结构中不可能形成的更多相互作用。这一在抗生素基准测试上的评估与我们在PDBbind测试集上对结合亲和力预测的基准测试结果一致补充表1其中DynamicBind一贯优于传统的对接方法和基于深度学习的刚性对接方法。这些结果表明DynamicBind凭借其结合亲和力预测能力在蛋白组水平的虚拟筛选应用中展示出显著的潜力。 总共12页的论文6页实验验证好太好怎么这么好真的很多话。但工作看起确实不错。
结论
DynamicBind将传统上分开的两个步骤——蛋白质构象生成和配体姿态预测——统一到一个框架中。作为一种端到端的深度学习方法它在抽样广泛的蛋白质构象变化方面比传统的分子动力学MD模拟快几个数量级。与需要预定义结合口袋的传统对接方法不同DynamicBind能够执行全局对接这一功能在结合口袋尚未被确定时变得至关重要。这些优势使DynamicBind能够虚拟筛选与隐蔽口袋结合的化合物。这类化合物很可能仅与目标蛋白质专一结合因此可能最小化副作用。此外DynamicBind可以预测新的药物候选物是否可能与非预期的蛋白质靶标结合或在通过表型筛选发现活性化合物时帮助识别结合靶标。
尽管DynamicBind在baseline中展示了最先进的性能但仍有改进的空间特别是在提高其对于与训练集中序列同源性低的蛋白质的泛化能力方面。作为一个数据驱动的模型它从冷冻电镜Cryo-EM方法的快速进步中受益匪浅。这些技术进步将扩大我们训练数据的多样性和全面性以更快的速度提供更多样化的蛋白质-配体复合物构象。还有潜力通过利用大量的非结构性结合亲和力数据来改进DynamicBind这些数据目前比结晶结构更为丰富。通过采用类似于AlphaFold的自我蒸馏方法我们可以通过整合之前只有亲和力数据可用的蛋白质-配体对的高信度预测来扩展我们的训练集。
总之DynamicBind提出了一种“动态对接”方法来研究蛋白质-配体相互作用将其与将蛋白质视为静态和计算要求高的分子动力学MD模拟的传统对接方法区分开来。它对大规模蛋白质动力学的处理能力对于发现药物分子特别是那些针对隐蔽口袋的药物分子具有特别重要的意义。此外DynamicBind生成的特定于配体的蛋白质构象可能提供有价值的洞察揭示配体对蛋白质的影响有助于阐明结构-功能关系增强对机制理解。
方法
总览
我们的模型是一个基于E(3)-等变的、扩散驱动的图神经网络采用了粗粒度表示。
E(3)-等变模型根据对3D空间中输入x进行旋转和宇称即空间的镜像反射操作来变换输出y。研究表明等变模型在结构的训练数据量可以少1000倍的情况下仍能在大块水的结构上获得更优的结果(Research has demonstrated that equivariant models can be trained with 1000 times less data while yielding superior results on the structures of bulk water)。尽管冷冻电子显微镜和晶体学方面取得了巨大进步现有的蛋白质-配体复合物数据库仍相对有限只扩展到数万的规模。因此需要一个高效的模型能够辨别最相关的信息并避免在移动或旋转整个结构时不成立的表面信息。实现SO(3)对称性描述三维空间中物体旋转不变性的一种方式的传统方法涉及仅使用或预测不变量如contact map。然而contact/distance map 并不总是与物理上可行的配置相关联。例如某个残基可能被预测与两个相距甚远的原子接触。此外接触图可能忽略了手性这在药物发现中是一个重要方面。
作为一个基于扩散的模型DynamicBind通过逐渐扭曲原生构象的过程进行训练使模型学会如何恢复正确的构象。扭曲原始配置通常涉及向原子添加平移位置的高斯噪声。由于化学键施加的键距约束和范德华力强加的排斥体积效应当扭曲相对较小时从这种扭曲中恢复是直接的。然而我们观察到仅添加高斯噪声不足以训练出能够预测从一个生物学上有意义的配置到另一个配置转换的模型。为了解决这个问题我们引入了一种形态变化它在晶体蛋白结构和AlphaFold预测的结构之间插值从而降低了诸如AlphaFold预测的构象和配体结合的全构象之间的过渡障碍。与其他训练评分函数的生成模型不同 S θ ( x , t ) ≈ ∇ l o g p t ( x ) S_θ(x,t) \approx \nabla logp_t(x) Sθ(x,t)≈∇logpt(x)直观上可以将它想象成一个指向原始数据分布高密度区域的指针告诉模型如何更新当前噪声数据以更接近原始分布本文扩散架构旨在将扰动的结构直接映射回原始构象类似于一致性模型。模型的输出表示为 f θ ( x t , t ) − φ ( x t , t ) f_θ(x_t,t)-φ(x_t,t) fθ(xt,t)−φ(xt,t)其中 φ ( x t , t ) φ(x_t,t) φ(xt,t)代表添加到原生构象的形态变化。
传统方法使用全原子表示明确建模每个原子的坐标。然而由于它们通过化学键的连接原子并不独立移动且局部几何是高度约束的——例如苯环通常是平的。为了减少这些非物理配置的自由度本文采用了蛋白质和配体的粗粒度表示。在模型中每个蛋白质残基由一个节点表示带有两个向量——坐标和方向以及侧链二面角。更多细节在“特征化”中提供。对于配体每个重原子由一个节点表示这些节点以外显到内在的方式变换其中扭转角的变化转换为笛卡尔坐标的变化。更多细节可以在“配体构象的变换”中找到。值得注意的是尽管是粗粒度表示所有非氢原子的坐标仍然可以以一对一的方式映射。
我们模型的输入是蛋白质和配体的当前构象。输出包括对 k l k^l kl个标量扭转角和每个配体的两个平移-旋转向量的预测更新以及对 k i p k^p_i kip个侧链标量二面角和每个蛋白质残基的主链的两个平移-旋转向量的更新。更多细节可以在“蛋白质构象的变换”中找到。此外模型产生两个标量输出一个是用来估计原生构象的程度通过cLDDTcontact-LDDT评估另一个是预测蛋白质和配体之间的结合亲和力。
特征化
在我们的模型中配体被表示为带属性的图 G l ( V l , E l ) G^l (V^l, E^l) Gl(Vl,El)其中每个节点 v i l ∈ V l v^l_i \in V^l vil∈Vl代表一个重原子边代表芳香键、单键、双键或三键。配体图的节点特征包括原子序数、手性、度数和形式电荷。除了键类型边的长度嵌入也被用作标量边特征。
蛋白质图被表示为 G p ( V p , E p ) G^p (V^p, E^p) Gp(Vp,Ep)其中每个节点 v i p ∈ V p v^p_i \in V^p vip∈Vp对应于 C α C_\alpha Cα 位置的一个残基。蛋白质图中的节点特征包括氨基酸类型、来自esm7的语言模型嵌入以及侧链二面角它们被表示为7 x 2维的零填充标量特征五个可旋转的chi角 [chi1, chi2, …, chi5] 和两个对称的chi角 [altchi1, altchi2] 对于每个氨基酸并且这些角度被转换为正弦和余弦值。为了确保给定结构的侧链角度的唯一性统一处理为 [max(chi1,altchi1), min(chi1,altchi1), max(chi2,altchi2), min(chi2, altchi2), chi3, chi4, chi5]。此外主链的方向被表示为两个单位向量特征即 x N − x C α ∥ x N − x C α ∥ \frac{x_N - x_{C_\alpha}}{\|x_N - x_{C_\alpha}\|} ∥xN−xCα∥xN−xCα和 x C − x C α ∥ x C − x C α ∥ \frac{x_C - x_{C_\alpha}}{\|x_C - x_{C_\alpha}\|} ∥xC−xCα∥xC−xCα。对于边键长用作标量特征。我们对氨基酸的特征化使模型能够推断出所有重原子的位置。
模型结构
DynamicBind 是一个图神经网络它使用等变和不变特征。它根据 e3nn 库中的定义使用不可约表示irreps的张量积来传播信息。 “不可约表示”Irreducible Representation通常缩写为 irrep是指一个群的表示这个表示不能被分解为更小的、非平凡的表示。群是一个数学结构它包含了一系列元素和一个满足某些条件的操作例如旋转、反射等。群的表示是一种将群的抽象元素映射到矩阵的方式使得群操作可以通过矩阵运算来实现。 节点和边的输入标量特征与扩散时间的正弦嵌入相连接然后由不同的多层感知机MLPs编码。对于蛋白质节点氨基酸的两个单位向量特征与新的标量表示结合形成交互层的初始特征。与 DiffDock 类似在图传播过程的每一步中配体和蛋白质图经历一次内部交互和一次外部交互。在配体的内部交互中每个配体原子的表示通过其他在 5 Å 距离内的配体原子进行更新。对于蛋白质每个氨基酸由在 15Å 距离内的其他氨基酸更新。为了减少模型的训练运行时间和内存使用每个残基允许的最大邻居数为 24。外部交互的边缘是基于是否有氨基酸在 ( 3 σ t r 12 ) (3σ_{tr} 12) (3σtr12) Å 的距离内从任何配体原子其中 σ t r σ_{tr} σtr 是当前扩散平移噪声的标准差。这个动态截止设计用来确保即使在配体远离受体时当 σ t r σ_{tr} σtr 较大时也能存在相互连接。 当链接图确定之后每个节点的消息是通过TensorProductLayer更新的。具体来说对每一个属于 C α C_\alpha Cα类别的节点 h a ← h a ⨁ c ∈ { l , r } B N ( c a , c ) ( 1 ∣ N a ( c ) ∣ ∑ b ∈ N a ( c ) Y ( r a b ) ⊗ ψ a b h b ) with ψ a b ψ ( c a , c ) ( e a b , h a 0 , h b 0 ) h_a \leftarrow h_a \bigoplus_{c \in \{l, r\}} BN^{(c_a,c)} \left( \frac{1}{\sqrt{|\mathcal N_a^{(c)}|}} \sum_{b \in N_a^{(c)}} Y(r_{ab}) \otimes_{\psi_{ab} }h_b\right) \quad \text{with} \quad \psi_{ab} \psi^{(c_a,c)}(e_{ab}, h_a^0, h_b^0) ha←hac∈{l,r}⨁BN(ca,c) ∣Na(c)∣ 1b∈Na(c)∑Y(rab)⊗ψabhb withψabψ(ca,c)(eab,ha0,hb0) 此处 h a h_a ha代表每个节点的特征 h a 0 h_a^0 ha0代表它的标量特征。 N a ( c ) \mathcal N_a^{(c)} Na(c) 代表 一个 c c c类别的节点 a a a的邻居包括配体或者蛋白。球谐函数用 Y Y Y 表示 B N BN BN 代表等变的批量归一化。模块 Ψ Ψ Ψ是一个包含可学习权重的多层感知器MLP这些权重基于边嵌入 e a b e_{ab} eab和标量特征 h a 0 h_a^0 ha0, h b 0 h_b^0 hb0来计算。 球谐函数Spherical Harmonics是一组在球面坐标系中定义的正交函数通常用于解决球对称的问题 在最终交互层之后节点表示被用来生成输出。为了生成 cLDDT、结合亲和力、配体的平移和旋转预测采用了每个配体原子与配体的几何中心的卷积 v 1 ∣ V l ∣ ∑ a ∈ V l Y ( r o a ) ⊗ ψ o a h a with ψ o a ψ ( e o a , h a 0 ) v \frac{1}{|\mathcal V^l|} \sum_{a \in \mathcal V^l} Y(r_{oa}) \otimes \psi_{oa} h_a \quad\text{with}\quad \psi_{oa} \psi(e_{oa}, h_a^0) v∣Vl∣1a∈Vl∑Y(roa)⊗ψoahawithψoaψ(eoa,ha0) 其中 e o a e_{oa} eoa 是配体的几何中心与配体节点 a 之间的边嵌入。输出 v v v 包含 144 个偶数标量2 个奇宇称向量和 2 个偶向量。这些标量用于预测 cLDDTD和结合亲和力的负对数A后者以浓度单位表示。 奇向量用于预测配体的移动或者平移偶向量用于预测配体的旋转 在这里 s t s_t st是扩散时间的正弦嵌入为了数值稳定性加上了 eps 1 0 − 12 10^{-12} 10−12。根据Jing等人的方法模型预测了配体每个可旋转键的标量扭转更新。对于键b扭转更新 T b l T^l_b Tbl 是通过将半径图上的每个原子与键的中心 o o o进行卷积生成的。 为了预测蛋白的构象变换我们需要对每个蛋白节点更新侧链 chis平移和旋转。这些操作都是由最后每个氨基酸的交互表征 h i h_i hi产生的 此处 T i p T^p_i Tip是一个五维度的标量输出代表了扭转角[chi1, chi2, …, chi5]的更新。
配体构象的变换
为了更新配体构象我们使用了一个统一的全局平移 t r l ∈ R 3 \mathbf{tr}^{l} \in \mathbb{R}^3 trl∈R3和旋转 R l ∈ R 3 × 3 \mathbf{R}^{l} \in \mathbb{R}^{3 \times 3} Rl∈R3×3。配体的所有原子将围绕配体的几何中心同时被平移和旋转几何中心计算公式为 x l ˉ 1 n ∑ i x i ′ \bar{\mathbf{x}^{l}} \frac{1}{n} \sum_{i} \mathbf{x}_i^{\prime} xlˉn1∑ixi′其中 n n n是配体重原子的总数 x i l \mathbf{x}_i^{l} xil表示第 i i i个原子的位置向量。具体来说变换后的位置向量 x l \mathbf{x}^{l} xl是通过 x l R l ( x l − x l ˉ ) x l ˉ t r l \mathbf{x}^{l} \mathbf{R}^{l}(\mathbf{x} ^{l}- \bar{\mathbf{x^{l}}}) \bar{\mathbf{x}^{l}} \mathbf{tr}^{l} xlRl(xl−xlˉ)xlˉtrl获得的。
在平移和旋转之外扭转角度对配体构象的修改也是至关重要的。然而修改扭转角可能会扰乱质心的位置。为了解决这个问题Corso等人展示了在更新扭转角后执行RMSD对齐可以确保扭转更新与平移更新是正交的因此可以将扭转更新的后果从平移更新和旋转平移更新中解耦。总的来说更新后的配体姿态是通过 x l RMSDAlign ( ( T 0 l ∘ … T k l ) ( x l ) , R l ( x l − x l ˉ ) x l ˉ t r l ) \mathbf{x}^{l} \text{RMSDAlign}((\mathbf{T}^{l}_{0} \circ\ldots \mathbf{T}^{l}_{k})(\mathbf{x}^l),\mathbf{R}^{l}(\mathbf{x}^{l} - \bar{\mathbf{x}^{l}}) \bar{\mathbf{x}^{l}} \mathbf{tr}^{l}) xlRMSDAlign((T0l∘…Tkl)(xl),Rl(xl−xlˉ)xlˉtrl)获得的其中 T k l \mathbf{T}^{l}_{k} Tkl是扭转旋转。
蛋白质构象的转换
根据AlphaFold我们使用C α _{\alpha} α作为残基节点来执行全局平移和旋转。此外模型还预测侧链扭转角的更新。对于180度旋转对称的侧链部分考虑对称性在推理阶段是不必要的但我们在训练期间引入对称侧链扭转特征以正确计算损失函数。由于C α _{\alpha} α的位置与侧链扭转角是独立的它不影响残基级别的平移和旋转。因此我们可以以任何顺序执行旋转平移和扭转旋转。最后每个蛋白质残基的更新构象表示为 x i p ( T i , 0 p ∘ … T i , k p ) R p ( x i p − x i , C α p ) x i , C α p t r i p \mathbf{x}^{p}_{i} (\mathbf{T}^{p}_{i,0} \circ \ldots \mathbf{T}^{p}_{i,k})\mathbf{R}^{p}(\mathbf{x}^p_{i} - \mathbf{x}^p_{i,C_{\alpha}}) \mathbf{x}^p_{i,C_{\alpha}} \mathbf{tr}^{p}_{i} xip(Ti,0p∘…Ti,kp)Rp(xip−xi,Cαp)xi,Cαptrip其中 T i , k p \mathbf{T}^{p}_{i,k} Ti,kp是第 i i i个残基的侧链扭转旋转。
训练和推理
在训练过程中输入的是通过对原始蛋白结构类形态转换后的替代结构配体结构是做了高斯加噪的。预期的输出是去噪的操作。输入蛋白质结构在时间 t t t被定义为 x i p ϕ ( x h o l o , t ) \mathbf{x}^{p}_{i} \phi(\mathbf{x}^{holo}, t) xipϕ(xholo,t)。具体来说对于第 i i i个氨基酸Kabsch算法用于计算C α _{\alpha} α的平移 t r i ∗ \mathbf{tr}^{*}_{i} tri∗和旋转 r o t i ∗ \mathbf{rot}^{*}_{i} roti∗从而可以对齐主链原子N−C α _{\alpha} α−C的主链使得整体结构与apo结构一致 考虑到扭转角之间的差异我们可以描绘出第 i i i个残基的构象变化 这里 T i , k ∗ T i , k a p o − T i , k h o l o T^*_{i,k} T^{apo}_{i,k}-T^{holo}_{i,k} Ti,k∗Ti,kapo−Ti,kholo 以弧度为单位 R i ∗ R^*_i Ri∗是 r o t i ∗ rot^*_i roti∗的旋转矩阵。在任何给定的时刻我们的目标是用一个因子扰动蛋白质结构这个因子表示为 u ( t ) u(t) u(t)使得扰动后的数据是在holo结构和apo结构之间的中间状态 其中 τ m i n p \tau^p_{min} τminp和 τ m a x p \tau^p_{max} τmaxp表示扩散噪音的参数。
为了克服训练和推理之间的分布漂移主要是由于推理部分使用的是RDKit生成的构象我们用与真实姿态 x g t x^{gt} xgt相匹配的构象 x 0 l x_0^l x0l替换了训练目标。在时间 t t t输入配体的姿态是一个随机扰动的构象 此处 x ˉ 0 l \bar x^l_0 xˉ0l是 x 0 l x^l_0 x0l的集合中心 p ( ω ) p(\omega) p(ω)是SO(3)上的各向同性高斯分布而 ω ^ \hat \omega ω^是通过随机采样生成的单位向量。
该网络通过八个损失函数来训练。总损失函数可以定义如下 其中 L t r l , L r o t l , L T l \mathcal L^l_{tr}, \mathcal L^l_{rot}, \mathcal L^l_{T} Ltrl,Lrotl,LTl是配体平移旋转和扭角的损失。 L t r p , L r o t p , L T p \mathcal L^p_{tr}, \mathcal L^p_{rot}, \mathcal L^p_{T} Ltrp,Lrotp,LTp则是蛋白残基的。 L A \mathcal L_{A} LA是亲和力损失 L D \mathcal L_{D} LD是contact-LDDT损失。用于计算真实cLDDT的距离差是 d ∣ d ( x 0 l , x h o l o ) − d ( x t l , x t p ) ∣ d|d(x_0^l,x^{holo})-d(x^l_t,x^p_t)| d∣d(x0l,xholo)−d(xtl,xtp)∣(更多关于cLDDT得分计算的细节可以在“评估指标”中找到。
如果一个旋转向量 u u u和另一个旋转向量 v v v方向相反且 ∥ u ∥ ∥ v ∥ 2 π ∥u∥∥v∥2π ∥u∥∥v∥2π那么 u u u代表与 v v v有着相同的旋转。所以我们在计算旋转损失时取正向和反向方向损失的最小值。扭曲角损失是通过计算预测值与加入的扭曲角噪声之间角度差的余弦来计算的。完整的训练程序可以在补充算法1中看到。
在推理过程中我们使用RDKit生成的配体结构和AlphaFold预测的蛋白质结构作为初始复合体构象。复合体结构通过20步更新。为了防止最终构象陷入局部最小值每一步都加入了一点随机噪声来去噪配体姿态。对于每对我们进行40次采样并根据cLDDT对可能的构象进行排序和评级。
我们还注意到预测的结合亲和力的加权平均值比实验测定的结合亲和力是更准确的估算器见补充表1。预测的cLDDT值用于完整的推断过程相关的推断过程可以在补充算法2中找到。
DynamicBind拥有63.67百万参数并且在八块Nvidia A100 80GB GPU上训练了5天。
评估指标
为了评估蛋白质和配体在预测复合物结构内的相互作用我们确定了分子间本质接触形成的程度。我们采用了类似于先前的局部距离差异测试LDDT评分的定义用于量化预测蛋白质结构的本质性。contact-LDDTcLDDT评分是通过考虑配体原子和蛋白质原子之间小于15埃的距离来计算的。距离差异是在基准真实结构和预测复合物结构之间确定的同时考虑了对称性。最终的cLDDT评分是从平均保守距离中得出的跨越四个容忍阈值0.5、1、2、和4埃。
为了评估预测蛋白质结构与结合口袋周围的本质蛋白质结构的偏差我们计算了口袋均方根偏差口袋RMSD。这是使用参考配体原子周围5埃内的蛋白质原子来执行的。最初预测的蛋白质结构与晶体蛋白质结构对齐。随后计算了预测口袋原子与晶体口袋原子之间的RMSD。
类似于AlphaFoldclash分数是通过计算配体原子和蛋白质原子之间小于4埃的所有距离的范德华重叠的均方根RMS来计算的。它按照以下方式计算 其中 N N N是需要考虑的距离的数量。
数据集构建
训练和测试数据集是建立在PDBbind2020数据库之上的该数据库包括19,443个晶体结构的蛋白质-配体复合物的策划集合每个都配有一个实验测量的结合亲和力。我们使用2019年前存放的结构进行训练和验证而2019年存放的用于测试。每个蛋白质都与对应相同蛋白质序列的AlphaFold预测结构对齐。对齐的AlphaFold结构和晶体结构用于生成训练样本和主要靶标的测试集通过morph-like插值来生成。Major Drug Targets (MDT)测试集是根据以下标准构建的PDBs存放在2020年或之后蛋白质属于四大主要药物靶标类别之一 - 激酶、GPCR、核受体和离子通道与晶体结构相比alphaFold-predicted protein structures具有超过2埃或口袋LDDT低于0.8的结构配体是分子量在200到650道尔顿之间的药物样小分子最多来自单一研究的10个PDBs被包含。这些标准确保了测试集是具有挑战性的并且初始蛋白质数据从本质上是完整的包括广泛的蛋白质靶标。此外它防止了少数蛋白质占据整个测试集因为某些研究存放了大量的PDBs这些PDBs是同一蛋白质共结晶的不同配体的结构。
Baseline
我们在PDBbind测试集303个配体-受体对和主要药物靶标MDT测试集599个配体-受体对上使用下面列出的不同对接方法进行了对接实验。对接配体是从共晶结构中提取出来的没有改变它们的原子坐标而对接受体结构则是由AlphaFold预测的。我们使用了一个有对称性意识的方法spyrmsd包中的symmrmsd函数来进行所有RMSD的计算。
Autodock VINA刚性对接
在Autodock Vina中配体从SDF格式转换为PDBQT格式是通过Meeko 2.0.0完成的。蛋白质准备是使用ADFR套件1.0中的“prepare_receptor”命令进行的。对接盒是围绕天然配体自动定义的所有六个面的默认缓冲区为4Å。盒子的中心是天然配体的质心。由于硼原子不是有效的AutoDock原子类型含有该原子的配体无法对接。因此PDBbind数据集中有301个配体-受体对MDT数据集中有597个配体-受体对在VINA刚性对接中有输出。
Autodock VINA柔性对接
与VINA刚性对接相比VINA柔性对接中有一个额外的柔性受体准备步骤。它是通过一个叫做“prepare_flexreceptor.py”的python脚本执行的该脚本可以在https://github.com/ccsb-scripps/AutoDock-Vina/tree/develop/example/autodock_scripts上找到。通过这一步蛋白质PDBQT格式文件被分成两个PDBQT格式文件一个用于刚性部分另一个用于柔性侧链。对于Vina Flex模式必须预先确定柔性侧链。我们识别了所有侧链原子与配体原子之间距离小于5Å的残基作为柔性残基。在这种模式下蛋白质主链保持刚性。配体准备和网格盒设置与VINA刚性对接一致。
GNINA刚性对接
GNINA的配体输入文件是PDBQT格式的是在添加了RDKit的氢原子后使用OpenBabel创建的。蛋白质输入文件是PDB格式文件。网格盒设置与VINA刚性对接一致。对于PDBbind数据集所有的配体-受体对都有对接输出。对于MDT数据集有1对没有输出因为原始PDB文件中的配体PDB ID: 8HMU没有完全解析缺少原子。
GLIDE
GLIDE是Schrödinger软件中的刚性蛋白质对接模块。配体的准备是通过使用LigPrep模块完成的。蛋白质准备是通过使用Protein Preparation Wizard模块完成的。网格文件是通过Receptor Grid Generation模块生成的内部盒子为10Å自动外部盒子围绕配体所有六个面的默认缓冲区为4Å中心位于配体的质心。然后进行了SP精度对接。PDBbind数据集中的一些配体是多肽无法通过LigPrep模块处理。此外与口袋原子严重冲突的配体在对接过程中没有输出姿态。因此PDBbind数据集中有266个配体-受体对MDT数据集中有472个配体-受体对在GLIDE刚性对接中有输出。
诱导拟合对接
Schrödinger软件中的诱导拟合对接IFD模块为用户提供了一个蛋白质柔性对接功能。与VINA和GNINA不同的是不仅残基侧链可以移动残基主链也可以稍微移动。配体准备和蛋白质准备与GLIDE刚性对接相同。搜索空间由默认参数定义内部盒子为10Å外部盒子为自动大小与配体类似大小并且中心位于天然配体的质心。与配体原子距离小于5Å的氨基酸残基被定义为柔性残基。对接过程按照标准协议进行生成多达20个姿态。总共有284个配体-受体对在PDBbind数据集中以及580个在MDT数据集中使用IFD模块成功对接。诱导拟合对接能够为PDBbind数据集中的更多配体-受体对成功输出姿态比GLIDE刚性对接多这表明该对接方法可以通过移动口袋残基来扩展口袋。
