摘要：降低小分子内能对设计有效的药物至关重要。扭转角特征分析(Torsion profile)围绕可旋转键按一定步长旋转360°后每个构象的张力能，确定能量极小点的取向以及相互转换的旋转能垒，在药物化学研究中用来研究化合物的最佳构象。本文综述了Torsion Profile在提高结合亲和力、大环分子设计、优选母核以及提高靶标选择性等药物化学的用途。并演示了如何同Flare实现基于QM的Torsion profile。

肖高铿/2020-01-30

1. Torsion profile与势能面扫描

降低小分子内能(张力能)对设计有效的药物至关重要1。扭转能分析(Torsion energy analysis)或扭转角特征分析(Torsion profile)估算围绕单个可旋转键按一定步长旋转360°后每个构象的张力能，确定能量极小点的取向以及相互转换的旋转能垒，在药物化学研究中用来研究化合物的最佳构象2。

有两种方式进行扭转角特征分析：1）刚性旋转结构的一部分作并计算每个角度的能量来快速确定近似扭转角特征；2）限制旋转键的扭转角并允许结构的其余部分松弛以产生最小的能量。这两个方法就是我们之前已经详细地讨论过的一维刚性与柔性势能面扫描3。鉴于扭转角特征分析(Torsion profile)仅对感兴趣的一个两面角进行研究，所以这是一维势能面扫描。而我们通常所说的势能面扫描可以是多变量的（多个两面角，距离等）扫描。

柔性扫描比刚性扫描可以获得更详细的扭转角特征，所有的这些计算可以采用力场（MM）或量子化学（QM）的方法来完成，后者需要更多的计算资源，同时可用于研究结构细节，这些细节是经典MM模型做不到的。本文仅对基于QM方法的扭转角特征分析(Torsion profile)进行初步地总结。

2. 势能面扫描在药物化学中的应用

2.1 用于优化结合亲和力：辉瑞(PFIZER)公司JAK3选择性共价抑制剂PF-06651600的设计
Tofacitinib,1 2

Figure 1. Tofacitinib与先导化合物25

辉瑞(PFIZER)公司JAK3选择性共价抑制剂PF-06651600的最初想法来源于结构生物学，Tofacitinib与JAK3激酶结构域的复合物晶体结构为与CYS909共价结合的抑制剂设计提供了灵感4。围绕嘧啶胺N-C键旋转可能会使酰胺取代基紧邻CYS909，适当设计的亲电试剂可以捕获到亲核的半胱氨酸（图2）。在托法替尼母核上接上丙烯酰胺得到化合物2(图1)。化合物2对JAK3的抑制作用中等（底物浓度为1mM时，IC50=1930nM），而对JAK家族的选择性较差，这表明可逆的（Ki）蛋白相互作用是结合亲和力的主导力量5。 X射线晶体结构分析支持了这一假设：发现吡咯并嘧啶铰链结合片段在活性位点排列整齐，而丙烯酰胺取代的哌啶高度无序（数据未显示）。由于JAK3在其家族成员中对ATP的亲和力最高，所以靶标抑制的Ki成分一定不能成为活性的主要驱动力，而只能通过与CYS909的共价相互作用来提供合适的选择性。

势能面扫描在Torsion Profile中的应用-墨灵格的博客

Figure 2.Tofacitinib与JAK3的结合模式(PDB:3LXK)。复合物晶体结构提示：围绕N-C轴旋转可以将酰胺放置在离CYS-909很近的位置。

在化合物2基础上进一步探索了构效关系，其中化合物3与4相比，两者与蛋白的相互作用模式没有什么变化，但是两者对JAK3的抑制率有巨大差异：IC50分别为56与3840nM(底物浓度为1mM)。这可以用化合物3、4的两种旋转异构体的构象分布来解释。如图3所示，构象A为酰胺靠近Cys909巯基的共价结合构象，构象B是与Tofacitinib活性构象一致。构象能是用Gaussian软件在M06-2X/6-31G*理论水平下计算，结果如图3所示。化合物3的A构象比B构象能量更低，为优势构象，在气相条件下能量低了~0.6kcal/mol, 导致构象A与B的构象分布比例大约为70:30；而化合物4的A、B构象的能量差异仅为0.02kcal/mol，两种构象在气相中的分布大致相等（50:50）。也就是说化合物4的A构象比例比3的更低，所以化合物4的活性(IC50=3840nM)比化合物3的活性(IC50=56nM)更低。

势能面扫描在Torsion Profile中的应用-墨灵格的博客
ID R Conformer A (kcal/mol) Conformer B (kcal/mol) IC50 (nM)
3 H 0.0 0.6 56
4 Me 0.0 0.002 3840

Figure 3. N-Me对构象A的能量影响很大：化合物3的A与B构象比例为70:30,而化合物4的A与B构象比例大约为50:50。

2.2 大环化分子设计：勃林格殷格翰(Boehringer-Ingelheim，BI)公司BI-4020先导化合物的设计

在《文献重现 | 大环EGFR激酶抑制剂BI-4020的先导化合物5与6的发现与设计》一文中，详细地介绍了勃林格殷格翰(BI)公司在BI-4020的发现过程中如何使用QM Torsion Profile技术识别优势构象6。

一方面，如图4所示，作者对联芳基两面角进行柔性扫描，并与X-Ray的生物活性构象进行比较，以确认生物活性构象是否在能量上是否有优势。

势能面扫描在Torsion Profile中的应用-墨灵格的博客

Figure 4. 对不同取代的联芳进行量子化学扭转能分析

另一方面，作者用DFT计算来辅助判断是否有必要进行大环化锁定生物活性构象。化合物5(图5)代表了具有选择性但是活性强度不足的化合物。作者结合二维核磁(2D-NMR)实验与DFT计算的玻尔兹曼权重表明，化合物5具有两种优势构象：活性与非活性构象（图5）。因此假设如果可以将化合物5约束在活性构象则可以提高其对EGFR的结合亲和力。作者借鉴了lorlatinib的设计思路7：通过大环化来降低配体的熵变以提高化合物的生物活性。虽然原文没有提供所采用DFT计算的方法与基组，但是此类计算是非常常规的计算，比如在我们的ECD计算里就经常计算玻尔兹曼权重，相信很容易重现该计算结果。

势能面扫描在Torsion Profile中的应用-墨灵格的博客

Figure 5. 先导化合物5的活性与非活性构象

总的来说，在BI的这篇JMC文章中，Gaussian的势能面扫描(DFT)用来帮助识别优势构象，并用来决策是否有必要进行大环化以锁定生物活性构象实现提高化合物的结合亲和力。

2.3 优选母核：罗氏(Roche)PDE10A抑制剂的设计

势能面扫描在Torsion Profile中的应用-墨灵格的博客

Figure 6.（a）具有吡啶-2-酮、哒嗪-4-酮或吡嗪-2-酮母核的吡唑抑制剂对人类PDE10A的IC50值。（b）抑制剂26与人类PDE10A（PDB 5i2r）的晶体结构。蛋白质残基和抑制剂分别以灰色和青色显示。（c）对连接两个中心杂芳基的二面角τ进行了柔性扫描，用Gaussian 9841在B3LYP/cc-pVDZ理论水平进行计算。垂直虚线表示在3-吡唑基哒嗪-4-酮母核与PDE10A共晶结构中观察到τ的范围（28-40°）。

我们在《真实世界的分子设计》介绍过罗氏的PDE10A项目8，如图6所示，化合物25与26等其它化合物有一样的结合模式，但是25的活性低了很多倍。结合模式无法解释活性差异，怀疑配体的张力能是潜在原因。量子力学（QM）计算确实揭示了两个杂环之间扭转角τ的能量分布存在重大差异。26和27（红色和蓝色曲线）的最小能量与蛋白质结合构象中的扭转角很好地吻合，而25的优势构象却扭曲得多（品红色曲线，在τ≈90°时最小）。在τ≈30°处有相当大的配体张力能（2-3 kcal/mol），这将转化为结合亲和力下降约100倍。随后合成的化合物证实了这一观察结果。例如，吡嗪-2-酮化合物28的二面角极小值在τ= 0°处，但该化合物仍然有效，这是因为在τ= 30°时，配体张力受到限制（0.3 kcal/mol）。由于计算的配体张力能与IC50值之间具有良好的定量一致性，因此进行了类似的二面能量扫描，以便优先考虑虚拟杂芳基的骨架修饰，然后再开始合成（结果未显示）。

总的来说，在本研究中，Gaussian的势能面扫描解释了化合物活性的差异，并可以用来指导骨架修饰。

2.4 化靶标选择性: 选择性JAK1抑制剂的设计

PF-04965842(图7，化合物25)是辉瑞公司开发的正在临床研究中的选择性JAK1抑制剂，虽然JAK1与JAK2的结合位点相似，但该化合物对两者的IC50分别为29nM与803nM9。有多个因素与该选择性有关，磺酰胺末端的丙基构象对该选择性也有较大的贡献。

势能面扫描在Torsion Profile中的应用-墨灵格的博客

Figure 7.JAK1抑制剂PF-04965842的化学结构，是选择性JAK1抑制剂

由于完全不同P-loop空间取向，迫使配体的两个两面角S(O2)−C−C−C与N−S(O2)−C−C为了最大化与P-loop在形状与静电上的互补性而采取了完全不同结合构象。为了更好理解两面角的构象能，作者对N-甲基丙基磺酰胺在B3LYP/6-31++G理论水平进行了量子力学扭转角分析（见图8a与b）。

势能面扫描在Torsion Profile中的应用-墨灵格的博客

Figure 8.两面角 (a) S(O2)−C−C−C 与(b) N−S(O2)−C−C 在B3LYP/6-31++G水平计算的扭转角特征分析。绿色与蓝色箭头分别为化合物25与JAK1、JAK2结合的构象。

在JAK1结构里，配体的丙基侧链采取了最低能构象，比如两面角S(O2)−C−C−C~175°；而在JAK2结构里则采取了一个高能构象（ΔE=1.2kcal/mol）,对应两面角为~95°。相反，两面角N−S(O2)−C−C在JAK1里（两面角=126°）为高能构象（ΔE=3.2kcal/mol）,而在JAK2里为低能构象（两面角=172°）。该计算表明，在两个配体结合状态下，配体采取了能够与P-loop高效发生范德华力和静电相互作用的构象以支付内能损失。尽管计算出JAK1的总能量损失较高，但配体似乎能够通过与P环形成更有利的范德华接触来克服 N−S(O2)−C−C 扭转造成的构象能损失，从而提高了JAK1的选择性。 这与观察到的烷基越大选择性越高的趋势是一致的，因为这些烷基在JAK1中比在JAK2中能更好的与P-loop相互作用，即化合物5（14倍）、6（大于19倍）、25 （28倍）和27（66倍）。

第4小节以化合物25为例，演示了如何用Flare进行QM-Based Torsion Profiling。

2.5 Cdc7抑制剂: 解释一个氮原子引起的巨大活性差异

Pennington等人[10]以Cdc7抑制剂为例，说明单一氮原子取代对化合物活性的巨大影响。

势能面扫描在Torsion Profile中的应用-墨灵格的博客

Figure 9. 单一氮原子取代引起化合物对Cdc7的IC50产生300倍的差异

如图9所示，化合物（1）一个苯环上的C被N取代为化合物（2）后，对Cdc7的活性提高了300倍。如2.3小节Kuhn等人[8,11]所述，当结合模式不能解释活性差异时，分子的张力能往往是潜在的原因。在B3LYP/6-31G*理论水平对这两个化合物进行两面角扫描(torsion profiling)分析表明：化合物(1)的优势低能构象偏向于~150°而化合物(2)的优势低能构象偏好于~0°。这种构象偏好可能原因是：化合物(1)在两面角=0°时的H7和H6'发生立体冲突，以及化合物（2）两面角为180°时N7和N2'之间存在孤对电子排斥。

势能面扫描在Torsion Profile中的应用-墨灵格的博客

Figure 10. 化合物(1)在B3LYP/6-31G*理论水平的Torsion profile分析结果

势能面扫描在Torsion Profile中的应用-墨灵格的博客

Figure 11. 化合物(2)在B3LYP/6-31G*理论水平的Torsion profile分析结果

我们用Flare/Custom Force Field对这两个化合物进行同样级别的理论计算(具体方法见第4小节)，化合物1与2的两面角特征分析结果分别如图10与11所示。如图10蓝色曲线所示，化合物(1)有两个玻尔兹曼统计分布差不多(50:50)的构象，其两面角分别大约为30°与150°。如图11蓝色曲线所示，化合物(2)在两面角大约为0°与150°有两个低能构象，但是150°的那个构象比起在0°的那个构象的构象能高不止8kcal/mol，这意味着化合物(2)几乎以单一0°构象为主(玻尔兹曼统计分布~100%)。总的来说，化合物(1)比化合物(2)活性构象的玻尔兹曼统计分布低的多（分别大约为50%与100%），这导致了化合物(1)比(2)更不利于与Cdc7结合的构象分布以及更差的Cdc7抑制活性。

图10与11的两面角扫描分析还发现AMBER/GAFF2力场(橘黄色曲线)不能很好的描述化合物：GAFF2力场对化合物1、2没有区分能力，在大约70°附近找在一个最低能构象。橘黄色的GAFF2力场势能曲线与蓝色的QM势能曲线差异很大，这告诫我们：AMBER/GAFF2分子力场不能足够准确地描述这两个分子，这也是我们为什么用QM方法的原因。经过Flare的自定义力场优化后，新的力场（灰色）与QM计算的势能面拟合得相当不错。这说明，在用分子力场对该类分子进行结构优化、构象分析、分子动力学模拟与FEP结合自由能计算时，需要非常小心，有必要对力场重新参数化以便能更好地描述化合物。

2.6 磺酰胺相关的构象效应：用Torsion profile解释单一键级差异引起的活性差异12

势能面扫描在Torsion Profile中的应用-墨灵格的博客

图12. FXa抑制剂5与6在结构上仅存在一个单双键的差别

化合物5与6(图12)为FXa抑制剂，它们在化学结构上仅存在一个单双键的差别，但是对FXa的Ki却存在20倍的差异，分别为367与17nM。为了理解化合物间的活性差异，GSK的科学家们用X-ray解释了复合物结构：化合物5、6与FXa结合时各个原子重合的非常好(见视频1)，磺酰胺侧链的丙烯基与丙基完全重合；化合物5与6与FXa的结合模式完成相同（见视频2）；蛋白的构象也基本没有差异。总的来说，从相互作用模式的角度讲，不能解释化合物5与6为何具有20倍的活性差异。

视频1. 化合物5与6在FXa结合位点里的叠合效果(PBD code分别为2JH5与2JH6)

视频2. 化合物5与6与FXa的相互作用模式比较(PBD code分别为2JH5与2JH6)

由于化合物5与6结构极其相似，两者之间不太可能存在与溶剂化、去溶剂化效应相关的能量差异；也不存在与蛋白结合相关的自由度损失差异。或者说，即使有，这两种差异不太可能带来20倍的活性差异。

鉴于此，作者认为化合物5、6的活性差异仅可能是由于构象效应产生。于是分别对化合物5、6的两面角(C=C)–(S=O) 与(C–C)–(S=O)进行了基于QM的Torsion profile分析。计算采用Gaussian软件在B3LYP/6-31G*理论水平下进行，化合物5、6分别用简化的分子模型M1与M2代替，结果如图13、14所示。

势能面扫描在Torsion Profile中的应用-墨灵格的博客

图13. 化合物5的模型分子M1在B3LYP6-31G*理论水平下的Torsion profile,绿色圆圈表示化合物5与蛋白结合时的两面角角度

如图13所示，当C=C与S=O共平面时时是有利的能量最低构象(绿色箭头指示位置），但是这个构象不是化合物5与蛋白结合时的构象。复合物里化合物5的C=C发生185°旋转，能量比最低能构象高了3kcal/mol。

势能面扫描在Torsion Profile中的应用-墨灵格的博客

图14. 化合物6的模型分子M2在B3LYP6-31G*理论水平下的Torsion profile,绿色圆圈表示化合物6与蛋白结合时的两面角角度

模型分子M2的两面角扫描结果如图14所示，化合物6与FXa的结合构象正是M2的最低能构象，结合构象与最低能构象间不存在能量惩罚。

总的来说，基于QM的两面角扫描分析结果与实验观察到的生物学活性一致：化合物6比化合物5的活性要高非常多。化合物5与化合物6大约20倍的活性差异是因为构象效应带来的。

3. 用Gaussian进行Torsion Profiling

通常所说的Torsion profile是对一个两面角进行柔性的势能面扫描。更详细地流程与方法参见：Gaussian教程 | 势能面扫描3。我们在文献重现 | 大环EGFR激酶抑制剂BI-4020的先导化合物5与6的发现与设计一文中，详细说明了如何用Gaussian进行扭转角特征分析，并提供了冻结两面角进行并行加速计算的示例6。

为了加速计算，我们可以设计一个两步走的策略：1）先在PM7或PM6半经验理论水平进行柔性扫描；2）将得到构象分别进一步用高级的理论水平比如APFD/6-311+g（2d，p)进行固定两面角的结构优化。这个两步的计算策略可以充分利用高性能计算进行计算，可以在几个小时获得计算结果。

4. 用Flare V4进行Torsion Profiling

Flare V4的自定义力场参数模块可以用来进行两面角扫描，下面视频以前面2.4小节提到的JAK1抑制设计为例演示了如何在Flare里完成QM-Based Torsion Profile计算。Flare支持云计算，每个构象分配到一个16核心的节点，因为35个构象可以几乎同时完成。计算完毕，生成一套AMBER/GAFF力场参数与Excel格式的Torsion Profile报告。

5. 总结

基于QM的扭转角特征分析(Torsion profile)已经大量地在药物化学研究中用来分析最佳构象、辅助药物分子设计。本文以罗氏、勃林格殷格翰、辉瑞等几个公司的药物发现为例，说明Torsion profile如何帮助进行结构优化。用Gaussian软件通过使用SCAN关键字可以方便地对任意两面角进行扭转角特征分析；还也可以通过冻结两面角的方式将扭转角特征分析拆解成多个单独的两面角限制的优化计算子任务，从而加速扭转角特征分析研究。

6. 文献

  1. Sellers BD, James NC, Gobbi A. A Comparison of Quantum and Molecular Mechanical Methods to Estimate Strain Energy in Druglike Fragments. J Chem Inf Model. 2017;57(6):1265-1275. doi:10.1021/acs.jcim.6b00614
  2. Muegge I, Bergner A, Kriegl JM. Computer-aided drug design at Boehringer Ingelheim. J Comput Aided Mol Des. 2017;31(3):275-285. doi:10.1007/s10822-016-9975-3
  3. 陈宇. Gaussian教程 | 势能面扫描. 墨灵格的博客. 2017-12-17. http://blog.molcalx.com.cn/2017/12/17/gaussian-tutorial-potential-energy-surface-scan.html
  4. Chrencik, J. E.; Patny, A.; Leung, I. K.; Korniski, B.; Emmons, T. L.; Hall, T.; Weinberg, R. A.; Gormley, J. A.; Williams, J. M.; Day, J. E.; Hirsch, J. L.; Kiefer, J. R.; Leone, J. W.; Fischer, H. D.; Sommers, C. D.; Huang, H. C.; Jacobsen, E. J.; Tenbrink, R. E.; Tomasselli, A. G.; Benson, T. E. Structural and thermodynamic characterization of the TYK2 and JAK3 kinase domains in complex with CP-690550 and CMP-6. J. Mol. Biol. 2010, 400, 413−433.
  5. Thorarensen A, E. Dowty M, Ellen Banker M, et al. Design of a Janus Kinase 3 (JAK3) Specific Inhibitor 1-((2 S ,5 R )-5-((7 H -Pyrrolo[2,3- d ]pyrimidin-4-yl)amino)-2-methylpiperidin-1-yl)prop-2-en-1-one (PF-06651600) Allowing for the Interrogation of JAK3 Signaling in Humans. J Med Chem. 2017;60(5):1971-1993. doi:10.1021/acs.jmedchem.6b01694
  6. 肖高铿.文献重现 | 大环EGFR激酶抑制剂BI-4020的先导化合物5与6的发现与设计. 墨灵格的博客. 2020-01-18. http://blog.molcalx.com.cn/2020/01/18/egfr-bi-4020.html
  7. Johnson TW, Richardson PF, Bailey S, et al. Discovery of (10 R )-7-Amino-12-fluoro-2,10,16-trimethyl-15-oxo-10,15,16,17-tetrahydro- 2H -8,4-(metheno)pyrazolo[4,3- h ][2,5,11]-benzoxadiazacyclotetradecine-3-carbonitrile (PF-06463922), a Macrocyclic Inhibitor of Anaplastic Lymphoma Kinase (ALK) and c. J Med Chem. 2014;57(11):4720-4744. doi:10.1021/jm500261q
  8. 真实世界的分子设计. 墨灵格的博客. http://blog.molcalx.com.cn/2020/01/23/real-world-perspective-on-molecular-design.html
  9. Vazquez, M. L.; Kaila, N.; Strohbach, J. W.; Trzupek, J. D.; Brown, M. F.; Flanagan, M. E.; Mitton-Fry, M. J.; Johnson, T. A.; TenBrink, R. E.; Arnold, E. P.; Basak, A.; Heasley, S. E.; Kwon, S.; Langille, J.; Parikh, M. D.; Griffin, S. H.; Casavant, J. M.; Duclos, B. A.; Fenwick, A. E.; Harris, T. M.; Han, S.; Caspers, N.; Dowty, M. E.; Yang, X.; Banker, M. E.; Hegen, M.; Symanowicz, P. T.; Li, L.; Wang, L.; Lin, T. H.; Jussif, J.; Clark, J. D.; Telliez, J.-B.; Robinson, R. P.; Unwalla, R. Identification of N -{ Cis -3-[Methyl(7 H -Pyrrolo[2,3- d ]Pyrimidin-4-Yl)Amino]Cyclobutyl}propane-1-Sulfonamide (PF-04965842): A Selective JAK1 Clinical Candidate for the Treatment of Autoimmune Diseases. J. Med. Chem. 2018, 61 (3), 1130–1152. https://doi.org/10.1021/acs.jmedchem.7b01598.
  10. Pennington, L. D.; Moustakas, D. T. The Necessary Nitrogen Atom: A Versatile High-Impact Design Element for Multiparameter Optimization. J. Med. Chem. 2017, 60 (9), 3552–3579. https://doi.org/10.1021/acs.jmedchem.6b01807.
  11. Kuhn, B.; Guba, W.; Hert, J.; Banner, D.; Bissantz, C.; Ceccarelli, S.; Haap, W.; Körner, M.; Kuglstatter, A.; Lerner, C.; Mattei, P.; Neidhart, W.; Pinard, E.; Rudolph, M. G.; Schulz-Gasch, T.; Woltering, T.; Stahl, M.; G. Rudolph, M.; Schulz-Gasch, T.; Woltering, T.; Stahl, M.; Rudolph, M. G.; Schulz-Gasch, T.; Woltering, T.; Stahl, M. A Real-World Perspective on Molecular Design. J. Med. Chem. 2016, 59 (9), 4087–4102. https://doi.org/10.1021/acs.jmedchem.5b01875.
  12. Senger, S.; Chan, C.; Convery, M. A.; Hubbard, J. A.; Shah, G. P.; Watson, N. S.; Young, R. J. Sulfonamide-Related Conformational Effects and Their Importance in Structure-Based Design. Bioorg. Med. Chem. Lett. 2007, 17 (10), 2931–2934. https://doi.org/10.1016/j.bmcl.2007.02.034.