分子動力學研究抑制劑APV與HIV-1蛋白酶的作用機制

時術華， 張少龍， 張慶剛

( 1. 山東建筑大學理學院, 濟南 250101; 2. 山東師范大學物理與電子科學學院, 濟南 250014)

分子動力學研究抑制劑APV與HIV-1蛋白酶的作用機制

時術華1， 張少龍2， 張慶剛2

( 1. 山東建筑大學理學院, 濟南 250101; 2. 山東師范大學物理與電子科學學院, 濟南 250014)

采用分子動力學模擬和結合自由能計算研究了抑制劑APV與HIV-1蛋白酶的作用機制. 研究結果表明范德瓦爾斯作用主控了APV與HIV-1蛋白酶的結合. 采用基于殘基的自由能分解方法計算了抑制劑-殘基相互作用，結果表明9個殘基Gly27、Ile32、Val47、Ile50、Ile84、Ala28′、Gly49′、Ile50′和Arg87′與APV產生了大于1.0 kcal/mol的強相互作用，而且證明CH-π，CH-O相互作用和極性作用是其結合的主要形式. 期待該結果可以為以HIV-1蛋白酶為靶標的抗艾滋病藥物設計提供理論上的指導.

分子動力學； 結合自由能； HIV-1蛋白酶； MM-PBSA

1 引 言

獲得性免疫缺陷綜合癥簡稱艾滋病, 是威脅人類生命健康的重大疾病之一. 近年來，死于艾滋病的人數呈逐年遞增趨勢, 共計2500多萬感染者死亡. 在艾滋病的生命周期中, HIV-1蛋白酶(PR)起著關鍵作用, 它能將gag和gag-pol表達的多聚蛋白分裂裝配為成熟并具有感染性的病毒顆粒[1, 2]，因此阻斷和抑制PR的活性功能是臨床上治療艾滋病的一條有效途徑.

PR由兩條相同的肽鏈組成, 每一條肽鏈均包含99個氨基酸, 結構上呈C2對稱性. PR的活性中心由保守殘基序列Asp-Thr-Gly構成, 并且其催化位點位于兩個亞基的接觸面[3, 4]. 兩個天冬氨酸Asp25和Asp25′在催化過程中起重要作用[5].

PR抑制劑已成為艾滋病聯合用藥治療方案的重要組成部分. 美國食品藥物管理局已頒布9種PR抑制劑, 其中多數藥物是以PR切割HIV前體蛋白上的活性位點為模板設計的[6]. 抑制劑Amprenaavir(APV)對PR有較強的抑制效果，其抑制常數為0.04-4.9 nM[7, 8]. 本文將從原子層次上識別APV與PR的作用機制，這對于以PR為靶標的用于治療艾滋病的高效抑制劑的研發具有重要意義.

近年來, 隨著計算理論的快速發展和計算機模擬技術的迅速提高, 分子動力學模擬和結合自由能計算成為研究抑制劑-蛋白質相互作用的重要工具[9-13]. MM-PBSA是一種基于經驗方程來快速計算抑制劑與蛋白質結合自由能的有效方法[14, 15]，該方法已經成功用于研究抑制劑-蛋白質和蛋白質-蛋白質相互作用[16-18]. 目前尚未發現APV與PR作用機制的分子動力學模擬和結合自由能計算研究，因此本文將采用MM-PBSA方法計算APV與PR的結合自由能, 研究控制APV與PR結合的主體力量，同時采用基于殘基的自由能分解方法計算抑制劑-殘基相互作用[19], 這些計算有利于從原子層次上闡明APV與PR復合體的結構親和能關系. 研究結果能為以PR為靶標的抗艾滋病藥物設計提供一定的理論啟示.

2 研究方法

取自蛋白質庫的晶體結構(3S45)用于動力學模擬的初始構象. APV與蛋白酶結合體中的結晶水分子保留在初始構象中，由Amber 12中的leap模塊添加晶體結構中缺失的氫原子[20]，蛋白質和結晶水的力場參數由Amber 12中的ff03力場產生[21]. 考慮天冬氨酸的質子化, 對殘基Asp25的氧原子O2執行了質子化. APV的力場參數源自GAFF力場[22]，使用Amber 12中的半經驗量子力學計算(AM1)給APV分配了BCC原子電荷，APV-PR結合體溶解在顯性的TIP3P的水盒子里，水盒子邊緣與復合體最近原子的距離是10.0?. 四個氯離子添加到由水和復合體組成的系統中來保證模擬系統呈現電中性.

采用Amber12中的sander模塊對整個體系執行兩步系統優化來消除一些原子間不合理的接觸：(1) 約束溶質、優化溶劑和中和離子. 約束力常數設為50 kcal/(mol·?2); (2) 無約束地優化整個系統. 在每一步優化中，首先執行2500步的最陡下降優化，接著執行2500步的共軛梯度優化；然后在500 ps內把系統從0 K加熱到300 K，隨后進行500 ps的常溫300 K、常壓1標準大氣壓條件下的動力學平衡；最后是15 ns無約束的分子動力學模擬. 模擬期間采用SHAKE方法限制所有與氫原子有關的化學鍵的伸縮[23]，模擬積分步長設為2 fs；采用PME方法用來計算長程的靜電相互作用；采用周期性邊界條件消除溶劑盒子的邊緣效應，非成鍵相互作用的截斷值為10.0 ?.

2.2 MM-PBSA方法計算結合自由能

采用單軌跡方案的MM-PBSA方法計算APV與PR的結合自由能，按照一定的時間間隔從動力學模擬平衡軌跡中抽取200個構象用于結合自由能分析，構象中刪掉水分子和氯離子，計算結合自由能的方程可表示為：

ΔG=ΔEMM+ΔGsol-TΔS

(1)

式中ΔEMM是氣相中的分子力學能，ΔGsol表示溶解自由能對APV結合的貢獻，TΔS表示熵變對結合自由能的貢獻. ΔEMM可進一步分解為：

ΔEMM=ΔEele+ΔEvdw

(2)

其中ΔEele和ΔEvdw分別表示氣相中的靜電相互作用能和范德瓦爾斯相互作用能. 溶解自由能ΔGsol可表示為：

ΔGsol=ΔGpb+ΔGsurf

(3)

其中ΔGpb和ΔGsurf分別表示極性溶解自由能和非極性溶解自由能，前者可使用Amber 12中的PBSA方法求解泊松-玻爾茲曼方程獲得，溶質和溶劑的電介常數分別設為1.0和80.0，ΔGsurf由如下經驗方程計算：

ΔGsurf=γSASA+β

(4)

式中γ和β值分別取為0.00542 kcal/(mol·?2)和0.92 kcal/mol.TΔS是由于抑制劑與蛋白酶結合前后自由度變化誘導的熵變對結合自由能的貢獻，通過使用簡振模分析和傳統的熱力學計算獲得.

2.2 分組患者HEART、MEWS評分比較 急診住院者HEART評分及MEWS評分均顯著高于留院觀察者（P＜0.05）；30 d死亡患者 HEART和 MEWS評分分值均較急診住院患者有進一步增高，且兩者分別與存活者比較，差異有統計學意義（P＜0.05）。見表1。

3 結果與討論

3.1 分子動力學的平衡

本文對APV-PR復合物系統執行了15 ns的分子動力學模擬. 為了評估動力學平衡的穩定性，使用Amber 12中的Ptraj程序計算了PR主鏈原子相對于晶體結構的均方根偏差(Root-mean-square deviation, RMSD)和蛋白結構的回旋半徑隨時間的變化關系(圖2A和B). 由圖2A觀察到，PR主鏈原子的RMSD大約在模擬開始4.5 ns后達到平衡. 結合體平衡后RMSD的平均值為1.22 ?，其漲落范圍小于0.49 ?，該結果表明PR與APV組成的系統的動力學平衡是可靠的.

圖1 分子結構：(A)抑制劑APV的結構；(B)APV與PR復合物的結構Fig.1 Structures of molecule: (A) the structure of APV; (B) the structure of the complex APV-PR

為了從整體結構上評估動力學模擬的穩定性，蛋白酶的回旋半徑隨模擬時間的漲落由圖2B給出.由圖2B觀察到，在動力學模擬開始大約4 ns后，整體結構的漲落趨向穩定，未出現異常情況. 該結果說明動力學模擬過程中，蛋白酶的結構是穩定的. 以上的RMSD和回旋半徑分析表明用于后續分析計算的動力學模擬軌跡的穩定性是可靠的.

3.2 結合自由能計算

本文采用MM-PBSA方法計算了APV與PR的結合自由能，并分析APV與蛋白酶作用的主體力量. 表1中列出了計算得到的結合自由能以及其貢獻成分.

圖2 (A) MD模擬中PR主鏈原子的RMSD；(B) MD模擬中PR結構的回旋半徑Fig.2 (A) RMSD of backbone atoms in PR during molecular dynamics simulation; (B) Gyroradius of PR during molecular dynamics

表1表明抑制劑APV與PR的結合自由能是-8.4 kcal/mol，這表明APV能夠與PR產生較強的相互作用. 由表1看到，抑制劑與蛋白酶產生的范德瓦爾斯作用能(ΔEvdw)是-66.1 kcal/mol，這表明范德瓦爾斯作用能有效地促進抑制劑與蛋白酶的結合. ΔGsurf與溶劑可及表面積相對應的非極性相互作用能, 其值為-7.0 kcal/mol, 該作用也能夠為APV與蛋白酶的結合也提供有利貢獻. 表1的結果表明抑制劑APV與PR的靜電相互作用能(ΔEele)為-40.6 kcal/mol，也有利于抑制劑與PR的結合，但完全被另一個大小為大小為57.37 kcal/mol極性溶解自由能(ΔGpb)抵消.

依據表1，熵變對結合自由能的貢獻(-TΔS)為24.35 kcal/mol，該成分消弱了抑制劑與蛋白酶的結合. 在有利的兩個相互作用成分中(ΔEvdw和ΔGsurf)，范德瓦爾斯相互作用能是非極性相互作用能的9倍之多. 從上面的分析可以得出結論，范德瓦爾斯相互作用是控制抑制劑APV與PR的結合的主要力量. 上述分析與幾個其他幾個工作組有關PR的研究結果相吻合.

表1 MM-PBSA計算所得到的能量(kcal/mol)a

aEele和Evdw分別表示靜電作用能和范德瓦爾斯作用能，Eint表示由鍵伸、鍵角彎折和二面角扭轉貢獻的內能，Egas= Eele+ Evdw+ Eint; Gsurf和Gpb分別是非極性溶劑化能和極性溶劑化能，Gsol是溶解自由能且Gsol= Gsurf+ Gpb；Gpbele= Eele+ Gpb; Gpbtot= Gsol+ Egas; TStra, TSrot和TSvib分別表示體系的平動、轉動和振動自由度熵變對自由能的貢獻；TStot= TStra+ TSrot+TSvib; ΔGbind= Gpbtot- TStot

3.3 結構-親和能關系

為了能從原子層次上更好地理解抑制劑APV與蛋白酶的作用機制，進而闡明它們的結合模式，本文采用基于殘基的自由能分解方法計算了抑制劑與蛋白酶殘基的相互作用能(圖3)，同時在圖4中顯示了PR的關鍵殘基與抑制劑APV的相對空間位置. 由圖3可以看出9個殘基Gly27、Ile32、Val47、Ile50、Ile84、Ala28′、Gly49′、Ile50′和Arg87′與抑制劑APV產生了大于1.0 kcal/mol的強相互作用. 這9個殘基為蛋白酶-抑制劑的結合提供了重要貢獻, 結構親和能關系分析如下.

圖3 抑制劑APV與PR各個分離殘基的相互作用譜Fig. 3 Inhibitor-residue interaction spectrum between the inhibitor APV and PR

從圖3可以觀察到，殘基Ala28′與APV產生了最強的相互作用，該作用能量是-2.76 kcal/mol. 其作用的主要來源是Ala28′的烷基與APV右端芳香環之間的CH-π相互作用. 在結構上，Ile50的烷基與抑制劑APV右端的芳香環相鄰近，產生較強的疏水性的CH-π相互作用；同時，Ile50的氮原子與APV的氧原子之間形成一個氫鍵(圖4), 這兩個相互作用為APV與蛋白酶的結合提供了-1.72 kcal/mol的貢獻. Ile50′與APV的相互作用為-2.03 kcal/mol，這個作用主要由Ile50′的烷基與抑制劑APV左端的芳香環形成的CH-π作用相提供. 依據圖4，殘基Ile32、Ile47和Ile84的烷基與抑制劑左端的芳香環臨近而產生CH-π作用，這三個作用分別為APV與蛋白酶的相互作用提供了-1.12、-1.12和-1.95 kcal/mol的能量貢獻. 在結構上，殘基Gly27和Gly49′的羰基氧原子與APV中間的苯環臨近，形成了有利于抑制劑結合的CH-O相互作用，這兩個相互作用分別為抑制劑結合貢獻了-1.20和-1.52 kcal/mol作用能. Arg87′與抑制劑APV產生較強的相互作用(-1.19 kcal/mol)，該相互作用主要源自Arg87′的電荷與APV的羰基產生的極性相互作用. 上述分析基本上吻合了幾個其他課題組的工作[24-26].

圖4 PR和抑制劑APV復合體中主要殘基的相對位置Fig. 4 Relative geometries of the key residues in PR-APV complex

4 結 語

本文對APV和PR的復合體系執行了15 ns的分子動力學模擬和結合自由能計算，研究了APV與PR的結合模式. 結果證明分子間的范德瓦爾斯作用主要驅動了抑制劑APV與PR的結合. 采用基于殘基的自由能分解方法計算了APV與蛋白酶各殘基的相互作用，結果表明9個殘基Gly27、Ile32、Val47、Ile50、Ile84、Ala28′、Gly49′、Ile50′和Arg87′與APV產生較強相互作用, 結構親和能關系分析證明CH-π，CH-O相互作用和極性作用驅動了APV與PR的結合. 該研究能為以PR為靶標的治療艾滋病的藥物設計提供重要理論啟示.

[1] Judd D A, Nettles J H, Nevins N,etal. Polyoxometalate HIV-1 protease inhibitors. A new mode of protease inhibition [J].J.Am.Chem.Soc., 2001, 123: 886.

[2] Tie Y, Wang Y F, Boross P I,etal. Critical differences in HIV‐1 and HIV‐2 protease specificity for clinical inhibitors [J].ProteinScience, 2012, 21: 339.

[3] Blum A, B?ttcher J, Heine A,etal. Structure-guided design of C 2-symmetric HIV-1 protease inhibitors based on a pyrrolidine scaffold? [J].J.Med.Chem., 2008, 51: 2078.

[4] Brik A, Wong C H. HIV-1 protease: mechanism and drug discovery [J].Org.Biomol.Chem., 2003, 1: 5.

[5] Perryman A L, Lin J H, McCammon J A. HIV‐1 protease molecular dynamics of a wild‐type and of the V82F/I84V mutant: Possible contributions to drug resistance and a potential new target site for drugs [J].ProteinScience, 2009, 13: 1108.

[6] Soriano V, Gomes P, Heneine W,etal. Human immunodeficiency virus type 2 (HIV‐2) in Portugal: Clinical spectrum, circulating subtypes, virus isolation, and plasma viral load [J].J.Med.Virol., 2000, 61: 111.

[7] Shen C H, Wang Y F, Kovalevsky A Y,etal. Amprenavir complexes with HIV‐1 protease and its drug‐resistant mutants altering hydrophobic clusters [J].FEBS.J., 2010, 277: 3699.

[8] Tie Y, Wang Y F, Boross P I,etal. Critical differences in HIV‐1 and HIV‐2 protease specificity for clinical inhibitors [J].ProteinScience, 2012, 21: 339.

[9] Chen J, Wang J, Zhu W,etal. A computational analysis of binding modes and conformation changes of MDM2 induced by p53 and inhibitor bindings [J].J.Comput.AidedMol.Des., 2013, 27: 965.

[10] 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- AEuropeanJournal, 2008, 14: 8704.

[11] Hu G D, Zhang S L, Liu X G,etal. Molecular dynamics and free energy calculation study progesterone binding to human progesterone receptor-ligand binding domain [J].J.At.Mol.Phys., 2010, 27: 333(in Chinese)[扈國棟, 張少龍, 劉新國, 等. 分子動力學模擬和自由能計算研究孕酮和孕酮蛋白受體的結合模式 [J]. 原子與分子物理學報, 2010, 27: 333]

[12] Cheng W Y, Liang Z Q, Zhang Q G,etal. Insight into p53-MDM2 interaction based on molecular dynamics simulation and molecular mechanics [J].J.At.Mol.Phys., 2012, 29: 393(in Chinese)[程偉淵, 梁志強, 張慶剛, 等. p53-MDM2 相互作用的分子力學和動力學研究 [J]. 原子分子物理學報, 2012, 29: 393]

[13] Chen J, Wang J, Xu B,etal. Insight into mechanism of small molecule inhibitors of the MDM2-p53 interaction: molecular dynamics simulation and free energy analysis [J].J.Mol.Graph.Model., 2011, 30: 46.

[14] Wang J, Morin P, Wang W,etal. Use of MM-PBSA in reproducing the binding free energies to HIV-1 RT of TIBO derivatives and predicting the binding mode to HIV-1 RT of efavirenz by docking and MM-PBSA [J].J.Am.Chem.Soc., 2001, 123: 5221.

[15] Wang W, Kollman P A. Free energy calculations on dimer stability of the HIV protease using molecular dynamics and a continuum solvent model1 [J].J.Mol.Biol., 2000, 303: 567.

[16] Yi C, Chen J, Zhu T,etal. Protonation of Asp25 of HIV-1 protease stabilizes its binding to the inhibitor GRL02031 [J].J.At.Mol.Phys., 2011, 28: 11(in Chinese)[伊長虹, 陳建中, 朱通, 等. Asp25質子化有利于HIV-1蛋白酶與抑制劑GRL02031結合 [J]. 原子分子物理學報, 2011, 28: 11]

[17] Chen J, Zhang D, Zhang Y,etal. Computational studies of difference in binding modes of peptide and non-peptide inhibitors to MDM2/MDMX based on molecular dynamics simulations [J].Int.J.Mol.Sci., 2012, 13: 2176.

[18] Shi S H, Chen J Z, Hu G D,etal. Molecular insight into the interaction mechanisms of inhibitors BEC and BEG with HIV-1 protease by using MM-PBSA method and molecular dynamics simulation [J].J.Mol.Struct.:THEOCHEM, 2009, 913: 22.

[19] Gohlke H, Kiel C, Case D A. Insights into protein-protein binding by binding free energy calculation and free energy decomposition for the Ras-Raf and Ras-RalGDS complexes [J].J.Mol.Biol., 2003, 330: 891.

[20] Case D A, Darden T A, Cheatham III T E,etal. AMBER 12[J].SanFrancisco:UniversityofCalifornia, 2012.

[21] Duan Y, Wu C, Chowdhury S,etal. A point-charge force field for molecular mechanics simulations of proteins based on condensed-phase quantum mechanical calculations [J].J.Comp.Chem., 2003, 24: 1999.

[22] Wang J, Wolf R M, Caldwell J W,etal. Development and testing of a general amber force field [J].J.Comput.Chem., 2004, 25: 1157.

[23] Coleman T G, Mesick H C, Darby R L. Numerical integration [J].Ann.Biomed.Eng., 1977, 5: 322.

[24] B?ttcher J, Blum A, Heine A,etal. Structural and kinetic analysis of pyrrolidine-based inhibitors of the drug-resistant Ile84Val mutant of HIV-1 protease [J].J.Mol.Biol., 2008, 383: 347.

[25] Chen J, Zhang S, Liu X,etal. Insights into drug resistance of mutations D30N and I50V to HIV-1 protease inhibitor TMC-114: Free energy calculation and molecular dynamic simulation [J].J.Mol.Model., 2010, 16: 459.

[26] Hou T, Yu R. Molecular dynamics and free energy studies on the wild-type and double mutant HIV-1 protease complexed with amprenavir and two amprenavir-related inhibitors: mechanism for binding and drug resistance [J].J.Med.Chem., 2007, 50: 1177.

Molecular dynamics insights into binding mode of inhibitor APV to HIV-1 protease

SHI Shu-Hua1, ZHANG Shao-Long2, ZHANG Qing-Gang2

(1. School of Science, Shandong Jianzhu University, Jinan 250101, China；2.College of Physics and Electronics, Shandong Normal University, Jinan 250014, China)

Molecular dynamics simulation and binding free energy calculation was performed to study the binding mode of inhibitor APV to HIV-1 protease. The results show that van der walls energy is the main force driving the binding of APV to HIV-1 protease. Residue-based free energy decomposition was adopted to calculate the inhibitor-residue interaction and the results suggest that nine residues Gly27, Ile32, Val47, Ile50, Ile84, Ala28′, Gly49′, Ile50′ and Arg87′ in HIV-1 protease produce strong interactions with APV, at the same time, these results also show that the CH-π， CH-O and polar interaction are the main ones between APV and HIV-1 protease. We expect that this study can provide significant contributions to the designs of anti-AIDS targeting HIV-1 protease.

Molecular dynamics; Binding free energy; HIV-1 protease; MM-PBSA

2014-02-11

國家自然科學基金(11274206); 山東省自然科學基金(ZR2011HM048,ZR2013AM014); 山東建筑大學博士啟動基金(XNBS1268)

時術華(1967—)，女, 山東文登人, 博士, 副教授, 主要研究領域為生物大分子的理論計算.E-mail: sdsfhf@sdjzu.edu.cn

103969/j.issn.1000-0364.2015.10.028

Q641.12

A

1000-0364(2015)05-0885-06