|
本文介绍一篇来自浙江大学侯廷军教授、康玉副教授和碳硅智慧联合发表在Chemical Science的论文《SDEGen: Learning to Evolve Molecular Conformations from Thermodynamic Noise for Conformation Generation》。该论文提出了一种将分子力学当中的随机动力学系统和深度学习当中的概率模型相结合的小分子三维构象生成模型：SDEGen。作者采用随机微分方程（Stochastic Differential Equation, SDE）模拟分子构象从热噪声分布到热平衡分布的过程，联合概率深度学习的最新DDIM（Denoising Diffusion Implicit Models）模型，不仅提高了模型生成构象的效率，并且在多项评测任务（包括构象生成质量、原子间距离分布和构象簇的热力学性质）上实现了精度的提升。如在构象生成质量上，其多样性指标优于传统方法22%，准确性指标优于传统方法40%；在热力学性质预测方面，将传统方法的精度提升了一个数量级，与量化计算的结果误差缩小至～2kJ/mol。除此之外，这篇文章还引入了晶体构象的比对实验和势能面分布实验，为构象生成任务的评测提供了更多维及更物理的视角。大量的实验表明，SDEGen不仅可以搜索到小分子晶体构象所在的势能面的势阱当中，还可以搜索到完整势能面上多个局域优势构象。同时，SDEGen模型计算效率极高，在分子对接、药效团识别、定量构效关系等药物设计任务中具有广泛的应用前景。
研究背景 分子的几何构象是一个分子内部所有原子的三维坐标所构成的集合，由于现实世界的热涨落，一个分子往往具有多个可能的构象。一个分子的构象族在统计物理上称之为系综。获得准确的分子三维构象对于热力学性质计算、3D-QSAR、药效团模型、分子对接等药物设计当中重要的任务息息相关。此外，在大部分应用当中，研究者不仅需要稳定的静态结构，还希望观察到更多的局部最优构象。所以利用物理方法进行分子构象生成是一个常规的手段。 传统的分子构象生成方法根据计算尺度可以分为从头算动力学（Ab Initio）、基于密度泛函理论（DFT）的动力学，以及基于拟合力场的分子动力学。从头算动力学和基于密度泛函的动力学由于计算量过大很难推广到大体系上，而基于力场的分子动力学虽然速度较前两者更快，但这种方法已经被证明是一种较为粗糙的近似。近些年来利用深度概率模型进行分子构象空间的探索已经取得了重大进展。目前的研究表明，在原子数小于9个的小分子数据集GEOM-QM9上，基于深度学习的方法大多比传统的构象生成方法ETKDG方法有更好的性能，但是在与药物设计相关的GEOM-Drugs数据集上表现仍不够优秀，有很大的提升空间。基于此，受到生成模型最新进展和随机动力学系统的启发，本文作者开发了SDEGen，一种基于随机微分方程（SDE）的深度生成模型。与传统的回归模型不同，SDEGen不仅可以生成一个能量最优的构象，还可以生成一系列与真实热力学环境相符的局部最优构象。 材料与方法 数据准备 本文使用了三个数据集用于评测SDEGen方法，即GEOM-QM9、GEOM-Drugs和ISO17。GEOM数据集当中的构象是由量子化学软件CREST计算的；GEOM-QM9的训练集当中包含40k个分子的200k个构象，测试集中有200个分子及其对应的22408个构象；GEOM-Drugs训练集当中包含40k个分子的200k个构象，测试集中有200个分子及其对应的14324个构象。ISO17数据集中的分子是从QM9数据集中最大的异构体集合中随机抽取，均有固定的化学式（C7O2H10）。其构象计算所利用的从头算动力学比半经验方法精度更高，并且由于采样的构象分布更为平滑，所以ISO17数据集用于评估原子间距离分布的评估任务。ISO17训练集包含167个分子的357621个构象，测试集包含30个分子的73071个构象。 物理内涵  图1：该图形象化地展示了SDEGen的物理内涵，即基于随机微分方程构建的随机动力学系统，将一个超球（随机噪声分布所代表的流形）上的相点映射到复杂的构象流形上面。 一个分子所处于的相空间维数是3N维，其中N为一个分子所含的原子数目。由于能量约束的存在，分子构象并不是均匀分布在相空间内，而是在高维相空间内张成一个低维流形。SDEGen的初始采样相当于在高维相空间随机采样，然后经过随机微分方程所表示的动力学系统演化到原始数据分布的低维流形上，形成热力学稳定的分子构象。这种相空间上的运动可以理解为力场指导下的运动，而该随机微分方程可以近似理解为一种动力学驱动的方法；也可以理解这个过程是一种从3N维空间的一个D维超球到D维复杂流形上的微分同胚映射。举一个形象的例子，如果将一团粒子置于水中，那么粒子一开始服从的分布是一个随机的噪声分布；随着水分子与粒子发生不断地随机碰撞，最后这一团粒子会演化到一个稳定的构象分布上去。 模型架构 SDEGen旨在学习给定数据分布扰动为随机噪声的过程，通过反转这个过程，随机噪声可以被平滑地退化为原始数据，这也就是SDEGen的采样步骤。具体而言，对x的分布加噪过程可以用如下的随机微分方程描述：  其中f(·,t)是x(t)的漂移系数，g(·)是x(t)的扩散系数，w是一个布朗运动。这个公式表示对原始数据分布不断添加噪声直到其演化为一个不含原始分布信息的高斯分布。其逆过程也是一个随机微分方程：  其中bar{w}是一个标准的从时间T演化到时间0的维纳过程。一旦确定了每一个时间步下的边缘分布的导数▽xlogPt(x)，即可数值求解出上述的逆向随机微分方程进行采样。所以最后的目标是训练一个网络sθ用于近似▽xlogPt(x)。 为了估计▽xlogPt(x)，需要训练一个时间相关的网络sθ(x,t)：  其中u是[0,T]是的均匀分布，pt(x)是x(t)的概率密度函数，p0t(x(t) | x(0))是x(0)到x(t)的转移核，λ(t)∈R > 0是一个权重函数。当网络sθ(x,t)训练好时，通过Euler-Maruyama或者Predictor-Corrector求解器可以数值解逆向随机微分方程(2)，实现对分子构象的采样。  图2：SDEGen的框架：在训练时，给定图G和构象R，1）从[0,1]均匀分布中抽取时间，利用高斯随机特征将时间信息编码到模型中；2）映射分子的边（E）和原子（A）特征，形成相应的嵌入（浅蓝色和稍深的蓝色箭头）；3）利用GNN模型将图的结构编码到模型中，通过去噪方式训练SDEGen网络（深蓝色箭头）。该程序相当于学习分子在给定时间内的随机动力学系统中的演化状态。 如图2所示，在SDEGen中，随机选择一个时间步长后，将该时间步下的热力学噪声加入到初始原子间距离中形成扰动的原子间距离~dt, 接着将其通过一个MLP层映射为距离特征。边的类型信息经过一个Embedding层之后与上述的距离特征点乘之后，获得与键类型相关的距离特征(~dt|E)。另一方面，分子内原子的类型也经过embedding层和GNN层获得边信息的特征，与(~dt|E)相加形成与分子键类型，原子类型相关的距离特征(~dt|G)。最后经过一个MLP层映射为一维的向量，计算得到带有原始噪声的损失。这个过程在不同的分子和[0,1]时间步内重复多次直至收敛。基于经验力场的能量优化被嵌入到模型的最后一步以微调去噪，以得到最后的分子构象。 实验设置 生成分子质量：在这项任务中，作者对测试集中的每个分子图生成两倍于其对标构象数量（量化计算所得）的构象，并计算生成的构象对标准构象之间的COV和MAT指标。其中COV和MAT是两个由均方根偏差（RMSD）导出的指标，RMSD是动力学模拟中分析中两个构象之间差异的衡量标准。  其中n是分子中非氢原子的数量，Φ是将两个构象对齐的函数，COV和MAT的指标定义如下：  其中Sg和Sr是生成的和参照的分子构象系综，δ是给定的RMSD阈值。一般来说，较高的COV分数代表生成构象多样性高，而较低的MAT分数代表生成的构象与量化计算所得的构象更为接近。 原子间距离分布：由于共价键的长度不足以完全复刻构象的三维几何信息，所以作者在第二个任务中测量的原子间距离包括键长（1-2连接）和1-3及1-4连接。这种考虑相当于衡量原子之间的直接相互作用，正是这种作用促使原子从随机分布弛豫到现实世界的热力学分布。在这项任务中，作者对每个测试分子采样1000个构象作为伪轨迹，然后用高斯核计算伪轨迹与从头算轨迹分布之间的最大均值差异MMD。具体而言，对于测试集中的每个分子，作者计算了原子间所有距离p(d|G)、成对距离p(dij,duv|G)和单个距离p(dij|G)的MMD。 热力学性质预测：一个系统的宏观热力学性质是所有可能的微观状态的加权和。对于测试集中的每个分子，作者利用PyScf软件（DFT（M06-2X/def2-TZVPP））来计算每个生成的构象和基准构象的电子能。然后统计生成构象簇和基准构象簇的热力学性质的平均绝对误差MAE。评测的指标包括平均能量、最低能量、最高能量、平均HOMO-LUMO间隙、最小间隙、最大间隙。 晶体构象和其他稳定构象：由实验所得的晶体构象是药物设计中的金标准。在大多数情况下，晶体构象处于自由能面上的局部极小值附近。从实际应用角度来讲，好的模型生成的构象应该覆盖自由能面上的多个局部极小值，以便将实验确定的自然晶体构象包括在生成的构象集合中。为了检验模型是否可以生成类-晶体构象，作者将Platinum数据集作为外部测试集。Platinum数据集包含4626个结构，从PDB当中的347k个共晶配体结构当中过滤得到。在GEOM-Drugs上训练的SDEGen模型生成Platinum数据集当中的分子构象，与原始晶体构象进行比较。为了直观地理解模型生成的构象集合和晶体构象的关系，作者从PDB的配体库中随机抽取一个只有两个可旋转键的分子做势能面扫描，然后将模型生成的构象和晶体构象绘制于扫描得到的势能面上。除此之外，为了检验模型捕捉势能面上的极小值的能力，作者还进行了一系列的案例研究。前半部分的案例研究为GEOM-QM9当中的12个具有两个可旋转键的分子的二面角扫描结果；后半部分的案例研究作者为GEOM-Drugs当中2个具有12个可旋转键分子的势能面。对于具有12个可选择键的复杂分布，作者利用半经验精度方法对分子做元动力学采样，之后采用TICA方法获得两个降维坐标，绘制出相应的势能面。 结果与讨论 构象生成质量：作者在该实验当中选取了5种AI模型（ConfGF、DMCG、DGSM、CGCF和ConfGFDist）以及1种经典构象生成算法（RDKit 中的ETKDG）作为基线。SDEGen在两个数据集上都显示出优异的性能，并在GEOM-Drugs上取得了SOTA的结果。与端到端的ConfGF相比，SDEGen在两个数据集的力场优化条件下在多样性上超出~5%，准确度上超出~10%。此外，与随机添加非结合边来增强长程相互作用的DGSM相比，SDEGen只通过加入1-2和1-3相互作用就取得了更好的COV和MAT分数，进一步证明了SDEGen卓越的构象生成能力。 表1：不同方法（联合分子力场弛豫）在GEOM-QM9和GEOM-Drugs数据集上的COV和MAT结果（QM9的阈值为0.5埃，Drugs的阈值为1.25埃）。  原子间距离分布：原子间距离不仅包含有共价键原子间的键长，还包括非共价键相互作用，即1-3键角和1-4二面角相互作用。如表2所示，尽管SDEGen的三维重建过程会降低对距离分布的估计精度，但其在所有的距离分布指标上仍然取得了令人印象深刻的结果。具体而言，SDEGen 在Single-Median、Pair-Median、All-Median和Pair-Mean指标上均优于SOTA方法ConfGF，在Single-Mean（0.3943 vs. 0.3684）和All-Mean（0.6249 vs. 0.6091）指标上与ConfGF相差不大。值得注意的是，尽管DMCG在第一个构象生成质量任务中表现优异，但在这项任务中的表现却不及SDEGen和ConfGF模型，可能DMCG在一定程度上存在过拟合现象。SDEGen与其他受DG启发的方法（即GraphDG和CGCF）的比较表明，SDEGen可以学习不同类型原子和化学键的平滑距离分布。对原子间距离分布的准确估计从另一个角度证明了SDEGen生成构象的合理性。 表2：与对标的量化构象集合相比，不同方法的原子间距离分布MMD的平均值和中位数。  热力学性质预测：在这项任务中，每个构象对应于一个微观的热力学状态，将这些微观状态依据玻尔兹曼权重集合起来就代表一个热力学系统。一般而言，在热力学计算中考虑的微观状态越具有代表性，相应的计算结果就越准确。作者评估了由SDEGen及其他基线模型（即RDKit、CGCF、ConfGF和DMCG）生成的构象集合的宏观热力学性质。计算结果显示（见原文），由SDEGen生成的构象的集合属性最接近（~2kJ/mol）量子化学计算得到的结果。在所有的结果中，SDEGen在这个任务上的表现大幅优于经典方法RDKit（~50kJ/mol），表明了SDEGen所代表的随机动力学系统对分子构象的热力学演化实现了量子精度级别的建模。 晶体构象和其他稳定构象：如图3B所示，SDEGen生成的构象以超过85%的概率命中Platinum数据集中的晶体结构（RMSD阈值1.5 Å）。此外，对于晶体构象的COV（覆盖率）并没有随着生成的构象数量的增加而明显增加，表明了模型良好的鲁棒性。作者进一步对一个有晶体构象的分子进行二面角扫描，如图3D所示。势能面上橙色的点代表SDEGen生成的构象，黄色点代表晶体构象。可以发现，由SDEGen生成的50个构象都落入了势能面的势阱当中，覆盖了几乎所有的局部最小值；同时，晶体构象也落入了SDEGen生成的构象所覆盖的势阱中。此外，作者还探索了含有12个可旋转键的类药分子体系。结果显示，SDEGen采样得到的样本点都集中在量子化学计算得到的优势构象附近；与RDKit生成的点相比，SDEGen生成的构象在势能面上的分布更接近于量子化学计算的结果。  图3：（A）生成的构象与代表性分子的晶体构象的比较图；（B）晶体构象的COV与RMSD阈值之间的关系，不同颜色的线代表不同的构象生成数量；（C）和（D）分别是利用DFT扫描分子的势能面的3D和2D表示。其中黄色的点代表晶体构象，橙色的点代表生成的构象。 结论 在本研究中，作者结合物理原理和前沿的深度生成模型，对分子构象的随机动力学演化进行建模，学到了原子从随机热噪声分布开始，最终弛豫到能量最优附近的过程。SDEGen在生成构象质量、原子间距离分布和热力学特性预测等三个实验上超过了大多数其他的基于AI的构象生成模型。此外，SDEGen模型与SOTA方法ConfGF相比，采样速度加速10倍左右，这对于在实际应用中的大规模虚拟筛选至关重要。与晶体结构相比，作者展示了SDEGen能够以80%的概率快速命中Platinum数据集中的小分子晶体结构中的构象。此外，势能面扫描结果显示，SDEGen生成的构象不仅接近于量化采样下的稳定构象，还可以生成完整势能面上多个局域优势构象。 参考资料 Zhang H, Li S, Zhang J, et al. SDEGen: Learning to Evolve Molecular Conformations from Thermodynamic Noise for Conformation Generation[J]. Chemical Science, 2023. 数据及代码链接： https://github.com/HaotianZhangAI4Science/SDEGen |
|
|
