1. 引言
近二十年来，各种病毒引起的传染性疾病频发，不但影响到人类的正常生活，也严重影响到了国家安全稳定。迄今，引起人畜共患病的冠状病毒有三种，即SARS-CoV、MERS-CoV和SARS-CoV-2。2020年初，新型冠状病毒SARS-CoV-2引起的新冠肺炎(COVID-19)席卷全球，给人类带来了巨大的生命和财产损失。目前，尚无针对COVID-19的特异性治疗方法。因此，迫切需要找到有效的抗病毒药物来对抗这种疾病。一种快速有效的药物发现方法是测试现有的抗病毒药物是否能有效地治疗其他相关的病毒感染。其中RNA聚合酶抑制剂瑞德西韦因其对SARS-CoV、MERS-CoV和埃博拉病毒的良好药效 [1] [2] 而被寄予厚望。最新报道表明：该药物在体外实验可以有效抑制SARS-CoV-2病毒 [3]，因此很快开展了临床用药试验。美国、欧洲和日本的联合科技团队发布的首份临床研究结果报告显示：53名新冠病毒感染者使用瑞德西韦后，其中68%的患者的呼吸苦难症状得以改善 [4]。但是，从创新靶向药物开发的科学严谨性与流程规范性和复杂性上考量，SARS-CoV-2并非完全等同于SARS-CoV、MERS-CoV及埃博拉病毒在内的任何一个病毒，因此不能仅仅寄希望于瑞德西弗这一种抗病毒药物，目前亟待发现更多可用于抑制SARS-CoV-2病毒RNA聚合酶的核苷类抑制剂。
针对SARS-CoV-2的药物筛选，北京望石智慧科技有限公司利用蛋白同源建模、计算化学虚拟筛选技术及基于人工智能的化合物筛选技术，发现二十余个核苷母核分子在末端三磷酸化后均可能对SARS-CoV-2病毒RNA聚合酶具有潜在的抑制作用，并向社会免费公开所有相关数据 [5]。本文基于望石智慧筛选出来的六种典型的RNA聚合酶抑制剂，借助量子化学计算将其几何构型和电子结构性质与瑞德西韦进行了对比分析，以期为筛选抗病毒药物提供一种新思路。
2. 计算方法
首先在M06-2X/6-31G (d)理论计算水平下对如表1所示的7种SARS-CoV-2抑制剂的初始结构进行几何优化，获得稳定分子构型。随后，在M06-2X/6-311 + G (d, p)水平下计算了7种药物分子的电子结构性质。此外，本文采用自洽反应场(Self-Consistent Reaction Field, SCRF)方法中的极化连续介质模型(Polarized Continuum Model, PCM) [6] [7] 模拟了水溶剂环境。本文所有计算都是采用量子化学计算程序包Gaussian 16 [8] 完成。分子构型和轨道图则用Gaussview软件 [9] 生成。采用Multiwfn程序包 [10] 结合VMD软件 [11] 获得静电势。分子对接是基于望石智慧公开的蛋白质模型(瑞德西韦与SARS-CoV-2 RNA聚合酶复合物结构)，使用SYBYL-X 2.0对接软件完成。
3. 结果与讨论
3.1. 几何构型
本文研究的7种RNA聚合酶核苷类抑制剂均来自于专利报道(见表1)。如图1所示，这7种RNA聚合酶抑制剂具有不同的核苷母核结构。与瑞德西弗不同，候选抑制剂化合物1~6并不包含磷酸骨架结构。
Figure 1. The geometric structures of seven SARS-CoV-2 inhibitors optimized at the M06-2X/6-31G (d) level
图1. 七种SARS-CoV-2抑制剂在M06-2X/6-31G (d)水平下优化的几何结构
Table 1. The Patent number and chemical formula for seven inhibitors
表1. 七种抑制剂对应的化学式及专利号
3.2. 电子结构
3.2.1. 偶极矩与电荷分布
药物分子的偶极矩可以影响药物分子和药物受体的相互作用。因此，偶极矩的研究分析对于药物研发具有重要的研究价值。偶极矩的大小可以反映分子极性的大小，偶极矩越大，分子极性越强，有利于药物分子在体内环境中的扩散与吸收，反之亦然。如表2所示，七种抑制剂的偶极矩数值从大到小顺序为：7.63 a.u. (化合物3) > 4.24 a.u. (化合物4) > 3.09 a.u. (化合物6) > 2.98 a.u. (瑞德西韦) > 2.93 a.u. (化合物2) > 2.32 a.u. (化合物1) > 1.65 a.u. (化合物5)。从整体上看，除了化合物3和化合物4有较大的偶极矩外，其他几种化合物的偶极矩均相对较小。由表2可知，化合物3和化合物4拥有比瑞德西韦更大的偶极矩，说明两种化合物较瑞德西韦而言更容易在体内环境中扩散与吸收；而化合物2和化合物6的偶极矩数值与瑞德西韦相差甚微，表明这两种化合物拥有和瑞德西韦相似的水溶性。
Table 2. The dipole moment (μ), vertical ionization energy (VIE), vertical electronic affinity (VEA), electronegativity (χ), hardness (η), softness (S), and electrophilic index (ω) of seven potential inhibitors at the M06-2X/6-311 + G (d, p) level
表2. 七种潜在抑制剂分子在M06-2X/6-311 + G (d, p)水平下计算得到的偶极矩(μ)、垂直电离能(VIE)、垂直电子亲和势(VEA)、电负性(χ)、硬度(η)、软度(S)和亲电指数(ω)
研究分子的电荷布居特征有助于详细分析分子的亲电反应和亲核反应的活性位点，以及分子间的相互作用。基于自然键轨道(Natural Bond Orbital, NBO)计算的自然布居分析(Natural Population Analysis, NPA)电荷是目前最为常用的电荷布居研究策略 [12]。因此，本文使用Gaussian 16程序包内置的NBO程序，在M06-2X/6-311 + G (d, p)水平下对七种新型冠状病毒抑制剂进行NBO计算，得到NPA电荷图，列于图2 (绿色区域表示正电荷，红色区域表示负电荷)。
如图2所示，七种SARS-CoV-2抑制剂的负电荷(红色区域)主要分布在O原子和N原子上。具体来说，六种化合物嘧啶环上相连的伯胺或仲胺结构的N原子拥有较大的负电荷，这一点与瑞德西韦相似。这些伯胺或仲胺结构的N原子上的孤对电子能够与嘧啶环产生p-π共轭，而且这种富电子的离域π键使之易与病毒RNA聚合酶的亲电部位产生相互作用。另外，由于O原子具有强电负性，对电子产生了吸引，使得这七种化合物中的呋喃环上的O原子、呋喃环上取代基侧链的甲醇基上的O原子均具有较大的负电荷。
六种SARS-CoV-2抑制剂的正电荷(绿色区域)主要集中在H原子和部分C原子上。根据图2观察到，这六种化合物呋喃环上羟基H原子的电荷与瑞德西韦同个部位H原子的电荷相差甚微。此外，还观察到这七种化合物呋喃环上羟基H原子的电荷大于其它部位H原子的电荷，推测可能是由于羟基中的O原子吸电子所致。当抑制剂进攻病毒RNA聚合酶时，这些带有较大正电荷的氢原子作为良好的亲电子体，可能更容易与RNA聚合酶中氨基酸残基上的N、O等原子结合，形成氢键。不同于其他五种化合物，由于化合物6呋喃环附近的O和N原子带有强电负性，对嘧啶环上的C原子有吸引作用，导致其嘧啶环上与O、N原子相连的C原子带有较大的正电荷。这一点同瑞德西韦是一样的，它们都是由于同时受到N原子和O原子吸引而产生较大的正电荷。但不同的是，瑞德西韦最大的正电荷分布在P原子上，而化合物6最大的正电荷处于C原子上。
3.2.2. 静电势分析
生物分子与药物之间的相互作用主要包括以静电相互作用主导的氢键、卤键等，而静电势(Electrostatic Potential, ESP)是与静电效应关系最密切的实空间函数，非常适合分析静电主导的弱相互作用 [13]。因此，对于静电势的分析有利于深入了解配合物分子结构中的电荷相关性质。由图3可知，七种化合物的嘧啶环、吡咯环上的N原子，呋喃环上的O原子均具有负的静电势分布，这与NPA电荷分析结果一致，表明上述的N原子和O原子均为给电子体。其中，化合物1、2、5、6的最大负电势分别为−62.90、−65.74、−47.00、−61.58 kcal/mol，主要集中在呋喃环的O原子中；化合物3和化合物4最大负电势分别为−91.41和−64.65 kcal/mol，分布在羰基附近；而瑞德西韦的最大负电位为−61.05 kcal/mol处于与P原子相连的O原子上。整体上，瑞德西韦和化合物6的最大负电位数值(−61.58 kcal/mol)仅相差0.53 kcal/mol，说明这两种抑制剂作用于病毒RNA聚合酶的位点可能是相同的。
Figure 2. NPA charge population diagram of SARS-CoV-2 inhibitors at the M06-2X/6-311 + G (d, p) level
图2. 七种SARS-CoV-2抑制剂在M06-2X/6-311 + G (d, p)水平下的NPA电荷分布图
Figure 3. The electrostatic potentials of SARS-CoV-2 inhibitors at the M06-2X/6-311 + G (d, p) level
图3. 七种SARS-CoV-2抑制剂在M06-2X/6-311 + G (d, p)水平下的ESP等势面图
包含瑞德西韦在内的7种抑制剂的静电势较正的区域主要集中在H原子和部分C原子上。其中化合物1、2、3、4、5的最大正电势分别为65.01、69.92、85.14、69.69、54.65 kcal/mol，主要分布在与呋喃环相连的羟基的H原子上。而化合物6最大正电位为62.78 kcal/mol，分布在与N原子相连的H原子上，这一点与瑞德西韦相同。因此我们推测当靠近病毒的RNA聚合酶时，这些与N原子相连的H原子可能更容易进攻RNA聚合酶氨基酸残基上的N、O原子等亲核位点，形成氢键。上述静电势的分布与NPA电荷分布的情况是一致的，为进一步筛选新冠病毒抑制剂提供了理论参考。
3.3. 前线分子轨道
根据前线分子轨道理论，能量最高的占据分子轨道(Highest Occupied Molecular Orbital, HOMO)与能量最低的未占据分子轨道(Lowest Unoccupied Molecular Orbital, LUMO)是决定体系发生化学反应的关键。因此，研究药物分子的前线轨道可以为探索作用机理和确定活性部位提供重要信息 [14]。本文得到瑞德西韦和其他六种抑制剂在水环境中的HOMO和LUMO轨道及其HOMO-LUMO能隙值如图4所示。
Figure 4. The frontier molecular orbitals of SARS-CoV-2 inhibitors at the M06-2X/6-311+G (d, p) level
图4. 七种SARS-CoV-2抑制剂在M06-2X/6-311 + G (d, p)水平下的前线分子轨道图
与瑞德西韦相似，六种化合物的HOMO轨道主要集中在吡咯环上。其中，化合物1、2、3、4可通过伯胺结构中N原子上的孤对电子与所连接的嘧啶环形成p-π共轭，因此它们的HOMO轨道上还分布于嘧啶环上。由于化合物5上的苯并氮杂䓬环对其HOMO轨道也有一定的贡献，所以该化合物的HOMO轨道与上述四种化合物略有差别。对于化合物6而言，其HOMO轨道的电子云虽未在嘧啶环上有分布，但苯环内的π电子对HOMO轨道有一定的贡献。因此，化合物6与瑞德西韦相比，其HOMO轨道电子云的分布并没有较大的区别。另外，化合物6的HOMO轨道还分布在仲胺结构中的N原子上。从图4可以看出，与嘧啶环相连的仲胺结构上的N原子能连接两侧杂环形成共轭，因此可以形成离域程度更高的分子轨道。
化合物2的LUMO轨道电子云虽然在嘧啶环也有分布，但主要还是集中在呋喃环上，而其他五种SARS-CoV-2抑制剂的LUMO轨道电子云主要集中在嘧啶环上，与瑞德西韦相似。其中，化合物1、4、5、6的LUMO电子云在嘧啶环上的分布几乎与瑞德西韦相同，主要分布在嘧啶环的π轨道之上。这就意味着当此类抑制剂靠近冠状病毒RNA聚合酶时，嘧啶环可能是与冠状病毒RNA聚合酶中富电子区域产生相互作用的主要部位。由图4可知，六种SARS-CoV-2抑制剂的能隙值均小于瑞德西韦。除了化合物1和化合物3外，其他四种新型冠状病毒抑制剂的能隙值在6.07~6.91 eV，数值较大。化合物6与瑞德西韦的能隙值分别为6.91 eV和7.18 eV，相差最小，表明化合物6可能具有和瑞德西韦相似的化学稳定性。
3.4. 概念密度泛函分析
众所知周，配合物的量子化学参数对其实际应用非常重要。因此，本文计算了七种潜在抑制剂的垂直电离能(Vertical Ionization Energy, VIE)、垂直电子亲和势(Vertical Electronic Affinity, VEA)以及采用概念密度泛函理论(Conceptual Density Functional Theory, CDFT)定义下的电负性(χ)、硬度(η)、软度(S)和亲电指数(ω)等理化性质 [15] [16] [17]，其计算公式如下：
(1)
(2)
(3)
(4)
由表2可以看出，化合物6的概念密度泛函参数数值与瑞德西韦最为相似，表明这两种抑制剂可能具有相似的生物活性。具体来说，化合物1和化合物3具有较低的垂直电离能，因此具有一定的还原能力。化合物2、4、5虽具有较大的垂直电离能，但硬度数值都较小，说明其稳定性不如瑞德西韦。化合物6的垂直电离能数值和硬度数值较大，稳定性较好。值得注意的是，化合物6和瑞德西韦电负性参数仅相差0.19 eV，软度仅相差0.01 eV−1。而且，化合物6和瑞德西韦的VEA分别为0.98 eV和0.90 eV，亲电指数分别为2.45 eV和2.47 eV，表明化合物6的亲电性能和瑞德西韦十分相似，因此也可能与SARS-CoV-2 RNA聚合酶富电区域结合，进而阻止病毒RNA合成，以抑制新冠状病毒复制。
3.5. 分子对接
分子对接是一种模拟小分子与受体之间静电、疏水、氢键、π-π堆积等相互作用的计算机模拟技术 [18] [19]。本文基于对接的准确性、打分函数、分子形态的相似性 [20] [21] [22] 等因素的考量，采用了SYBYL软件上的“Surflex-dock”功能模块对七种抑制剂与SARS-CoV-2 RNA聚合酶模型进行分子对接。对接时所用的三维聚合酶模型结构来源于望石智慧公开的蛋白质模型 [4]。在对接预处理过程中，抽离出原复合物中的瑞德西韦，然后对其他六种抑制剂进行分子对接研究，所得复合物结构如图5所示。
由图5可见，七种抑制剂都位于SARS-CoV-2 RNA聚合酶的活性口袋内。观测配体的形状，可发现七种抑制剂对接后的构象与望石智慧提供的晶体结构较为接近。由图5可知，六种抑制剂对接整体打分值(Total Score, TS)皆小于瑞德西韦。打分数值与配体–受体的解离平衡常数的负对数(−logKd)相关，得分越高说明药物分子与受体的结合越稳定 [23]。从图5可知化合物6的得分与瑞德西韦接近，评分较高，说明化合物6与SARS-CoV-2病毒RNA聚合酶之间的亲和力较大，且与瑞德西韦同该靶点结合能力相当。因此，可以推断化合物6也可以与新冠病毒的RNA聚合酶有效结合，从而抑制病毒的复制过程。
4. 结语
本论文基于密度泛函理论，针对瑞德西韦和其他六种潜在SARS-CoV-2抑制剂的电子结构、前线分子轨道及理化性质等进行了系统的比较。理论计算结果表明：化合物6的电荷分布、轨道性状、概念密度泛函参数以及与SARS-CoV-2 RNA聚合酶对接情况都与瑞德西韦最为接近，说明化合物6也能与新
Figure 5. Docking maps and total scores of seven inhibitors to SARS-CoV-2 RNA polymerase
图5. 七种抑制剂与SARS-CoV-2 RNA聚合酶模型分子对接结构图和打分值
冠病毒的RNA聚合酶进行特异性结合，从而抑制病毒的复制过程，进而起到治疗新冠肺炎的作用。本论文提出了一种通过量子化学计算获得药物分子的电子结构性质和概念密度泛函参数等，辅助进一步筛选临床前候选药物分子的新思路。
致谢
作者感谢北京望石智慧科技有限公司2020年2月5日向科研人员共享已公开专利中部分涵盖的RNA聚合酶核苷类抑制剂的可编辑结构化数据。
基金项目
福建省科技创新联合资金项目(2017Y9122)和福建省高校杰出青年科研人才计划(2018B032)。
NOTES
*通讯作者。