第一性原理计算2019-nCoV病毒分子3CL水解酶结构

新型冠状病毒肺炎2019-nCoV3CL水解酶的分子结构已经被中国科学家公布，这对病毒药物开发有着重要的作用。3CL水解酶是病毒的一个关键蛋白质，病毒需要利用它来复制RNA。如果把3CL水解酶的功能抑制住，就可以大概率阻隔病毒的复制，为抗病毒药物研发指明方向。一个蛋白质要起作用，需要与其它蛋白质、DNA或RNA 相结合。类似钥匙和锁的关系，蛋白质的结合部分会有一个口袋（pocket），如果在这个pocket 上放一个小分子，那么其他蛋白质就不可能再与它结合，这个蛋白质就失效了。如果一个小分子只与这个病毒蛋白质结合，而与我们身体的其它蛋白质结合效果不强，那么它的毒性就很小，这正是靶向药物的原理。

所以要研制小分子靶向药物，关键之一就是要寻找合适的小分子，看它与靶向蛋白质之间的结合能，这其实是一个物理问题。在现有的药物设计与计算中，经常使用的是经典力场方法。经典力场方法假设每个分子键有一个弹簧一样的力，原子与原子之间还有点电荷的相互作用力等等。经典力场的优点是计算速度快，缺点是精度不高，往往达不到几个kcal/mol 的生化精度要求。经典力场方法很难描述电荷的极化以及它对结合能的影响。

近年来基于量子力学的第一性原理计算有了很大的发展，密度泛函的方法已在材料模拟领域广泛应用。使用新的SCAN+rVV10泛函，已经可以用来计算生物分子电子结构以及它的范德瓦尔力。新的计算算法、计算软件以及新的GPU计算芯片也使得生物分子尺度规模的计算变为现实。基于量子力学的第一性原理计算比经典力场方法的计算精确很多。那么问题来了，我们可不可以利用第一性原理计算，来帮助设计药物小分子呢？将第一性原理计算与经典力场方法计算相结合，取长补短，用来解决一些现在无法解决的问题。

我们以2019-nCoV病毒的水解酶为例，看看我们如何来计算这个体系。非常感谢上海科技大学的饶子和/杨海涛课题组，他们公布了2019-nCoV 冠状病毒3CL水解酶(Mpro)的高分辨率的晶体结构图，我们从他们那里要来了PDB file(相当于xyz file)。首先需要把H 原子补上，这个蛋白质晶体结构在VESTA 上看起来是这样的：

如上图所示，结构非常大，一个晶胞有上万个原子。这个结构使用传统的计算方法是无法计算的，好在真正的drug和它的pocket并不大。晶体结构有四重对称性并由四个蛋白分子组成。如果只关注其中一个蛋白质，如下图：

只取上图中方框里的原子（包括pocket 和drug)，那么我们的系统只有大约800个原子，如果范围取大一点，能有1600个原子（这些包括切割表面的passivation atoms)。这就是我们可以计算的尺度了。在上面的结构中ligand结构如下：

首先算一下结合能，就是pocket与药物分子之间的吸收能。这个只需要计算整个体系的能量、pocket的能量和药物分子的能量，然后相减就可以得到。这里我们是用平面波方法计算的，使用的是NCPP中的SG15赝势，截断能为50Ryd。

有意思现象的是charge transfer due tobinding，极化产生的电子电荷重新分布，这正是经典力场计算不了的部分。我们本次计算出来的结合能是负的，也就是说pocket与ligand不愿意相互结合。这可能有很多原因，比如我们没有考虑水的作用，这个在一定程度上可以用DFT中solvent model模块来解决。另一方面，我们没有考虑configuration entropy，这个必须通过经典力场thermodynamicintegration来计算。最后一个可能是我们需要对原子进行结构优化（relax)，原子位置relax花的时间就更长了，大约需要十二个小时左右。下面是relax过程中能量的变化。我们可以看到能量变化非常大。

在下图中我们可以看到，relax带来的原子移动还是非常大的。在relax的过程中，为了避免pocket 散架，我们固定了切割表面的原子。

最后，如果体系再大一些，DFT还能算吗？我们用了1600 原子的体系，最后的结果如下图。这只是一个自洽计算，DFT计算已经有点吃力了，但还能计算，大约需要三个小时。

需要特别说明的是，我们上面做的只是一些测试，并没有说哪一个ligand更好。我们想知道第一性原理计算对药物设计是否有用。通过测试，我们认为它已经可以起到一定的作用，起码可以起到一些辅助设计的作用来加速抗病毒药物的筛选。在这个全世界抵抗2019-nCoV病毒的非常时期，我们也希望在科研上贡献一份力量。我们自己并不是药物研究者，如果有药物研究者认为我们的计算工作有用，欢迎与我们联系（www.lonxun.com），我们愿意同大家合作工作加速抗病毒药物的研发工作。

欢迎大家交流转发。