【分子构象生成】mETKDG-几何方法辅助的环状化合物的构象生成

苏黎世联邦理工大学物理化学院的Sereina Riniker 教授课题组于2020年在Journal of Chemical Information and Modeling上发表了题为“Improving Conformer Generation for Small Rings and Macrocycles Based on Distance Geometry and Experimental Torsional-Angle Preferences”的工作，讨论了他们在环状化合物构象生成方向的研究成果，利用现有知识对环状化合物的认识来构建特殊的几何关系，从而达到更准确采样的目的。

——背景介绍——

生成分子三维结构的能力是药物发现的一个重要助力，因其可能是3D-QSAR，分子对接或分子动力学研究的起点。目前市面上已经有很多构象生成器被开发出来供研究者使用，包括商业产品（一些提供免费的学术许可证）如Omega, ConfGen, MacroModel and Prime, CSD Conformer Generator, LowModMD, BRIKARD, ForceGen 3D, TCG 和 Conformator, 当然还有一些开源的软件工具比如说 CONFECT, BALLOON和RDkit中的构象生成器ETKDG。

构象生成器ETKDG是一种随机搜索方法（也是本文作者开发的），它利用距离几何和从实验晶体结构中得到的知识。它已被证明可以为无环的柔性分子产生良好的构象。本文的工作建立在ETKDG的基础上，改进了含有小的或大的脂肪族（即非芳香族）环的分子的构象生成。首先，我们设计了额外的扭转角势来描述小的脂肪族环，并调整先前开发的无环键势，以促进大环的采样。然而，由于大环的自由度较大，取样的构象空间比小分子要宽得多，这给构象生成器带来了挑战。因此，我们引入不同的启发式来限制宏循环的搜索空间，并使抽样偏向于更实验相关的结构。具体地说，我们展示了使用椭圆几何和可定制的库仑相互作用作为启发式。改进的mETKDG的性能在不同的大环和环肽的测试集上得到了证明。

——方法——

Ø  ETKDG的做法流程：

RDkit中最早的分子构建方法名为DG，虽然DG表现得很好，但生成的构象可能会扭曲芳香环和sp2中心，以及二面角值超出了实验晶体结构中观察到的范围，所以作者在之前ETKDG工作中进一步针对二面角进行了的优化。（还引入了一些规则性的约束来限制芳香环在同一平面内，增加sp2杂化的翻转情况，并且约束了三键体系的键角）

图一：ETKDG的流程图

1. 首先先要居于拓扑结构确定一系列的距离界矩阵：包括1-2（相当于距离），1-3（角度由角度决定），1-4（由二面角决定），1-5（距离）；

2.  然后利用结构搜索算法来找到匹配的二面角样式：首先作者根据CSD数据库中的大部分化合物，统计出了387种分布相互独立的二面角种类，然后将结构中的二面角与之进行匹配，得到二面角信息的先验知识，用来更新界矩阵；

3. 然后就是根据搜索到的结构信息重新设定一系列的界矩阵；

4. 由于生成的界矩阵可能会存在空间上的不连续的问题，作者利用三角不等式来平滑化界矩阵；

5. 根据这些相互作用界矩阵，随机生成符合这一系列约束的距离矩阵（笔者认为构象生成的随机多样性体现在这一步）；

6. 对界矩阵做降维处理，找到排名前三的特征向量做投影，作为三维坐标输入，使每个维度上的方差尽可能的大；

7. 对生成的原子坐标根据拓扑关系进行优化。

当然步骤6可以替换为采用随机生成的坐标，而不用PCA，这样增加了随即多样性，但是在后续的测评中，准确程度有所下降。

Ø  mETKDG的新增内容

作者在后续的研究中发现，在环状化合物体系中在RDKit中边界矩阵的生成过程中，根据所涉及的环系原子的数量，设置1,4界矩阵时存在不同分布，即环系原子数以8为分界线，分界线及其以下二面角的分布会存在偏向性，分界线以上二面角的分布就与线性分子基本无异，这一点可能是由于环张力随着环大小的改变引起的。但是在ETKDG中，这些特异性的分布并没有被体现。

因此作者针对环状分子分两种情况进行探究，即原子数小于8和原子数大于8的情况，我们暂时分别代称为小环、大环。

首先作者采用之前的做法，重新收集了CSD数据库中有关环状化合物的数据，形成了环状化合物中有偏的105种二面角的分布（分为大环和小环分别进行统计）。

图二：环状化合物二面角的分布（截取部分为例）粉色代表小环，蓝色代表大环

 

图二也可以直观体现出大环和小环分布的差异，对于小环体系，作者直接将有偏的二面角分布加入ETKDG算法中，就拿到了很好的结果。

但是对于大环体系，因其分子柔性大，所以界矩阵的分布与线型分子并没有明显的差异，但是环状化合物的结构稳定性往往比线型分子要好一些，这是源于环状分子特殊的几何拓扑结构，但是这一点要怎么体现在构象生成之中呢？

作者抓住大环分子的几何特点，启发式地想出了一个解决方案：

1. 对于大环，这里作者仍然使用了非环体系的距离界矩阵的参数；

2. 将三维的结构投影至二维形成二维图形，然后使用椭圆的几何结构来描述和约束坐标的生成：

          

图3：大环分子（左）和其经过PCA降维后拟合出的椭圆（右）

利用周长和离心率参数来控制环状结构的“圆度”，给定第一个原子的坐标为椭圆的一个顶点，剩下的原子按照顺序排列，利用非线性变换拟合出一个椭圆方程，然后利用这些距离约束信息来修正界矩阵的数据，再输入模型。

图四：不同离心率的大环化合e1=0.08, e2=0.65,e3=0.92

3. 加入CPCI（Custom Pairwise Coulombic Interactions），由于之前的小环体系中的电荷相互作用不是很重要（主要为脂肪环），所以没有考虑真正的静电相互作用。这里由于大环体系中会存在较多的杂原子之间的极性相互作用（例如环肽中的分子内氢键），所以显式地加入原子两两之间的静电相互作用来优化结构。

此外，模型会考虑到大环化合物的长度来生成不同数目的构象，因为环系原子越多，分子柔性越大，需要的构象就更多。

 

——结果与分析——

对于小环分子，作者利用不同的阶段的模型对于生成的构象与晶体构象（CSD数据集）的rRMSD和rTFD进行了分析：

图五：A：不同生成构象topn的rTFD和rRMSD对比，B: 不同方法生成六元环分子的二面角分布

从图五不难看出，ETKDG+sr（即ETKDG加入小环二面角矫正）的生成效果最好。

这里重点看一下大环体系的效果，对于大环体系，作者首先构建了三种测试集mc-OE-set（CSD+BIRD+PDB的大环结构）、mc-D3R-set（BACE数据库中的所有大环ligand）、mc-PEP-set（decapeptide 和 cyclosporine A (CsA) 50 ），然后在这三种测试集上分别测试了mETKDG的性能：

图六：mETKDG在mc-OE-set上生成任务的表现，RC代表不加入椭圆的几何约束直接使用随机坐标，Eccen表示加入椭圆几何条件，flatness表示结构骨架的平滑度（越接近1越好，因为测试集本身的平滑度就很接近1）

图七：mETKDG在mc-D3R-set上生成任务的表现

对于环肽测试集，作者做了更为细致的测评，还分析了后续进行力场优化后的结构rRMSD的分布情况：

图八：mETKDG在环肽数据集上的表现，橙色是只用mETKDG加入椭圆约束和CPCI生成构象的rRMSD分布；蓝色是mETKDG加入椭圆约束和CPCI生成构象后经过OpenForceField优化后的结构rRMSD的分布；绿色是mETKDG只加入加入椭圆约束生成构象的rRMSD分布

——总结——

由于大环化合物，尤其是环肽的结构数据较少，因此很难采取目前热门的深度学习的方法对其进行建模和构象生成，本文提供了一种简便易行的方法高效地生成较为准确的大环化合物构象。但是我们可以从测试结果中看出来（作者在也在最后提到），生成构象的准确度和一些更加细致的氢键相互作用有些依赖于后续的力场优化，这也是后续这方面工作的一个上升空间。

 

Reference

Wang S , Witek J , Landrum G A , et al. Improving Conformer Generation for Small Rings and Macrocycles Based on Distance Geometry and Experimental Torsional-Angle Preferences[J]. 2020.

Riniker S , Landrum G . Better Informed Distance Geometry: Using What We Know To Improve Conformation Generation[J]. Journal of Chemical Information & Modeling, 2015:2562-74.