王加磊，孟現(xiàn)美，張慶剛，張懌慈，張少龍

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

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

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

圖1 HIV－1蛋白酶的結(jié)構(gòu)Fig.1 Crystal structure of HIV－1 protease

1 理論與計算方法

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

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

GPU計算是一種大規(guī)模并行計算的新方法，該方法通過芯片內(nèi)數(shù)百個甚至上千個處理核心的通信和協(xié)作，以最高可超過傳統(tǒng)方法百倍的速度解決復(fù)雜的計算問題。

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

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

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

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

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

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

2 計算過程

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

3 結(jié)果與分析

3.1 動力學(xué)軌跡分析

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

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

圖2 分子片段4DX在蛋白酶表面Exo位結(jié)合示意圖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 異位抑制劑體系擋板區(qū)的運動Fig.4 Illustation of flap motion in allosteric inhibitor system

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

圖6為兩個體系擋板區(qū)骨架碳原子相對于初始結(jié)構(gòu)的RMSD隨時間變化示意圖?？梢钥闯霎愇灰种苿w系中的擋板區(qū)在整個模擬時間段無明顯起伏，保持在0.8 ?上下浮動，較為穩(wěn)定。活性位抑制劑體系中，擋板區(qū)在整個模擬時間段內(nèi)起伏變化劇烈，比較不穩(wěn)定。

圖7為兩個體系抑制劑相對于初始結(jié)構(gòu)的RMSD隨時間變化示意圖?？梢钥闯鲈陂_始的30 ns內(nèi)，由于兩體系未達到平衡狀態(tài)，RMSD變化劇烈。達到平衡后異位抑制劑體系的RMSD值整體低于活性位抑制劑體系，因此4DX與Exo位的結(jié)合有利于抑制劑在活性位點與蛋白質(zhì)結(jié)合，異位抑制劑體系穩(wěn)定性好。

圖5 分子動力學(xué)過程中骨架碳原子相對初始結(jié)構(gòu)的均方根偏差隨時間變化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)是衡量蛋白質(zhì)某個殘基在分子動力學(xué)過程中偏離其平均結(jié)構(gòu)的指標(biāo)，反映了蛋白質(zhì)殘基的靈活性和柔性度。B因子越大，說明相應(yīng)殘基靈活性大﹑柔性度高［15］。

我們重點對Exo位和擋板區(qū)氨基酸殘基B因子進行分析，擋板區(qū)由序列43/43′～58/58′組成的肽段構(gòu)成，Exo 位由序列 14/14′～18/18′、38/38′～42/42′和59/59′～65/65′構(gòu)成。它們在圖 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對稱的同質(zhì)二聚體，其Exo位、flap都為兩處，且C2對稱。本研究中4DX結(jié)合Exo－1位。圖8看出異位抑制劑體系中肽鏈1比肽鏈2穩(wěn)定;兩體系比較看出異位抑制劑體系的Exo－1位和除48號殘基外的flap－1的B因子比活性位抑制劑體系小，這些殘基波動幅度小、柔性弱、穩(wěn)定;異位抑制劑體系鏈2上Exo－2位無4DX結(jié)合，其B因子比活性位抑制劑體系大，因此殘基柔性大，不穩(wěn)定;flap－2的B因子明顯比活性位抑制劑體系小，因此殘基柔性弱，穩(wěn)定。上述分析看出未結(jié)合4DX的Exo－2位柔性大，不穩(wěn)定，而4DX對整個擋板區(qū)和與其結(jié)合的Exo－1位影響比較明顯，這些殘基活動范圍明顯減少，剛性更強，更穩(wěn)定，使活性位抑制劑不易脫落。

3.3 結(jié)合自由能分析

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

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

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

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

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

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

4 結(jié)語

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

［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］時術(shù)華，扈國棟，陳建中，等.分子動力學(xué)模擬研究質(zhì)子化態(tài)在HIV－1 protease－Indinavir復(fù)合物中的作用［J］.化學(xué)學(xué)報，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］侯廷軍，徐筱杰.基于分子動力學(xué)模擬和連續(xù)介質(zhì)模型的自由能計算方法［J］.化學(xué)進展，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蛋白酶與抑制劑結(jié)合中橋接水分子作用的QM/MM研究［J］.山東科學(xué)，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］孫德福，王加磊，張少龍，等.艾滋病毒蛋白酶異位抑制劑體系的分子動力學(xué)研究［J］.山東科學(xué)，2012，25(2):20－25.