HIV-1蛋白酶與抑制劑BEG 相互作用的分子對接計(jì)算研究

伊長虹，梁志強(qiáng)，王偉，張少龍，張慶剛，陳建中

（1.山東交通學(xué)院理學(xué)院，濟(jì)南250023；2.山東師范大學(xué)物理與電子科學(xué)學(xué)院，濟(jì)南250014）

1 前 言

人類免疫缺陷病毒蛋白酶 （HIV-1protease）是HIV-1病毒生命周期中的關(guān)鍵酶之一.HIV-1蛋白酶的作用就是負(fù)責(zé)分解病毒顆粒的蛋白前體，促使病毒的成熟，使之變成具有感染性的病毒顆粒，HIV-1病毒與抑制劑BEG 的三維空間結(jié)構(gòu)如圖1所示，它是一個(gè)對稱的同質(zhì)二聚體，每一條單體分別由99個(gè)氨基酸殘基組成.圖1中顯示的殘基Asp25和Asp25′是HIV-1蛋白酶的催化中心，參與底物肽鏈的分解.氨基酸序列的殘基段45-55和45′-55′處于HIV-1蛋白酶的柔性中心.在大多數(shù)的HIV-1蛋白酶-抑制劑的復(fù)合體中，結(jié)晶水分子Wat301位于蛋白酶的殘基Ile50／Ile50′和抑制劑以及底物之間，通過4 個(gè)氫鍵形成了抑制劑與HIV-1蛋白酶相互作用的橋梁［1，2］.

圖1 HIV-1蛋白酶與抑制劑BEG 復(fù)合體的結(jié)構(gòu) （抑制劑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

國內(nèi)外許多課題組都在努力研究抑制劑與HIV-1蛋白酶的相互作用機(jī)制，如果能夠非常清晰地闡明抑制劑與HIV-1 蛋白酶的相互作用機(jī)制，必定能夠在治療艾滋病藥物的合理化設(shè)計(jì)方面起到非常重要的作用，也必定能加速藥物研發(fā)的過程，縮短藥物的研發(fā)周期.目前，許多的實(shí)驗(yàn)工作已經(jīng)揭示了一些抑制劑與HIV-1蛋白酶的相互作用機(jī)制［3-4］.盡管實(shí)驗(yàn)上取得了很大成功，但由于兩個(gè)重要事實(shí)的存在，仍限制了藥物的研發(fā)，（1）HIV-1蛋白酶具有很強(qiáng)的變異能力，限制了藥效的發(fā)揮；（2）缺乏HIV-1蛋白酶與抑制劑相互作用的原子層次的了解.而分子計(jì)算生物學(xué)可以彌補(bǔ)這個(gè)缺陷，與實(shí)驗(yàn)研究形成互補(bǔ)，可以提高藥物的研發(fā)的效率.

結(jié)合自由能計(jì)算是計(jì)算機(jī)模擬研究蛋白質(zhì)與抑制劑相互作用強(qiáng)度的重要工具.這種方法不僅能在原子和分子的層次上給出蛋白質(zhì)和抑制劑的相互作用細(xì)節(jié)，而且還能便于闡明蛋白質(zhì)-抑制劑復(fù)合體的結(jié)構(gòu)親和能關(guān)系.MM／PBSA方法已經(jīng)成功地用于研究HIV-1 蛋白酶和抑制劑的相互作用機(jī)制［5-10］.因此，我們將應(yīng)用MM／PBSA 方法研究抑制劑BEG和HIV-1蛋白酶的相互作用機(jī)制.圖1給出了HIV-1蛋白酶和抑制劑BEG 復(fù)合體的三維結(jié)構(gòu)，圖2給出了抑制劑的BEG的分子結(jié)構(gòu).

分子對接是研究蛋白質(zhì)與抑制劑相互作用的另一個(gè)計(jì)算模擬的重要工具.分子對接不僅能較準(zhǔn)確地探測出抑制劑與蛋白質(zhì)作用的靶點(diǎn)，而且還能描述抑制劑與蛋白質(zhì)間力的作用性質(zhì)、氫鍵和范德華作用的匹配性，給出抑制劑與蛋白質(zhì)的結(jié)合自由能，研究抑制劑與蛋白質(zhì)的結(jié)合模式.國內(nèi)外的幾個(gè)課題組使用柔性的分子對接研究了抑制劑和HIV-1蛋白酶識別的分子基礎(chǔ)［11-14］.

我們將使用MM／PBSA 方法和分子對接計(jì)算HIV-1蛋白酶與抑制劑BEG 的結(jié)合自由能，研究抑制劑與蛋白酶的結(jié)合模式，探索驅(qū)動抑制劑與蛋白酶結(jié)合的主要力量，為治療艾滋病藥物的合理化設(shè)計(jì)提供理論上的指導(dǎo).

圖2 HIV-1蛋白酶抑制劑BEG 的分子結(jié)構(gòu)Fig.2 Molecular structure of HIV-1protease and inhibitor BEG

2 研究方法

2.1 分子動力學(xué)模擬

動力學(xué)模擬的初始構(gòu)象取自蛋白質(zhì)庫的晶體結(jié)構(gòu)（1D4I）.復(fù)合物中的結(jié)晶水分子保留在初始構(gòu)象中，晶體結(jié)構(gòu)中缺失的氫原子由Amber 12中的leap 模塊添加［15］，復(fù)合物的力場參數(shù)由Amber 12中的ff03力場產(chǎn)生［16］.抑制劑BEG 的力場參數(shù)源自GAFF 力場［17］，BEG 的原子電荷由Amber 12中的半經(jīng)驗(yàn)的量子力學(xué) （AM1）的計(jì)算分配獲得，BEG-PR 復(fù)合體溶解在顯性的TIP3P的水盒子里，水盒子邊緣與復(fù)合體最近原子的距離是10.0?.為保證模擬系統(tǒng)呈電中性，四個(gè)氯離子添加到由水和復(fù)合體組成的系統(tǒng).

為了消除原子間一些不合理的接觸，采用Amber 12中的sander模塊對復(fù)合體執(zhí)行兩步的系統(tǒng)優(yōu)化：（1）約束溶質(zhì).優(yōu)化溶劑和中和離子，約束力常數(shù)設(shè)為100kcal／ （mol·?2）；（2）無約束地優(yōu)化整個(gè)系統(tǒng).每一步優(yōu)化均先執(zhí)行1000步的最陡下降優(yōu)化，接著執(zhí)行3000步的共軛梯度優(yōu)化；然后在100ps內(nèi)把系統(tǒng)從0K 加熱到300K，隨后進(jìn)行100ps的常溫300K，常壓1標(biāo)準(zhǔn)大氣壓條件下的動力學(xué)平衡；最后是10ns的無約束分子動力學(xué)模擬.模擬期間采用SHAKE 方法限制所有含氫原子化學(xué)鍵的伸縮，模擬積分步長為2fs，PME 方法用來計(jì)算長程的靜電相互作用，應(yīng)用周期性邊界條件以消除溶劑盒子的邊緣效應(yīng)，非成鍵相互作用的截?cái)嘀禐?0.0?.

2.2 MM-PBSA 方法計(jì)算結(jié)合自由能

采用單軌跡方案的MM-PBSA 方法計(jì)算BEG與HIV-1 蛋白酶的結(jié)合自由能，結(jié)合自由能計(jì)算使用的構(gòu)象以一定的間隔取自動力學(xué)模擬軌跡，構(gòu)象中刪掉水分子和氯離子，結(jié)合自由能由如下方程計(jì)算：

式中ΔEMM是氣相中的分子力學(xué)能，ΔGsol表示溶解自由能對分子結(jié)合的貢獻(xiàn)，TΔS表示熵變對結(jié)合自由能的貢獻(xiàn).，ΔEMM由兩部分組成：

其中ΔEele和ΔEvdw分別表示氣相中的靜電相互作用能和范德華相互作用能，溶解自由能ΔGsol可分解為如下兩項(xiàng)：

其中ΔGpb和ΔGsurf分別表示極性溶解自由能和非極性溶解自由能，前者可使用Amber 12中的pbsa方法求解泊松-玻爾茲曼方程獲得，溶質(zhì)和溶劑的電介常數(shù)分別設(shè)為1.0和80.0，后者由如下經(jīng)驗(yàn)方程計(jì)算：

式中γ和β值分別取為0.00542kcal／ （mol·?2）和0.92kcal／mol，TΔS是由于自由度變化導(dǎo)致的熵變對結(jié)合自由能的貢獻(xiàn)，該項(xiàng)使用簡振模和傳統(tǒng)的熱力學(xué)計(jì)算.

2.3 分子對接方法

分子對接是依據(jù)配體與受體的 “鎖鑰原理”，模擬小分子配體與受體 （諸如蛋白質(zhì)等生物大分子）的相互作用.配體與受體相互作用主要包括靜電相互作用，氫鍵作用、疏水作用和范德華作用等.通過對空間匹配和能量匹配的計(jì)算，預(yù)測配體與受體間的結(jié)合模式和親和力，從而獲得最佳的構(gòu)象結(jié)構(gòu).空間匹配的計(jì)算通常采用格點(diǎn)計(jì)算和片段生長的方法，而能量計(jì)算采用局部搜索、模擬退火和遺傳算法等方法［18，19］.本文主要采用Autodock對接程序計(jì)算HIV-1 蛋白酶與抑制劑的結(jié)合自由能，探測它們的結(jié)合模式，為藥物的合理化設(shè)計(jì)提供一種指導(dǎo).

3 結(jié)果和分析

3.1 分子動力學(xué)的平衡

本文對BEG-PR 復(fù)合物系統(tǒng)執(zhí)行了10ns的分子動力學(xué)模擬，為了評價(jià)系統(tǒng)動力學(xué)平衡的穩(wěn)定性，我們用Amber 12程序中的Ptraj工具計(jì)算了HIV-1 蛋白酶主鏈原子相對于晶體結(jié)構(gòu)的RMSD 隨時(shí)間的變化關(guān)系 （圖3）.由圖3 觀察到，蛋白酶主鏈原子的RMSD 在前4ns之內(nèi)波動比較大，說明系統(tǒng)在消除原子間不合理的接觸，而系統(tǒng)在6ns后達(dá)到了平衡.系統(tǒng)平衡后RMSD的平均值為1.24?，漲落范圍低于0.6?，這表明HIV-1蛋白酶與BEG 組成的系統(tǒng)的動力學(xué)平衡是可靠的.這表明動力學(xué)模擬過程中，蛋白酶的結(jié)構(gòu)是穩(wěn)定的.上述分析表明用于后加工分析的動力學(xué)模擬軌跡的穩(wěn)定性是可靠的.

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

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

3.2 采用MM／PBSA 方法計(jì)算結(jié)合自由能

基于動力學(xué)模擬比較浪費(fèi)時(shí)間，我們只使用MM-PBSA 方法和優(yōu)化后的單個(gè)構(gòu)象計(jì)算抑制劑BEG 與HIV-1 蛋白酶的結(jié)合自由能，以便考察什么性質(zhì)的力驅(qū)動兩者的相互作用，結(jié)合自由能各成分的貢獻(xiàn)均列在表1中.

據(jù)表1，抑制劑BEG 與HIV-1蛋白酶的結(jié)合自由能是-22.25kcal／mol，表明抑制劑能與HIV-1蛋白酶有比較強(qiáng)的結(jié)合能力.范德華作用能（VDW）是-61.41kcal／mol，這個(gè)作用成分促進(jìn)了抑制劑與蛋白的結(jié)合.與溶劑可及表面積相對應(yīng)的非極性成分相互作用能 （PBSUR）是 -6.92kcal／mol，這個(gè)作用成分也為抑制劑的結(jié)合提供了較小的有利貢獻(xiàn).雖然抑制劑與HIV-1蛋白酶的靜電相互作用 （ELE）是-56.64kcal／mol，它也有利于抑制劑的結(jié)合，但是這種有力的作用成分被更強(qiáng)的極性溶解自由能所中和.從表一中可以看到熵變對對結(jié)合自由能 （TSTOT）的貢獻(xiàn)消弱了抑制劑與蛋白酶的結(jié)合.在有利的兩個(gè)相互作用成分中，范德華相互作用能是非極性相互作用的10倍.因此得出結(jié)論范德華相互作用能驅(qū)動了抑制劑BEG 與HIV-1蛋白酶的結(jié)合.

ELE和VDW 分別表示靜電作用能和范德華作用能，INT表示由鍵伸、鍵角彎折和二面角扭轉(zhuǎn)貢獻(xiàn)的內(nèi)能，Gas=ELE+VDW+INT；PBSUR 和PBCAL分別是非極性溶劑化能和極性溶劑化能，PBSOL溶解自由能且PBSOL=PBSUR+PBCAL；PBELE=ELE+PBCAL；PBTOT=PBSOL+GAS；TSTRA，TSTOT 和TSVIB分別表示體系的平動、轉(zhuǎn)動和振動對熵的貢獻(xiàn)；TSTOT=TSTRA+TSTOT+TSVIB；ΔGbind=PBTOT-TSTOT.

3.3 分子對接計(jì)算結(jié)合自由能

分子對接是一種基于經(jīng)典力場和打分函數(shù)評估抑制劑與蛋白質(zhì)的結(jié)合，能較好地探測藥物的靶標(biāo)，合理地預(yù)測和分析抑制劑與蛋白質(zhì)的結(jié)合模式，能為藥物的合理化設(shè)計(jì)提供理論指導(dǎo).這一部分我們將利用分子對接計(jì)算BEG 與HIV-1蛋白酶的結(jié)合自由能，以便能評估抑制劑與HIV-1蛋白酶的結(jié)合效果.

分子對接的計(jì)算結(jié)果如表2顯示，其中計(jì)算得到的抑制常數(shù)KI是5.82nM，分子間的相互作用能，內(nèi)能和扭轉(zhuǎn)能分別是-14.93，-0.88 和3.84 kcal／mol，未成鍵的排斥能為-0.74kcal／mol.依據(jù)表2可以看到BEG 與HIV-1蛋白酶的結(jié)合自由能為-11.23kcal／mol，這與實(shí)驗(yàn)值的結(jié)果 （-12.29 kcal／mol）［20］吻合的非常好.這個(gè)計(jì)算結(jié)果表明分子對接軟件的所使用的力場比Amber 12力場更適合于BEG-HIV-1蛋白酶系統(tǒng).

表2 分子對接計(jì)算的結(jié)合自由能 （kcal·mol-1 ）Table 2 Binding free energies calculated using the molecular docking method

3.4 抑制劑BEG 與HIV-1蛋白酶結(jié)合模式的分析

為了能更好地闡明抑制劑BEG 與HIV-1 蛋白酶的結(jié)合模式，我們采用MM／GBSA 方法計(jì)算抑制劑-殘基相互作用譜，如圖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 產(chǎn)生較強(qiáng)相互作用的相關(guān)殘基Fig.5 Key residues in HIV-1protease-BEG complex

據(jù)圖4，A 鏈中的殘基Leu23、Val82 和Ile84與BEG 的相互作用能分別為-1.23、-1.56和-2.08kcal／mol，從圖5的分子結(jié)構(gòu)看，這些相互作用能極好的吻合了三個(gè)殘基的烷基與BEG 苯環(huán)間的CH-π相互作用.殘基Gly49與BEG 的相互作用能為-2.08kcal／mol，從圖5 中的結(jié)構(gòu)看，這個(gè)能量主要來源于Gly49的羰基與苯環(huán)和五元環(huán)間的CH-O 相互作用.依據(jù)圖5 可以看到，Ile50的烷基與BEG 的苯環(huán)相互靠近，這形成了比較強(qiáng)烈的CH-π 的相互作用，該能量值為-3.28kcal／mol.圖4表明B鏈中的Asp25′與BEG有一個(gè)1.02kcal／mol的正的相互作用能，該能量主要因?yàn)锳sp25′的兩個(gè)帶負(fù)電荷的氧原子與BEG骨架上的兩個(gè)氧原子相距較近，從而產(chǎn)生了強(qiáng)烈的排斥作用.Gly27′與BEG 的相互作用能為-1.96kcal／mol，該能量主要來源于Gly27′與抑制劑BEG 間的CH-O.從圖5 的結(jié)構(gòu)中可以觀察到，B鏈中殘基Ala28′的烷基與BEG 的苯環(huán)相鄰近，因此能形成CH-π相互作用，從圖4獲知該相互作用能量值大約為-1.49kcal／mol.依據(jù)HIV-1蛋白酶的C2 對稱性，B 鏈中的三個(gè)殘基Gly49′、Ile50′和Ile84′與抑制劑BEG的相互作用分析類似于A 鏈中相對應(yīng)的三個(gè)殘基.總之，從上面的分析可以看到CH-π和CH-O 相互作用驅(qū)動了抑制劑BEG與HIV-1蛋白酶的結(jié)合.

4 結(jié) 論

本文用分子對接方法和MM／PBSA 方法計(jì)算抑制劑BEG 與HIV-1蛋白酶的結(jié)合自由能，結(jié)果表明范德華相互作用主宰了BEG 與HIV-1蛋白酶的相互作用.用基于殘基的相互作用能分解方法計(jì)算了抑制劑-殘基的相互作用，結(jié)果表明CH-π和CH-O 相互作用驅(qū)動了抑制劑BEG 與HIV-1蛋白酶的結(jié)合.我們期望這個(gè)研究能為艾滋病治療藥物的合理化設(shè)計(jì)提供理論上的指導(dǎo).

［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.（原子與分子物理學(xué)報(bào)），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，F(xiàn)eng 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.