# 辉瑞新冠病毒新药COVID19_restrained_Gromacs_ABFE **Repository Path**: wangbinyuhao/abf_restriction ## Basic Information - **Project Name**: 辉瑞新冠病毒新药COVID19_restrained_Gromacs_ABFE - **Description**: 7SI9 是新上市药物辉瑞新冠病毒特效药的常温下的晶体结构。这个项目尝试使用Gromacs对该结构进行绝对结合自由能（ABFE）计算，以希望得到一个结合自由能的数值与实验值进行比对。 - **Primary Language**: Unknown - **License**: MIT - **Default Branch**: master - **Homepage**: None - **GVP Project**: No ## Statistics - **Stars**: 0 - **Forks**: 2 - **Created**: 2022-05-13 - **Last Updated**: 2022-05-13 ## Categories & Tags **Categories**: Uncategorized **Tags**: None ## README 绝对结合自由能（Absolute binding free energy， ABFE） 是指小分子药物与靶点蛋白结合后相对于结合之前的能量变化，是蛋白靶点，溶剂包括离子，小分子三者之间能量的变化。近年来，ABFE的关注度越来越高，因为此方法不但可以比较骨架类似的小分子对靶点的结合强弱，也可用来比较骨架不同的小分子。在所有的ABFE计算策略中，基于炼金术化学转变的方法是基本方法之一，也是本文讨论的范畴。通过缓慢地，顺序的关闭共晶结构中小分子和口袋的相互作用，直到完全“消除”小分子，通过充分的分子动力学模拟采样，就可以反推出小分子完全存在时的结合自由能。 请注意区分，我们计算的是一个过程的结合自由能，这本质上仍是一个自由能的差值，而不是蛋白和小分子结合后的那个“状态”的自由能，所谓”绝对结合自由能“的叫法可能会带来歧义，特别是对专业物理学背景的人而言，因为蛋白和小分子结合后的那个体系的绝对自由能可能是完全难以计算的，理论上几乎都难以计算，这是天文数量级别的计算量。即使物理学，也只能计算有限个理想气体模型的”绝对自由能“，绝大多数人们所讨论的都是“自由能”的差。在药物设计领域，人们经常还会提到一个“相对结合自由能”的时候，所谓绝对结合自由能可能容易被误解。 在缓慢消除小分子作用力的过程中，小分子药物与蛋白口袋的结合不可避免的变弱，以至于难以维持一个相对固定的姿态，在我们施加的分子动力学所处的温度和压强条件下甚至会漂移到溶剂中导致无法近一步采样。为了解决这个问题，我们需要采用一个所谓的“分子间的约束” （intermolecular restrain）来维持小分子和口袋的相对位置，因此引入了“一个距离，两个角度，三个二面角”的约束参数，人为强迫小分子在全过程的模拟中都要保持这六个参数。这解决了小分子的漂移问题，但同时带来了模拟失真的问题，为了矫正这个约束带来的自由能贡献，最后的模拟结果会加上一个基于经验公式（通过以上6个参数计算出一个能量值）的能量矫正，以消除约束带来的不真实的结合自由能贡献。 # 以下为试图计算新冠病毒靶点7si9与辉瑞药物结合自用能的步骤，但此处着重强调分子间约束的参数设置，不会列出所有具体模拟所需的步骤，仅列出第一步作为例子。本文可以作为[github/quantaosun](https://github.com/quantaosun/Gromacs-ABF/blob/main/Gromacs_ABF.ipynb)的补充，因为在github中并没有细致讨论分子间限制的生成方法。注意区分分子内的restrain与分子间的restrain的区别和不同目的。 # Gromacs_ABFE_restrain，分子间约束的选取与参数生成。 限制势的选取有多种可能，基本原则是选取刚性较好，在模拟过程中本来就不太容易发生相对偏移的位置。小分子最好选择靠近质量中心位置的原子，蛋白最好选择氨基酸的主链。 如果你不知道如何选择，可以考虑选取能形成Pie-pie作用的临近氨基酸主链。如果没有pie-pie 可考虑形成氢键的靠近小分子中部的部分。 使用不同的限制势参数可能会影响最终的计算结果，但理论层面这种影响是可以被矫正的。 Restrain 1:  ![输入图片说明](restrain%201restrain1.png) ``` Ligand: a:4712; b:4685; c: 4707 Protein: A:2541; B:2526; C:2525 r(a-A) = 3.00 angle baA =159.2 angle aAB =111.8 dihedral cbaA =-115.5 dihedral baAB =87.9 dihedral aABC =3.9 ``` Add the following block to the end of the Gromacs topology file ``` [ intermolecular_interactions] [ bonds ] ; ai aj type bA kA bB kB 4712 2541 6 3.00 0.0 3.00 4184.0 [ angles ] ; ai aj ak type thA fcA thB fcB 4685 4712 2541 1 159.2 0.0 159.2 41.84 4712 2541 2526 1 111.8 0.0 111.8 41.84 [ dihedrals ] ; ai aj ak al type thA fcA thB fcB 4707 4685 4712 2541 2 -115.5 0.0 -115.5 41.84 4685 4712 2541 2526 2 87.9 0.0 87.9 41.84 4712 2541 2526 2525 2 3.9 0.0 3.9 41.84 ``` ------------------------------------ # 自由能模拟具体步骤例子，第一步 生成第一步文件夹内文件的代码（python）: 共晶清洗与结构分离，7SI9并不包含 BME和PO4等小分子，但这不影响我们使用下面的代码。 ``` complex = "7si9" #@param {type:"string"} ligand = "4WI" #@param {type:"string"} pdb = complex + ".pdb" LIGNAD = ligand + ".pdb" #!/home/aistudio/external-libraries/bin/pdbfixer '{pdb}' --ph=7 --replace-nonstandard --add-residues !wget wget https://files.rcsb.org/download/'{pdb}' > prot.pdb !grep -v HOH '{pdb}' > prot_clean.pdb !grep '{ligand}' prot_clean.pdb > '{ligand}'.pdb !grep HETATM '{ligand}'.pdb > LIG.pdb !grep -v '{ligand}' prot_clean.pdb > prot_clean2.pdb !grep -v PO4 prot_clean2.pdb > prot_clean3.pdb !grep -v BME prot_clean3.pdb > prot_clean4.pdb !grep ATOM prot_clean4.pdb > prot_clean5.pdb !grep -v REMARK prot_clean5.pdb > prot_clean6.pdb # here we drop all co-factors to have only protein itself #substitute all ligand name to "LIG" ``` 小分子称谓重命名从4WI到LIG，以方便后面的代码统一反复使用。 ``` Original_ligand = "4WI" #@param {type:"string"} !sed -i "s/$Original_ligand/LIG/g" LIG.pdb ``` 转换小分子格式，并对小分子加H （请提前安装Open Babel） ``` #@title Add hydrogen to LIG.pdb, and Modify small molecule bond orders #ligand_name = "LIG" #@param {type:"string"} #ligand_NAME = ligand_name + ".pdb" !obabel -ipdb LIG.pdb -omol2 -O LIG.mol2 -h # convert the format and add hydrogen atoms #!sed -i "s/'{ligand_NAME}'/'{ligand_name}'/" LIG.mol2 !cat LIG.mol2 !wget http://www.mdtutorials.com/gmx/complex/Files/sort_mol2_bonds.pl !perl sort_mol2_bonds.pl LIG.mol2 LIG_fix.mol2 !cat LIG_fix.mol2 # pay attention to the bond order section ``` 小分子的拓扑文件将需要用第三方网站Cgenff生成，然后将蛋白和小分子的结构重新整合到一起，具体 后续步骤不再粘贴在此，请参考[github/quantaosun](https://github.com/quantaosun/Gromacs-ABF/blob/main/Gromacs_ABF.ipynb)的步骤,该例子，需要在谷歌Collab中运行。 ![输入图片说明](image.png) # 自由能计算说明，github链接里面的计算流程没有采用[Gromacs2016——ABFE教程](http://www.alchemistry.org/wiki/Absolute_Binding_Free_Energy_-_Gromacs_2016)相同的分子力场，理论上这些力场的表现应当是差不多的，在此认为这不会对最终的计算结果产生较大影响。