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.