HIV-1蛋白酶与抑制剂BEG 相互作用的分子对接计算研究

伊长虹，梁志强，王伟，张少龙，张庆刚，陈建中

（1.山东交通学院理学院，济南250023；2.山东师范大学物理与电子科学学院，济南250014）

1 前 言

人类免疫缺陷病毒蛋白酶 （HIV-1protease）是HIV-1病毒生命周期中的关键酶之一.HIV-1蛋白酶的作用就是负责分解病毒颗粒的蛋白前体，促使病毒的成熟，使之变成具有感染性的病毒颗粒，HIV-1病毒与抑制剂BEG 的三维空间结构如图1所示，它是一个对称的同质二聚体，每一条单体分别由99个氨基酸残基组成.图1中显示的残基Asp25和Asp25′是HIV-1蛋白酶的催化中心，参与底物肽链的分解.氨基酸序列的残基段45-55和45′-55′处于HIV-1蛋白酶的柔性中心.在大多数的HIV-1蛋白酶-抑制剂的复合体中，结晶水分子Wat301位于蛋白酶的残基Ile50／Ile50′和抑制剂以及底物之间，通过4 个氢键形成了抑制剂与HIV-1蛋白酶相互作用的桥梁［1，2］.

图1 HIV-1蛋白酶与抑制剂BEG 复合体的结构 （抑制剂BEG 使用线状显示，氨基酸残基Asp25 和Ile50使用球棒方式表示，水分子Wat301使用球状显示）Fig.1 Structures of HIV-1protease complexed with inhibitor BEG，BEG is shown in a line mode，Asp25and Ile50are shown in a ball-and-stick mode，Wat301is shown in a ball mode

国内外许多课题组都在努力研究抑制剂与HIV-1蛋白酶的相互作用机制，如果能够非常清晰地阐明抑制剂与HIV-1 蛋白酶的相互作用机制，必定能够在治疗艾滋病药物的合理化设计方面起到非常重要的作用，也必定能加速药物研发的过程，缩短药物的研发周期.目前，许多的实验工作已经揭示了一些抑制剂与HIV-1蛋白酶的相互作用机制［3-4］.尽管实验上取得了很大成功，但由于两个重要事实的存在，仍限制了药物的研发，（1）HIV-1蛋白酶具有很强的变异能力，限制了药效的发挥；（2）缺乏HIV-1蛋白酶与抑制剂相互作用的原子层次的了解.而分子计算生物学可以弥补这个缺陷，与实验研究形成互补，可以提高药物的研发的效率.

结合自由能计算是计算机模拟研究蛋白质与抑制剂相互作用强度的重要工具.这种方法不仅能在原子和分子的层次上给出蛋白质和抑制剂的相互作用细节，而且还能便于阐明蛋白质-抑制剂复合体的结构亲和能关系.MM／PBSA方法已经成功地用于研究HIV-1 蛋白酶和抑制剂的相互作用机制［5-10］.因此，我们将应用MM／PBSA 方法研究抑制剂BEG和HIV-1蛋白酶的相互作用机制.图1给出了HIV-1蛋白酶和抑制剂BEG 复合体的三维结构，图2给出了抑制剂的BEG的分子结构.

分子对接是研究蛋白质与抑制剂相互作用的另一个计算模拟的重要工具.分子对接不仅能较准确地探测出抑制剂与蛋白质作用的靶点，而且还能描述抑制剂与蛋白质间力的作用性质、氢键和范德华作用的匹配性，给出抑制剂与蛋白质的结合自由能，研究抑制剂与蛋白质的结合模式.国内外的几个课题组使用柔性的分子对接研究了抑制剂和HIV-1蛋白酶识别的分子基础［11-14］.

我们将使用MM／PBSA 方法和分子对接计算HIV-1蛋白酶与抑制剂BEG 的结合自由能，研究抑制剂与蛋白酶的结合模式，探索驱动抑制剂与蛋白酶结合的主要力量，为治疗艾滋病药物的合理化设计提供理论上的指导.

图2 HIV-1蛋白酶抑制剂BEG 的分子结构Fig.2 Molecular structure of HIV-1protease and inhibitor BEG

2 研究方法

2.1 分子动力学模拟

动力学模拟的初始构象取自蛋白质库的晶体结构（1D4I）.复合物中的结晶水分子保留在初始构象中，晶体结构中缺失的氢原子由Amber 12中的leap 模块添加［15］，复合物的力场参数由Amber 12中的ff03力场产生［16］.抑制剂BEG 的力场参数源自GAFF 力场［17］，BEG 的原子电荷由Amber 12中的半经验的量子力学 （AM1）的计算分配获得，BEG-PR 复合体溶解在显性的TIP3P的水盒子里，水盒子边缘与复合体最近原子的距离是10.0Å.为保证模拟系统呈电中性，四个氯离子添加到由水和复合体组成的系统.

为了消除原子间一些不合理的接触，采用Amber 12中的sander模块对复合体执行两步的系统优化：（1）约束溶质.优化溶剂和中和离子，约束力常数设为100kcal／ （mol·Å2）；（2）无约束地优化整个系统.每一步优化均先执行1000步的最陡下降优化，接着执行3000步的共轭梯度优化；然后在100ps内把系统从0K 加热到300K，随后进行100ps的常温300K，常压1标准大气压条件下的动力学平衡；最后是10ns的无约束分子动力学模拟.模拟期间采用SHAKE 方法限制所有含氢原子化学键的伸缩，模拟积分步长为2fs，PME 方法用来计算长程的静电相互作用，应用周期性边界条件以消除溶剂盒子的边缘效应，非成键相互作用的截断值为10.0Å.

2.2 MM-PBSA 方法计算结合自由能

采用单轨迹方案的MM-PBSA 方法计算BEG与HIV-1 蛋白酶的结合自由能，结合自由能计算使用的构象以一定的间隔取自动力学模拟轨迹，构象中删掉水分子和氯离子，结合自由能由如下方程计算：

式中ΔEMM是气相中的分子力学能，ΔGsol表示溶解自由能对分子结合的贡献，TΔS表示熵变对结合自由能的贡献.，ΔEMM由两部分组成：

其中ΔEele和ΔEvdw分别表示气相中的静电相互作用能和范德华相互作用能，溶解自由能ΔGsol可分解为如下两项：

其中ΔGpb和ΔGsurf分别表示极性溶解自由能和非极性溶解自由能，前者可使用Amber 12中的pbsa方法求解泊松-玻尔兹曼方程获得，溶质和溶剂的电介常数分别设为1.0和80.0，后者由如下经验方程计算：

式中γ和β值分别取为0.00542kcal／ （mol·Å2）和0.92kcal／mol，TΔS是由于自由度变化导致的熵变对结合自由能的贡献，该项使用简振模和传统的热力学计算.

2.3 分子对接方法

分子对接是依据配体与受体的 “锁钥原理”，模拟小分子配体与受体 （诸如蛋白质等生物大分子）的相互作用.配体与受体相互作用主要包括静电相互作用，氢键作用、疏水作用和范德华作用等.通过对空间匹配和能量匹配的计算，预测配体与受体间的结合模式和亲和力，从而获得最佳的构象结构.空间匹配的计算通常采用格点计算和片段生长的方法，而能量计算采用局部搜索、模拟退火和遗传算法等方法［18，19］.本文主要采用Autodock对接程序计算HIV-1 蛋白酶与抑制剂的结合自由能，探测它们的结合模式，为药物的合理化设计提供一种指导.

3 结果和分析

3.1 分子动力学的平衡

本文对BEG-PR 复合物系统执行了10ns的分子动力学模拟，为了评价系统动力学平衡的稳定性，我们用Amber 12程序中的Ptraj工具计算了HIV-1 蛋白酶主链原子相对于晶体结构的RMSD 随时间的变化关系 （图3）.由图3 观察到，蛋白酶主链原子的RMSD 在前4ns之内波动比较大，说明系统在消除原子间不合理的接触，而系统在6ns后达到了平衡.系统平衡后RMSD的平均值为1.24Å，涨落范围低于0.6Å，这表明HIV-1蛋白酶与BEG 组成的系统的动力学平衡是可靠的.这表明动力学模拟过程中，蛋白酶的结构是稳定的.上述分析表明用于后加工分析的动力学模拟轨迹的稳定性是可靠的.

图3 MD 模拟中HIV-1蛋白酶主链原子的均方根偏差Fig.3 Root-Mean-Square-Deviation of the backbone atoms from HIV-1protease

表1 MM-PBSA 计算所得到的能量 （kcal·mol-1 ）Table 1 Binding free energies calculated using MM-PBSA method

3.2 采用MM／PBSA 方法计算结合自由能

基于动力学模拟比较浪费时间，我们只使用MM-PBSA 方法和优化后的单个构象计算抑制剂BEG 与HIV-1 蛋白酶的结合自由能，以便考察什么性质的力驱动两者的相互作用，结合自由能各成分的贡献均列在表1中.

据表1，抑制剂BEG 与HIV-1蛋白酶的结合自由能是-22.25kcal／mol，表明抑制剂能与HIV-1蛋白酶有比较强的结合能力.范德华作用能（VDW）是-61.41kcal／mol，这个作用成分促进了抑制剂与蛋白的结合.与溶剂可及表面积相对应的非极性成分相互作用能 （PBSUR）是 -6.92kcal／mol，这个作用成分也为抑制剂的结合提供了较小的有利贡献.虽然抑制剂与HIV-1蛋白酶的静电相互作用 （ELE）是-56.64kcal／mol，它也有利于抑制剂的结合，但是这种有力的作用成分被更强的极性溶解自由能所中和.从表一中可以看到熵变对对结合自由能 （TSTOT）的贡献消弱了抑制剂与蛋白酶的结合.在有利的两个相互作用成分中，范德华相互作用能是非极性相互作用的10倍.因此得出结论范德华相互作用能驱动了抑制剂BEG 与HIV-1蛋白酶的结合.

ELE和VDW 分别表示静电作用能和范德华作用能，INT表示由键伸、键角弯折和二面角扭转贡献的内能，Gas=ELE+VDW+INT；PBSUR 和PBCAL分别是非极性溶剂化能和极性溶剂化能，PBSOL溶解自由能且PBSOL=PBSUR+PBCAL；PBELE=ELE+PBCAL；PBTOT=PBSOL+GAS；TSTRA，TSTOT 和TSVIB分别表示体系的平动、转动和振动对熵的贡献；TSTOT=TSTRA+TSTOT+TSVIB；ΔGbind=PBTOT-TSTOT.

3.3 分子对接计算结合自由能

分子对接是一种基于经典力场和打分函数评估抑制剂与蛋白质的结合，能较好地探测药物的靶标，合理地预测和分析抑制剂与蛋白质的结合模式，能为药物的合理化设计提供理论指导.这一部分我们将利用分子对接计算BEG 与HIV-1蛋白酶的结合自由能，以便能评估抑制剂与HIV-1蛋白酶的结合效果.

分子对接的计算结果如表2显示，其中计算得到的抑制常数KI是5.82nM，分子间的相互作用能，内能和扭转能分别是-14.93，-0.88 和3.84 kcal／mol，未成键的排斥能为-0.74kcal／mol.依据表2可以看到BEG 与HIV-1蛋白酶的结合自由能为-11.23kcal／mol，这与实验值的结果 （-12.29 kcal／mol）［20］吻合的非常好.这个计算结果表明分子对接软件的所使用的力场比Amber 12力场更适合于BEG-HIV-1蛋白酶系统.

表2 分子对接计算的结合自由能 （kcal·mol-1 ）Table 2 Binding free energies calculated using the molecular docking method

3.4 抑制剂BEG 与HIV-1蛋白酶结合模式的分析

为了能更好地阐明抑制剂BEG 与HIV-1 蛋白酶的结合模式，我们采用MM／GBSA 方法计算抑制剂-残基相互作用谱，如图4所示.我们只在图中显示了作用能大于1.2kcal／mol的残基.

图4 抑制剂BEG 与HIV-1蛋白酶残基相互作用谱Fig.4 Interaction spectrum for the inhibitor／HIV-1 protease binding complex calculated using GBSA method

图5 与抑制剂BEG 产生较强相互作用的相关残基Fig.5 Key residues in HIV-1protease-BEG complex

据图4，A 链中的残基Leu23、Val82 和Ile84与BEG 的相互作用能分别为-1.23、-1.56和-2.08kcal／mol，从图5的分子结构看，这些相互作用能极好的吻合了三个残基的烷基与BEG 苯环间的CH-π相互作用.残基Gly49与BEG 的相互作用能为-2.08kcal／mol，从图5 中的结构看，这个能量主要来源于Gly49的羰基与苯环和五元环间的CH-O 相互作用.依据图5 可以看到，Ile50的烷基与BEG 的苯环相互靠近，这形成了比较强烈的CH-π 的相互作用，该能量值为-3.28kcal／mol.图4表明B链中的Asp25′与BEG有一个1.02kcal／mol的正的相互作用能，该能量主要因为Asp25′的两个带负电荷的氧原子与BEG骨架上的两个氧原子相距较近，从而产生了强烈的排斥作用.Gly27′与BEG 的相互作用能为-1.96kcal／mol，该能量主要来源于Gly27′与抑制剂BEG 间的CH-O.从图5 的结构中可以观察到，B链中残基Ala28′的烷基与BEG 的苯环相邻近，因此能形成CH-π相互作用，从图4获知该相互作用能量值大约为-1.49kcal／mol.依据HIV-1蛋白酶的C2 对称性，B 链中的三个残基Gly49′、Ile50′和Ile84′与抑制剂BEG的相互作用分析类似于A 链中相对应的三个残基.总之，从上面的分析可以看到CH-π和CH-O 相互作用驱动了抑制剂BEG与HIV-1蛋白酶的结合.

4 结 论

本文用分子对接方法和MM／PBSA 方法计算抑制剂BEG 与HIV-1蛋白酶的结合自由能，结果表明范德华相互作用主宰了BEG 与HIV-1蛋白酶的相互作用.用基于残基的相互作用能分解方法计算了抑制剂-残基的相互作用，结果表明CH-π和CH-O 相互作用驱动了抑制剂BEG 与HIV-1蛋白酶的结合.我们期望这个研究能为艾滋病治疗药物的合理化设计提供理论上的指导.

［1］ Böttcher J，Blum A，et al.Structural and kinetic analysis of pyrrolidine-based inhibitors of the drug-resistant Ile84Val mutant of HIV-1protease［J］.J.Mol.Biol.，2008，383：347.

［2］ Wang W，Kollman P A.Computational study of protein specificity：the molecular basis of HIV-1protease drug resistance［J］.Proc.Natl.Acad.Sci.，2001，98：14937.

［3］ Erickson J，Neidhart D J，VanDrie J，et al.Design，activity，and 2.8 Acrystal structure of a C2symmetric inhibitor complexed to HIV-1 protease［J］.Science，1990，249：527.

［4］ Whiting M，Muldoon J，Lin Y C，et al.Inhibitors of HIV-1protease by using in situ click chemistry［J］.Angew.Chem.Int.Ed.，2006，45：1435.

［5］ Chen J，Yang M，Hu G，et al.Insights into the functional role of protonation states in the HIV-1 protease-BEA369complex：molecular dynamics simulations and free energy calculations［J］.J.Mol.Model.，2009，15：1245.

［6］ Chen J，Zhang S，Liu X，Zhang Q.Insights into drug resistance of mutations D30N and I50V to HIV-1protease inhibitor TMC-114：Free energy calculation and molecular dynamic simulation［J］.J.Mol.Model.，2010，16：459.

［7］ Chen J Z，Yang M Y，Yi C H，et al.Molecular dynamics simulation and free energy calculations of symmetric fluoro-substituted diol-based HIV-1protease inhibitors ［J］.J.Mol.Struct.：THEOCHEM，2009，899：1.

［8］ Yi C H，Zhang Q G.Study of binding free energies between HIV-1protease and its inhibitors［J］.Acta Chim.Sinica，2010，68：2029（in Chinese）

［9］ Yi C H，Chen J Z，Zhu T，Zhang Q G，Protonation of Asp25of HIV-1protease stabilizes its binding to the inhibitor GRL02031［J］.J.of A.and Mol.Phy.（原子与分子物理学报），2011，28：11（in Chinese）

［10］ Shi S H，Chen J Z，Hu G D，et al.Molecular insight into the interaction mechanisms of inhibitors BEC and BEG with HIV-1protease by using MMPBSA method and molecular dynamics simulation［J］.J.Mol.Struct.：THEOCHEM，2009，913：22.

［11］ Gehlhaar D K，Verkhivker G M，Rejto P A，et al.Molecular recognition of the inhibitor AG-1343by HIV-1 protease：conformationally flexible docking by evolutionary programming ［J］.Chem.Biol.，1995，2：317.

［12］ Schaffer L，Verkhivker G M.Predicting structural effects in HIV-1 protease mutant complexes with flexible ligand docking and protein side-chain optimization［J］.Proteins.，1998，33：295.

［13］ Vieth M，Cummins D J.DoMCoSAR：a novel approach for establishing the docking mode that is consistent with the structure-activity relationship.Application to HIV-1protease inhibitors and VEGF receptor tyrosine kinase inhibitors ［J］.J.Med.Chem.，2000，43：3020.

［14］ Wu E L，Han K L，Zhang J Z H.Selectivity of Neutral／Weakly Basic P1 Group Inhibitors of Thrombin and Trypsin by a Molecular Dynamics Study［J］.Chemistry-A European J，2008，14：8704.

［15］ Case D A，Darden T A，Cheatham III T E，et al.AMBER 12［J］.University of California，San Francisco，2012.

［16］ Duan Y，Wu C，Chowdhury S，et al.A pointcharge force field for molecular mechanics simulations of proteins based on condensed-phase quantum mechanical calculations ［J］.J.Comput.Chem.2003，24：1999.

［17］ Wang J，Wolf R M，Caldwell J W，et al.Development and testing of a general amber force field［J］.J.Comput.Chem.，2004，25：1157.

［18］ Li X R，Wei D Q，Feng Y，et al.Simulation of design of Protein structures［M］.Beijing：Chemical Industry Press，2011：11-156.

［19］ Morris G M，Goodsell D S，Halliday R S，et al.Automated docking using a Lamarckian genetic algorithm and an empirical binding free energy function［J］.J.Comput.Chem.，1998，19：1639.

［20］ Jimmy L，David P，Seved L，et al.Symmetric fluoro-substituted diol-based HIV protease inhibitors Ortho-fluorinated and meta-fluorinated P1／P1′-benzyloxy side groups significantly improve the antiviral activity and preserve binding efficacy［J］.Eur.J.Biochem.，2004，271：4594.