HIV－1蛋白酶異位抑制劑體系的長時間分子動力學模擬

王加磊，孟現美，張慶剛，張懌慈，張少龍

(山東師范大學物理與電子科學學院，山東 濟南 250014)

HIV－1蛋白酶(HIV－1 protease)是艾滋病病毒生命過程必需的一種酶，蛋白酶抑制劑藥物作為高效抗逆轉錄病毒療法(HAART，俗稱雞尾酒療法)的重要組成部分已經在幫助艾滋病人存活方面發揮了很大作用［1］。美國食品藥品監督管理局(Food and Drug Administration，FDA)目前批準的9種蛋白酶藥物都是以活性位(active site)為靶標進行抑制的。但是，艾滋病毒已進化出對現有藥物具有耐藥性的變異形式，美國超過10%的新艾滋病病例是由病毒的耐藥株引起的，因此設計新機制的蛋白酶藥物迫在眉睫［2］。A.Perryman等［3］2010年發現了HIV－1蛋白酶上兩個在活性位之外的新結合位點(allosteric binding sites，異位結合位點):Exo位和Outside/top of flap位(圖1)，以及結合在這兩個異位結合位點上的化合物2－甲基環己醇(4DX)和6－吲哚甲酸，它們既能結合到新位點又能穩定活性位抑制劑的構象，雖然尚未表現出抑制作用，但可以成為開發針對耐藥性艾滋病病毒新藥的出發點。這一發現奠定了開發一類新型抗艾滋病毒藥物的基礎，沿此路線設計的新藥可以加強現有療法的療效，治療這種疾病的耐藥類型，減慢病毒的耐藥性進化。

HIV－I蛋白酶的結構如圖1所示，它由各包含99個氨基酸殘基的兩條肽鏈組成，活性位由Asp25/25′－Thr26/26′－Gly27/27′殘基組成，Asp25/25′在催化中起關鍵性的作用［4］，擋板區(flap)位于活性位上面并將活性部位遮蓋。flap對抑制劑或底物的進出起很大作用。

圖1 HIV－1蛋白酶的結構Fig.1 Crystal structure of HIV－1 protease

1 理論與計算方法

分子動力學模擬可以得到體系中各原子運動隨時間變化的細節，是理解蛋白質－配體相互作用的有力工具。勢函數不準確(力場問題)和模擬時間有限是分子動力學模擬成功應用受到的主要限制。

Amber ff12SB力場是目前最新的改進的非極化力場，它在ff99SB基礎上加入新的主鏈和側鏈扭轉勢，主鏈骨架扭轉勢改進后更精確，側鏈的改進是通過擬合更精確的量子從頭算數據獲得的，與實驗數據符合得比較好［5］。

GPU計算是一種大規模并行計算的新方法，該方法通過芯片內數百個甚至上千個處理核心的通信和協作，以最高可超過傳統方法百倍的速度解決復雜的計算問題。

本研究應用最新的Amber ff12SB力場在NVIDIA GTX670 GPU上對HIV－1蛋白酶異位抑制劑體系進行了長時間分子動力學模擬計算，得到了比以前更精確的結果。

結合自由能是蛋白質與抑制劑結合的重要指標，結合自由能的正負反映了蛋白質與抑制劑結合的可能性，結合自由能的大小反映了結合緊密程度［6］。MM－PB/GBSA(molecular mechanics Poisson－Boltzmann/Generalized Born surface area)是計算溶液環境中結合自由能的流行方法。受體和配體的結合自由能為:

MM－PB/GBSA方法把各類的自由能(complex，ligand，receptor)分解為氣相能量焓，溶解自由能和熵項［7］:

其中Ebat為成鍵的鍵、角和扭轉角項能的和，Evdw和Ecoul為非成鍵的范德華能和靜電相互作用能，Gsolv，polar為溶解自由能的極化項，Gsolv，nonpolar為溶解自由能的非極化項［8］。Ebat、Evdw、Ecoul之和為完整的氣相力場能，即 MM 部分。

溶解自由能極化項Gsolv，polar可通過Poisson－Boltzmann(PB)方程求解，也可用廣義波恩(GB)模型求解，溶解自由能非極化項Gsolv，nonpolar通常用經驗公式得到:

其中，SASA是溶劑可接觸表面積(solvent－accessible surface)。另外，溶解自由能非極化項Gsolv，nonpolar在更精確的方法下可以分解為排斥項和吸引項［9］。

2 計算過程

分子動力學模擬使用Amber12和AmberTools12軟件包。初始構象取自蛋白質庫(PDB ID:3KF0)，PDB中包含 HIV－1 PR、3TL、4DX、DMSO 和 BME，其中二甲基亞砜(Dimethyl Sulfoxide，DMSO)和 β－巰基乙醇(Beta－Mercaptoethano，BME)在HIV－1蛋白酶的表達、提純和結晶中作為有機溶劑起溶解和稀釋作用［3］，在本研究中予以忽略;提取PDB中HIV－1 PR+3TL+4DX作為異位抑制劑體系，提取HIV－1 PR+3TL作為活性位抑制劑體系;由于結晶水在HIV－1蛋白酶與抑制劑的結合中起了非常重要的作用［10］，PDB中的全部結晶水予以保留。將得到的結構在AmberTools12的LEaP模塊中補全缺失的氫原子并添加ff12SB力場參數。抑制劑3TL和分子片段4DX的部分電荷用sqm模塊中的AM1－BCC方法計算。兩個體系周圍加有9個氯離子以保持體系中性。使用SHAKE算法進行鍵約束計算［11］。將兩個體系質心置于截角八面體盒子中心，盒子壁到溶質原子的距離為10 ?，加入水分子，水分子采用TIP3P模型。接著采用Amber12的PMEMD模塊對兩個體系進行能量最小化［12］，能量最小化先是采用25000步最陡下降法，繼而是25000步共軛梯度法，然后用500000步將兩個體系溫度從0 K升到300 K，并在此溫度下進行2000000步動力學平衡。最后，對兩個體系分別進行100 ns的分子動力學模擬，運動方程積分步長(時間步長)2 fs，每隔1 ps記錄一次坐標文件，用于后續的分析。

3 結果與分析

3.1 動力學軌跡分析

X射線晶體衍射實驗發現分子片段2－甲基環已醇(4DX)結合在蛋白酶的非活性位點Exo位，Exo是由一對環區(殘基14/14′～18/18′和 38/38′～42/42′)和一 β 鏈(殘基 59/59′～65/65′)形成如圖 2所示的蛋白酶表面凹槽［3］。2－甲基環已醇的羥基與殘基Gly17形成氫鍵，它的甲基環已基(the methycylclohexyl group)與Lys14和Leu63存在疏水性結合。Exo位的38/38′～42/42′和 59/59′～65/65′序列緊鄰蛋白酶的擋板區(flap)，flap的變化活動與Exo位為反相關關系，當flap打開時，Exo位凹槽會被壓縮，反之flap關閉時Exo位凹槽就會擴大，類似于杠桿系統［3］。

動力學運動軌跡由不同時刻體系中所有原子的坐標組成［13］，其最直觀的呈現是將其用VMD軟件制作成電影。觀察軌跡文件的電影可見在100 ns模擬過程中4DX一直結合在Exo位(圖2)。由于組成Exo位的Lys14、Leu63和Gly17殘基與4DX存在疏水結合和氫鍵，這三個殘基構成Exo位凹槽整體框架(圖3)，本研究以這三個殘基α碳原子構成的三角形邊長距離來衡量Exo位凹槽大小，計算50 ns時刻時異位抑制劑體系三角形邊長為 9.09 ?、6.23 ?、10.18 ?，活性位抑制劑體系為 8.40 ?、6.15 ?、9.58 ?，比較可知異位抑制劑體系的三角形邊長比活性位抑制劑體系長，Exo位凹槽開口更大，同時比較兩體系得異位抑制劑體系中擋板區活動方向如圖4所示，即擋板區關閉更牢固，這說明flap活動與Exo位為反相關關系。

圖2 分子片段4DX在蛋白酶表面Exo位結合示意圖Fig.2 Illustration of PR surface Exo site binded molecular fragment 4DX

圖3 Exo位的Lys14、Leu63、Gly17殘基示意圖Fig.3 Illustration of residues Lys14，Leu63 and Gly17 at Exo site

圖4 異位抑制劑體系擋板區的運動Fig.4 Illustation of flap motion in allosteric inhibitor system

系統是否達到動態平衡通常用體系的溫度﹑密度﹑總能量和均方根偏差 (Root－mean－square deviation，RMSD)來衡量［14］。兩個體系相對初始結構的RMSD如圖5所示，可以看出在100 ns分子動力學模擬過程中兩個體系之間蛋白質骨架碳原子的均方根波動差異明顯，異位抑制劑體系的RMSD從開始到10 ns緩慢上升至1.0 ?，并在隨后的 90 ns內保持在 1.0 ? 上下浮動。活性位抑制劑體系的RMSD在開始時快速躍升到1.2 ?，并在后面模擬時間段內保持在1.2 ?上下浮動。兩個體系在30 ns以后達到平衡。整體上觀察，異位抑制劑體系的RMSD比活性位抑制劑體系要小，這說明異位抑制劑體系更穩定。

圖6為兩個體系擋板區骨架碳原子相對于初始結構的RMSD隨時間變化示意圖。可以看出異位抑制劑體系中的擋板區在整個模擬時間段無明顯起伏，保持在0.8 ?上下浮動，較為穩定。活性位抑制劑體系中，擋板區在整個模擬時間段內起伏變化劇烈，比較不穩定。

圖7為兩個體系抑制劑相對于初始結構的RMSD隨時間變化示意圖。可以看出在開始的30 ns內，由于兩體系未達到平衡狀態，RMSD變化劇烈。達到平衡后異位抑制劑體系的RMSD值整體低于活性位抑制劑體系，因此4DX與Exo位的結合有利于抑制劑在活性位點與蛋白質結合，異位抑制劑體系穩定性好。

圖5 分子動力學過程中骨架碳原子相對初始結構的均方根偏差隨時間變化Fig.5 The curves of time vs.RMSD of the backbone C atoms relative to their initial structures in MD

3.2 蛋白酶殘基B因子分析

B因子(B－factor)是衡量蛋白質某個殘基在分子動力學過程中偏離其平均結構的指標，反映了蛋白質殘基的靈活性和柔性度。B因子越大，說明相應殘基靈活性大﹑柔性度高［15］。

我們重點對Exo位和擋板區氨基酸殘基B因子進行分析，擋板區由序列43/43′～58/58′組成的肽段構成，Exo 位由序列 14/14′～18/18′、38/38′～42/42′和59/59′～65/65′構成。它們在圖 8 表示為:Exo－1 位(14～18、38 ～42、59 ～65);Exo－2 位(113 ～ 117、137 ～ 141、158 ～ 164);flap－1(42 ～ 58);flap－2(141 ～ 157)。HIV－1蛋白酶因為是C2對稱的同質二聚體，其Exo位、flap都為兩處，且C2對稱。本研究中4DX結合Exo－1位。圖8看出異位抑制劑體系中肽鏈1比肽鏈2穩定;兩體系比較看出異位抑制劑體系的Exo－1位和除48號殘基外的flap－1的B因子比活性位抑制劑體系小，這些殘基波動幅度小、柔性弱、穩定;異位抑制劑體系鏈2上Exo－2位無4DX結合，其B因子比活性位抑制劑體系大，因此殘基柔性大，不穩定;flap－2的B因子明顯比活性位抑制劑體系小，因此殘基柔性弱，穩定。上述分析看出未結合4DX的Exo－2位柔性大，不穩定，而4DX對整個擋板區和與其結合的Exo－1位影響比較明顯，這些殘基活動范圍明顯減少，剛性更強，更穩定，使活性位抑制劑不易脫落。

3.3 結合自由能分析

圖8 兩體系殘基的B因子Fig.8 B－factor graph of residues of the two systems

本研究用MM－PB/GBSA方法分別計算了兩個體系中活性位抑制劑TL－3與HIV－1蛋白酶的結合自由能，結果如表1、2所示。表中ΔGELE表示分子力場中靜電能，ΔGVDW表示分子力場中范德華能，ΔGGAS表示氣相能，ΔGPBCAL和ΔGPBSUR表示PB計算方法中極性溶解能和非極性溶解能，ΔGPBSOL表示PB計算方法中極性和非極性的溶解能總和，ΔGGB和ΔGGBSUR表示GB計算方法中極性的溶解能和非極性的溶解能，ΔGGBSOL表示GB計算方法中極性和非極性溶解能總和，ΔGPBTOT和ΔGGBTOT分別表示PB、GB模型中最終的結合自由能［12］。其中，氣相能為分子力場中靜電能、范德華能和內能三者之和。由于分子力場中內能在單軌跡近似方法中總是為 0［6］，表中未列出此項。

從表1、2可以看出，MM－PB/GBSA計算中異位抑制劑體系中蛋白酶與抑制劑的結合自由能為－85.78 kcal/mol和 － 105.89 kcal/mol，活性位抑制劑體系中結合自由能為 － 79.45 kcal/mol和－92.03 kcal/mol，異位抑制劑體系的結合自由能數值大小比活性位抑制劑體系大，這說明異位抑制劑體系中蛋白酶與抑制劑的結合更牢固，即4DX與蛋白酶和活性位抑制劑TL－3形成的復合結構更穩定。

從表1、表2還可以看出，對結合自由能貢獻較大的是靜電能、范德華能和極化溶解能。靜電能、范德華能和非極化溶解能為負值，有利于抑制劑與蛋白酶的結合，極化溶解能為正值，不利于抑制劑與蛋白酶的結合。

表1 MM－PBSA計算結合自由能(kcal/mol)Table 1 Binding free energy of MM－PBSA(kcal/mol)

表2 MM－GBSA計算結合自由能(kcal/mol)Table 2 Binding free energy of MM－GBSA(kcal/mol)

4 結語

本研究對HIV－1蛋白酶的活性位抑制劑體系和異位抑制劑體系分別進行了100 ns的分子動力學模擬，并用MM－PB/GBSA方法計算了兩體系的活性位抑制劑TL－3與HIV－1蛋白酶結合自由能。MM－PBSA計算方法中異位抑制劑體系中蛋白酶與抑制劑的結合自由能為－85.78 kcal/mol，活性位抑制劑體系中結合自由能為－79.45 kcal/mol;MM－GBSA計算方法中異位抑制劑體系中蛋白酶與抑制劑的結合自由能為－105.89 kcal/mol，活性位抑制劑體系中結合自由能為－92.03 kcal/mol。4DX結合在HIV－1蛋白酶表面Exo位凹槽，擴大了凹槽體積，使擋板區關閉更牢固，同時，兩體系擋板區B因子的對比表明異位抑制劑體系擋板區活動范圍減少，剛性變強，更穩定。這些分子動力學結果證實異位抑制劑體系使活性位抑制劑TL－3與HIV－1蛋白酶的結合更牢固，因此體系更穩定。這為開發艾滋病新型抑制劑提供了有價值的線索。

［1］JASKOLSKI M，TOMASSELLI A G，SAWYER T K，et al.Structure at 2.5－A resolution of chemically synthesized human immunodeficiency virus type 1 protease complexed with a hydroxyethylene－based inhibitor［J］.Biochemistry，1991，30(6):1600－1609.

［2］WATKINS T，RESCH W，IRLBECK D，et al.Selection of high－level resistance to human immunodeficiency virus type 1 protease inhibitors［J］.Antimicrob Agents Chemother，2003，47(2):759 －769.

［3］PERRYMAN A L，ZHANG Q，SOUTTER H H，et al.Fragment－based screen against HIV protease［J］.Chem Biol Drug Des，2010march，75(3):257 －268.

［4］時術華，扈國棟，陳建中，等.分子動力學模擬研究質子化態在HIV－1 protease－Indinavir復合物中的作用［J］.化學學報，2009，67(24):2791 －2797.

［5］LINDORFF－LARSEN K，PIANA S，PALMO K，et al.Improved side－chain torsion potentials for the Amber ff99SB protein force field［J］.Proteins，2010，78(8):1950 －1958.

［6］侯廷軍，徐筱杰.基于分子動力學模擬和連續介質模型的自由能計算方法［J］.化學進展，2004，16(2):153－158.

［7］MILLER B R III，McGee T D Jr，SWAILS J M，et al.MMPBSA.py:An efficient program for end－state free energy calculations［J］.J Chem Theory Comput，2012，8(9):3314 –3321.

［8］CASE D A，CHEATHAM T E Ⅲ，DARDEN T，et al.The amber biomolecular simulation programs［J］.J Comput Chem，2005，26(16):1668－1688.

［9］TAN C H，TAN Y H，LUO R.Implicit nonpolar solvent models［J］.J Phys Chem B，2007，111(42):12263 －12274.

［10］張玉新，張少龍，張慶剛.HIV－1蛋白酶與抑制劑結合中橋接水分子作用的QM/MM研究［J］.山東科學，2010，23(1):1－5.

［11］ADCOCK S A，McCAMMON J A.Molecular dynamics:Survey of methods for simulating the activity of proteins［J］.Chem Rev，2006，106(5):1589－1615.

［12］PEARLMAN D A，CASE D A，CALDWELL J W，et al.AMBER，a package of computer programs for applying molecular mechanics，normal mode analysis，molecular dynamics and free energy calculations to simulate the structural and energetic properties of molecules［J］.Computer Physics Communications，1995，91(1 －3):1 －41.

［13］MOORE J P，STEVENSON M.New targets for inhibitors of HIV－1 replication［J］.Nat Rev Mol Cell Biol，2000，1(1):40 －49.

［14］de CLERCQ E.The design of drugs for HIV and HCV［J］.Nat Rev Drug Discov，2007，6(12):1001 －1018.

［15］孫德福，王加磊，張少龍，等.艾滋病毒蛋白酶異位抑制劑體系的分子動力學研究［J］.山東科學，2012，25(2):20－25.