HTLV—1蛋白酶與印地那韋識別的分子動力學模擬

常珊等

摘要：對HTLV-1蛋白酶和抑制劑印地那韋的復合物進行了11 ns的分子動力學模擬，從分子整體、氫鍵、結合自由能及運動性等方面進行了分析。結果表明，動力學模擬使得體系更為松弛，但大部分氫鍵仍得到保持；疏水相互作用有利于HTLV-1蛋白酶和印地那韋的結合，關鍵殘基主要分布在R10、L30～V39、V56～F67和W98～I1004個區域；結合自由能預測值與試驗數據吻合較好。印地那韋結合僅能部分抑制HTLV-1蛋白酶的Flap、top-site以及兩個exo-site位點柔性，解釋了印地那韋的HIV-1 PR抑制活性為何略大于HTLV-1 PR。該研究為基于受體結構的抗HTLV-1藥物研發打下了一定的結構和理論基礎。

關鍵詞：HTLV-1蛋白酶；印地那韋；分子動力學模擬；分子識別；結合自由能

中圖分類號：O641.3；R969.1 文獻標識碼：A 文章編號：0439-8114（2015）12-3034-07

DOI：10.14088/j.cnki.issn0439-8114.2015.12.062

Molecular Dynamics Simulation of the Recognition between HTLV-1

Protease and Indinavir

CHANG Shan1，LIANG Li2，LIU Wei2，HU Jian-ping2，3，GOU Xiao-jun2

（1. Institute of Bioinformatics and Medical Engineering， Jiangsu University of Technology， Changzhou 213001， Jiangsu China；

2. Key Laboratory of Medicinal and Edible Plants Resources Development， Chengdu University， Chengdu 610106， China；

3. College of Chemistry， Leshan Normal University， Leshan 614004， Sichuan China）

Abstract： One 11 ns molecular dynamics simulation （MD） for the compocind system of Human T-cell leukemia virus type 1（HTLV-1） protease （PR） and its inhibitor indinavir was performed. The mutual recognition of the two molecules was analyzed from exploring the whole molecule，hydrogen bond，binding free energy and mobility，etc.After dynamics simulation，the system keeps more relaxed compared with the crystal structure， and most of the hydrogen bonds were still maintained. Both energy decomposition and free energy calculation showed that，hydrophobic interaction facilitates the association of HTLV-1 PR and indinavir，and key residues were mainly distributed in the four regions including R10，L30～V39，V56～F67 and W98～I100.Additionally，binding free energy prediction was in good agreement with the experimental data. Indinavir can only partially attenuate the flexibility of the flap，top-site and two exo-site in the PR system，which investigated why indinavir possess more HIV-1 PR inhibition ability than that of HTLV-1 PR. The study will provide some structural and theoretical basis for anti-HTLV-1 drug research based on the receptor structure.

Key words： HTLV-1 protease； Indinavir； molecular dynamics simulation； molecular recognition； binding free energy

I型人類T細胞白血病病毒（Human T-cell leukemia virus type 1，HTLV-1）參與了眾多種類癌癥，如成人T細胞白血病、熱帶痙攣性脊髓病以及HTLV-1相關的皮炎等[1，2]。目前，世界上約2 000萬人感染了該病毒，臨床治療策略主要是采用細胞毒性化療及干擾素、齊多呋定用藥，開展抗HTLV-1的新型藥物研發比較迫切[3]。

HTLV-1蛋白酶（Protease，PR）是同源蛋白二聚體，每個鏈由125個氨基酸殘基，屬于天冬氨酸蛋白酶家族[4，5]。HIV-1蛋白酶被認定為抗AIDS藥物研發的一線靶點已有30年歷史，目前FDA已經批準10多種HIV-1蛋白酶抑制劑藥物上市[6]。而HTLV-1 PR和HIV-1 PR的序列同源性為28%，其中抑制劑結合區同源性更高，達到了45%。兩類PR盡管具有相似的整體折疊，但是有不同的抑制劑結合特異性[7]。大部分HIV-1 PR抑制劑沒有體現出抑制HTLV-1 PR的潛力，目前僅有HIV-1蛋白酶抑制劑印地那韋能夠在較低的微摩爾級層次上抑制HTLV-1 PR活性，抑制常數為3.5 μm[8]。

圖1給出了印地那韋的分子結構，分子式為C36H47N5O4，由主鏈（C10-C11-C12-C13-C21-N4）和側鏈部分（比如吡啶基、哌嗪基、叔丁基，苯基，茚基）組成。綜合考慮HTLV-1 PR和HIV-1 PR的結構類似性和印地那韋的雙重抑制活性，本研究利用分子動力模擬方法研究了印地那韋和HTLV-1 PR的相互作用，將有助于后續的印地那韋類抗HTLV-1 PR抑制劑的設計。

1 模擬方法

1.1 分子動力學模擬

基于HTLV-1蛋白酶結合抑制劑印地那韋的復合物晶體[9]結構，僅保留A鏈和B鏈接的氨基酸部分以及底物印地那韋，以此作為分子動力學（Molecular dynamics，MD）模擬的初始構象。MD模擬應用AMBER程序[10]，采用的力場是AMBER力場，蛋白質的力場參數均基于試驗值擬合[11]，模擬溫度為300 K，溶劑采用TIP3P水模型[12]。在溶質印地那韋外圍加上1.0 nm的去頭八面體水盒子，總共加入9 913個水分子，含水后體系的總原子數為33 453。在MD模擬之前，對體系進行兩次能量優化：①約束溶質[約束力常數為2.09×103 kJ/（mol·nm2）]，用最陡下降法優化5 000步，再用共軛梯度法優化5 000步；②去約束后，再進行5 000 步最陡下降法優化和5 000步共軛梯度法優化，收斂條件為能量梯度小于4.18×10-4 kJ/（mol·nm）。

MD 模擬分為兩步：首先進行1 ns約束溶質的MD 模擬[（約束力常數為41.8 kJ/（mol·nm2）]，溫度從0 K逐步升高到300 K；再進行10 ns的無約束恒溫MD模擬，用VMD 圖像顯示軟件實時跟蹤體系的構象。采用SHAKE 算法[13]約束鍵長，MD模擬的積分步長設為2 fs，非鍵相互作用截斷半徑設為1.2 nm，每隔1 ps采集1次構象，共采集11 000個構象。

1.2 MM-PBSA自由能計算

用MM-PBSA（Molecular mechanics/Poisson-Boltzmann surface area）方法[14，15]計算了HTLV-1 PR和底物印地那韋的絕對結合自由能，該方法通過從復合物的MD模擬軌跡中每隔一定的時間間隔提取出體系結構，計算出平均結合自由能，計算公式為：

ΔGbind=ΔH-TΔS≈ΔEMM+ΔGsol-TΔS

式中，ΔGbind是結合自由能，ΔH是焓變，TΔS是絕對溫度和熵變的乘積；ΔEMM是真空下分子內能差，包括分子非極性能差（VDWIN）和分子靜電能（ELEIN）差；ΔGsol是溶劑化自由能，ΔEMM可通過分子力學方法計算得到，ΔGsol由極性溶劑化自由能差（ELEPB）和非極性溶劑化自由能（VDWPB）差組成，極性部分通過求解有限差分Poisson-Blotzmann（PB）方程求得[16-18]，非極性部分通過估算溶劑可接近表面積（Surface area，SA）[19]擬合得到。TΔS的計算采用了正則模方法[20]。

1.3 能量分解

能量分解采用MM/GBSA（Molecular mechanics/generalized Born surface area）方法[14，15]，將每個殘基的能量貢獻近似分為用分子力學（MM）方法計算的真空分子內能、用廣義波恩（GB）模型[21，22]計算的極性溶劑化能，以及用LCPO模型[23，24]計算的非極性溶劑化能，并且把能量分解到殘基的主鏈原子和側鏈原子上。其中，真空分子內能是由極性部分（靜電相互作用）和非極性部分（范德華相互作用）組成的。LCPO模型的核心是非極性溶劑化能與分子溶劑可接近表面積正相關，通過MM/GBSA能量分解，可以考察蛋白質中的主要殘基對于結合底物分子所做的貢獻。本研究基于HTLV-1蛋白酶與抑制劑印地那韋的復合物體系，MD 模擬平衡后6.5～11.0 ns 內的軌跡（共4.5 ns），每隔100 ps采樣一次，共取45個構象，通過計算獲得所有構象的平均能量分解結果。

1.4 構象成簇

對體系MD 模擬所獲得的11 000 個分子構象進行了成簇分析，分析采用MMTSB工具[25]，計算過程依據Daura等[26]報道的策略。逐個計算各個構象間的RMSD值，隨后建立RMSD矩陣（N×N，N為構象數）。設定一個RMSD域值，如體系兩個構象之間的RMSD小于該域值，則將它們歸為一簇，最后將每簇中能量最低構象作為該簇的代表性結構。本研究的成簇分析則是建立HTLV-1 PR-印地那韋復合物體系整體RMSD矩陣，即計算各個體系間的整體RMSD值，以分析基于整體構象的成簇分布，進而依據構象所占概率獲得分子整體所體現出的優勢構象。

2 結果與分析

2.1 體系的整體動力學

2.1.1 勢能及回轉半徑 圖2給出了HTLV-1蛋白酶與印地那韋復合物體系的勢能和回轉半徑（Radius of gyration）隨時間的變化。從圖2A可以看出，在約束溶質的MD模擬初級階段（即0～1 ns之間），勢能平均值和標準偏差分別是-405 811×105 kJ/mol和6 477×105 kJ/mol，偏差率為1.6%。在去約束的MD模擬階段（即1～11 ns之間），勢能維持在-（414 201±681）×105 kJ/mol，偏差率為0.2%。在長時間全含水MD模擬中，體系勢能維持穩定，模擬過程及軌跡可靠。從圖2B可以看出，HTLV-1蛋白酶結構最初保持約束狀態（0～1 ns之間），回轉半徑數值穩定，均值為1.895 nm；在去約束后的1.0～5.5 ns之間，回轉半徑漲落幅度變大，均值為1.899 nm；經過5.5～6.5.0ns的回轉半徑（均值1.913 nm）快速增加階段，蛋白質已經去掉了晶體中固有的大部分約束和碰撞，結構更松弛。在6.5～11 ns之間，回轉半徑維持在1.930 nm附近漲落。

2.1.2 RMSD分析 圖3給出了復合物RMSD隨時間的變化。從圖3可以看出，HTLV-1 PR在0～1 ns的約束MD模擬過程中，RMSD非常穩定，均值維持在0.01 nm；在1.0～5.5 ns區間對應的數值為（0.13±0.01） nm，5.5～6.5 ns為RMSD快速增加區間，在6.5～11.0 ns區間，對應的數值為（0.19±0.01） nm。A鏈和B鏈的RMSD走勢稍有差別，A鏈在1.0～5.5 ns區間的RMSD要低于B鏈，5.5 ns之后，二者走勢趨同。就印地那韋而言，從模擬開始到3.0 ns之間，RMSD迅速提高到0.3 nm，隨后逐漸保持平穩，數值為（0.22±0.02） nm，略大于蛋白質部分，可見小分子有較高的運動性。與圖2B所示的回轉半徑結果相比，蛋白質部分RMSD和回轉半徑數值快速增加的時間節點相同，均為5.5～6.5 ns之間，蛋白質RMSD走勢與回轉半徑有較強相關性。另外，蛋白質A鏈與B鏈的運動性略有不同，相對而言A鏈略高于B鏈。

2.1.3 RMSF分析 方均根漲落（Root mean square deviation，RMSF）是體系原子相對于平均構型的方均根漲落，反映了體系在整個模擬中原子的平均變化情況。RMSF數值大說明對應位置的柔性大，反之相反。圖4給出了HTLV-1蛋白酶二聚體的A鏈和B鏈的Cα原子的RMSF及柔性分布。從圖4A可以看出，HTLV-1蛋白酶的兩條鏈RMSF分布類似，這里將RMSF值大于0.4 nm和小于0.2 nm的區域分別定義為柔性區和剛性區。可見，柔性區有S46～N53，G60～F67和F80～T83；剛性區有A29～V39，A99～A105。圖4B直觀給出了體系柔性的區域分布，柔性大的區域都是遠離結合部位的無規卷曲或β片；剛性區則是直接參與識別抑制劑IND，也包括附近剛性較大的α螺旋區。另外，A鏈的柔性略高于B鏈。

在X-ray試驗中，蛋白質的B因子（B-factor）體現了晶體中原子電子密度的模糊度（Diffusion），模糊度實際反映了蛋白質分子在晶體中的構象狀態。B因子越高，模糊度越大，相應部位的構象柔性越大。將基于X-ray試驗獲得的B因子和MD模擬獲得的RMSF進行了相關性分析。HTLV-1蛋白酶A鏈和B鏈的Cα原子B因子及RMSF數據的相關性，相關性值分別為0.66和0.91（圖5A、圖5B）。在116個樣本數前提下，相關性非常高，可見蛋白質A鏈和B鏈的柔性分布類似，也表明用X-ray和MD方法來表述蛋白質柔性分布的有效性。蛋白質A鏈及B鏈的B因子和RMSF數值的相關性，相關系數分別為0.56和0.71（圖5C、圖5D）。綜上，基于RMSF分析得到的蛋白質柔性分布結果可信，模擬所得軌跡可以用于后續分子間識別的探討。

2.2 分子識別

2.2.1 氫鍵分析 根據圖3的RMSD分布，分析了7.0～9.0 ns之間的氫鍵分布。氫鍵采用幾何判據[27]，供體和受體之間距離小于0.35 nm，供體-氫-受體的夾角大于135°，氫鍵占有率表示在整個模擬過程中氫鍵構象出現的幾率。從表1可以看出，參與抑制劑印地那韋和蛋白酶識別的有6個氫鍵，氫鍵給體涉及到抑制劑的O2、O4、N1和N4原子，貫穿小分子的大部分結構。氫鍵是促進印地那韋與HTLV-1蛋白酶特異性識別的重要作用力。X-ray試驗結果[28]表明，R10、D32、G34、D36、A59、W98能與印地那韋形成氫鍵，在全含水MD優化后，大部分氫鍵得到了保持。

圖6給出了HTLV-1蛋白酶B鏈D32-OD2與印地那韋-O2-H18間最穩定氫鍵的距離及角度隨時間的變化。在7～9 ns間，受體-供體距離和受體-氫-供體角度維持穩定，漲落幅度較小，在2 000個構象中，共有1 889個構象均保持有印地那韋-O2-H18-D32-OD2氫鍵（占有率為94.45%）。

2.2.2 結合自由能計算 基于HTLV-1蛋白酶-印地那韋復合物7～9 ns的MD模擬軌跡，刪除水分子和4個平衡氯離子，然后從該軌跡中每隔200 ps取一次構象，共11個構象用于自由能計算，最終結果通過取平均得到。表2列出了各能量項對結合自由能的貢獻，用MM-PBSA方法計算得到的結合自由能模擬值-32.87 kJ/mol。印地那韋抑制HTLV-1蛋白酶的抑制常數Ki為3.5 μmol/L，相應的ΔGbind試驗值為-31.33 kJ/mol，自由能計算值和試驗值非常吻合。

在形成蛋白酶-印地那韋復合物時，體系的范德華相互作用（VDWIN）和非極性溶劑化能（VDWPB）都有利于底物的結合；體系的靜電相互作用則是不利于底物的結合，主要表現在盡管真空條件下分子靜電能（ELEIN）有利于底物的結合，但是極性溶劑化能（ELEPB）極其不利于二者識別，并直接反轉了ELEIN有利于結合的趨勢。體系熵效應（TS）是不利于底物結合的，這與分子間識別基本都是熵減效應一致。

2.2.3 能量分解及結合模式 為了研究HTLV-1蛋白酶與抑制劑印地那韋的相互作用，同時考慮到用GB模型計算靜電省時間的因素，用GBSA方法計算了蛋白酶結合印地那韋所涉及到的重要殘基能量貢獻。

表3給出了HTLV-1蛋白酶中與抑制劑印地那韋識別密切相關殘基的能量分解結果（能量小于-0.5 kJ/mol），表中的右上角撇號表示B鏈。從表3可以看出，利于HTLV-1蛋白酶與印地那韋結合的關鍵殘基主要分布在4個區域：R10、L30～V39、 V56～F67、W98～I100。其中A鏈和B鏈有利于結合印地那韋的關鍵殘基分別是10和11個。其中在A鏈和B鏈都出現的關鍵殘基有7個，具體為R10、G34、A35、V56、G58、W98和I100。盡管蛋白質A鏈和B鏈空間上中心對稱，小分子也結合在正中間，但是由于小分子結構并不對稱，具有一定的空間異構性，A鏈及B鏈對于結合印地那韋的貢獻類似，但不盡相同，該結果與之前RMSD（圖3）及RMSF（圖4）分析結果一致。B鏈的F67也是促進蛋白質與小分子識別的關鍵殘基，這個殘基在晶體結構[28]中并沒有被提到。

基于上面的能量分解結果，圖7給出了HTLV-1蛋白酶與印地那韋的識別模式。圖中藍色和綠色stick模型分別表示關鍵殘基，蛋白質整體用flat-ribbon模型表示。

2.3 分子的運動性

用MMTSB軟件包[25]對MD模擬軌跡11 000個構象進行了成簇分析。圖8給出了基于相互RMSD對比的復合物體系整體構象成簇分析結果，其中RMSD閾值設為0.13 nm。MD模擬構象共分為5簇，2.5 ns之前基本為第1簇，共2 460個構象，所占概率為22.36%；2.5～6.3 ns之間為第2簇，共3 484個構象，占比31.67%；第3簇則主要是出現在9.0 ns后，構象數和占比分別為1 447和13.15%；第4簇和第5簇的時間區間基本為6.0～9.0 ns之間。值得一提的是，第3簇、第4簇和第5簇構象有較為明顯的交叉躍遷現象，尤其第4簇和第5簇之間。

圖9給出了HITLV-1蛋白酶與抑制劑印地那韋復合物體系的每簇代表性構象的疊落結果。代表性結構指的是每個體系各簇中能量最低構象，即第1～5簇代表性構象代碼分別為1 754、4 220、10 998、8 410和7 745。圖9A給出了基于蛋白質主鏈Cα原子重疊的5個代表性構象疊落結果。蛋白質A，B鏈部分的Flap、top-site以及兩個exo-site位點疊落較差，相應區域的柔性非常大。藥物篩選研究試驗[29]表明，flap、top-site、exo-site等長loop區較高的柔性是抑制劑與蛋白酶結構不夠穩定的結果。印地那韋抑制劑與HTLV-1蛋白酶抑制劑的結合力要弱于與HIV-1蛋白酶的識別，這能夠部分解釋loop區柔性為何較高。

另外，在代表性構象中首先保證小分子哌嗪環（N1-C1-C2-N3-C8-C9）重疊，然后進行疊落操作，結果如圖9B所示。蛋白質相應區域都是中心對稱，由于印地那韋的不對稱導致其與蛋白酶結合和運動具有一定的復雜度。印地那韋的吡啶、哌嗪環和叔丁基主要與B鏈結合，運動幅度較小，這與之前的RMSD和RMSF分析一致；印環基則主要靠近A鏈，疊落效果最差，運動幅度也最高；而苯環基主要靠近B鏈的R10殘基，柔性居中。

3 結論

基于HTLV-1蛋白酶和抑制劑印地那韋的復合物長時間顯含水的分子動力學模擬，分析了二者的相互作用機制。體系的回轉半徑和方均根偏差分析表明，與晶體結構相比，HTLV-1蛋白酶-印地那韋復合物構象在水溶液中逐漸松弛，且A鏈柔性略高于B鏈。蛋白質A鏈及B鏈的B因子和RMSF數值有很高的相關性，相關系數分別達到0.56和0.71，表明柔性分析合理。柔性大的區域都是遠離結合部位的無規卷曲或β片；剛性區則是直接參與識別抑制劑印地那韋，也包括附近剛性較大的α螺旋區。能量分解表明，HTLV-1蛋白酶與印地那韋結合的關鍵殘基主要分布在R10、L30～V39、V56～F67和W98～I100等4個區域。用MM-PBSA方法獲得HTLV-1蛋白酶和抑制劑印地那韋的結合自由能計算值（-32.87 kJ/mol）與試驗值（-31.33 kJ/mol）吻合較好，其中，體系的范德華相互作用和非極性溶劑化能都有利于二者結合運動性分析表明，蛋白質A，B鏈部分的Flap、top-site以及兩個exo-site位點疊落較差，相應區域的柔性非常大；而印地那韋抑制劑的印環基則主要靠近柔性較大的A鏈，疊落效果最差，運動幅度也最高。模擬結果為基于HTLV-1蛋白酶受體結構的印地那韋類抑制劑設計有一定的指導意義。

參考文獻：

[1] LI M， LACO G S， JASKOLSKI M， et al. Crystal structure of human T cell leukemia virus protease， a novel target for anticancer drug design[J]. Proc Natl Acad Sci USA， 2005， 102： 18332-18337.

[2] MCGILL N K，VYAS J， SHIMAUCHI T， et al. HTLV-1-associated infective dermatitis： Updates on the pathogenesis[J]. Exp Dermatol， 2012，21：815-821.

[3] ISHIKAWA T. Current status of therapeutic approaches to adult T-cell leukemia[J]. Int J Hematol， 2003， 78： 304-311.

[4] WLODAWER A. Crystal structures of inhibitor complexes of human T-cell leukemia virus （HTLV-1） protease[J]. J Mol Biol， 2010， 401： 626-641.

[5] SHUKER S B， MARIANI V L， HERGER B E， et al. Understanding HTLV-1 protease[J]. Chem Biol，2003，10：373-380.

[6] MEHELLOU Y， DE CLERCQ E. Twenty-six years of anti-HIV drug discovery： Where do we stand and where do we go [J]. J Med Chem， 2010， 53： 521-538.

[7] KADAS J， WEBER I T， BAGOSSI P， et al. Narrow substrate specificity and sensitivity toward ligand-binding site mutations of human T-cell leukemia virus type 1 protease[J]. J Biol Chem， 2004，279： 27148-27157.

[8] BAGOSSI P， KADAS J， MIKLOSSY G， et al. Development of a microtiter plate fluorescent assay for inhibition studies on the HTLV-1 and HIV-1 proteinases[J]. J Virol Methods， 2004， 119： 87-93.

[9] MAREN K， HOLGER S， WIBKE E， et al. Structural basis for HTLV-1 protease inhibition by the HIV-1 protease inhibitor indinavir[J]. J Med Chem， 2014， 57： 6266-6272.

[10] WANG J M， CIEPLAK P， KOLLMAN P A. How well does a restrained electrostatic potential （RESP） model perform in calculating conformational energies of organic and biological molecules[J]. J Comput Chem， 2000， 21（12）： 1049-1074.

[11] WANG J M， WOLF R M， CALDWELL J W， et al. Development and testing of a general amber force field[J]. J Comput Chem， 2004， 25（9）： 1157-1174.

[12] JORGENSEN W L， CHANDRASEKHAR J， MADURA J D， et al. Comparison of simple potential functions for simulating liquid water [J]. J Chem Phys， 1983， 79 （2）： 926-935.

[13] RYCKAERT J P， CICCOTTI G， BERENDSEN H J C. Numerical integration of the cartesian equations of dynamics of n-alkanes[J]. J Comput Phys， 1977， 23（3）： 327-341.

[14] KOLLMAN P A， MASSOVA I， REYES C， et al. Calculating structures and free energies of complex molecules： Combining molecular mechanics and continuum models[J]. Acc Chem Res， 2000， 33： 889-897.

[15] WANG W， DONINI O， REYES C， et al. Biomolecular simulations： Recent developments in force fields， simulations of enzyme catalysis， protein-ligand， protein-protein， and protein-nucleic acid noncovalent interactions[J]. Annu Rev Biophys Biomol Struct， 2001， 30： 211-243.

[16] HONIG B， NICHOLLS A. Classical electrostatics in biology and chemistry[J]. Science，1995，268：1144-1149.

[17] HECHT J L， HONIG B， SHIN Y K， et al. Electrostatic potentials near the surface of DNA： Comparing theory and experiment[J]. J Phys Chem， 1995， 99： 7782-7786.

[18] SANNER M F， OLSON A J， SPEHNER J. Reduced surface： An efficient way to compute molecular surfaces[J]. Biopolymers， 1996， 38： 305-320.

[19] SITKOFF D， SHARP K A， HONIG B. Accurate calculation of hydration free-energies using macroscopic solvent models[J]. J Phys Chem， 1994， 98： 1978-1988.

[20] KOTTALAM J， CASE D A. Langevin modes of macromolecules-applications to crambin and DNA hexamers[J]. Biopolymers， 1990， 29： 1409-1421.

[21] SIMONSON T. Macromolecular electrostatics： Continuum models and their growing pains[J]. Curr Opin Struct Biol， 2001， 11（2）： 243-252.

[22] BASHFORD D， CASE D A. Generalized born models of macromolecular solvation effects[J]. Ann Rev Phys Chem， 2000， 51： 129-152.

[23] STILL W C， TEMPCZYK A， HAWLEY R C， et al. Semianalytical treatment of solvation for molecular mechanics and dynamics[J]. J Am Chem Soc， 1990， 112（16）： 6127-6129.

[24] WEISER J， SHENKIN P S， STILL W C. Approximate atomic surfaces from linear combinations of pairwise overlaps （LCPO）[J]. J Comput Chem， 1999， 20（2）： 217-230.

[25] MICHAEL F， John K， Charles L， et al. MMTSB tool set： Enhanced sampling and multiscale modeling methods for applications in structural biology[J]. J Mol Graph Model， 2004， 22： 377-395.

[26] DAURA X， GADEMANN K， JAUN B， et al. Peptide folding： When simulation meets experiment[J]. Angew Chem Int Ed， 1999， 38： 236-240.

[27] HU J P， GONG X Q， SU J G， et al. Study on the molecular mechanism of inhibiting HIV-1 integrase by EBR28 peptide via molecular modeling approach[J]. Biophy Chem， 2008， 132： 69-80.

[28] MAREN K， HOLGER S， WIBKE E D. Structural basis for HTLV-1 protease inhibition by the HIV-1 protease inhibitor indinavir[J]. J Med Chem， 2014， 57： 6266-6272.

[29] PERRYMAN A L， ZHANG Q， SOUTTER H H， et al. Fragment-based screen against HIV protease[J]. Chem Biol Drug Des， 2010， 75： 257-268.