【精品论文】HIV-1 蛋白酶异位抑制剂体系的长时间分子.doc

精品论文hiv-1 蛋白酶异位抑制剂体系的长时间分子动力学模拟王加磊，孟现美，张庆刚，张怿慈，张少龙5（山东师范大学物理与电子科学学院，山东 济南 250014）摘要：采用新开发的 ff12sb 力场在 nvidia cuda gpu 上对 hiv-1 蛋白酶的活性位抑制剂体系和异位抑制剂体系分别进行了 100 纳秒的长时间分子动力学模拟,并用 mm-pb/gbsa方法计算了活性位抑制剂 tl-3 与 hiv-1 蛋白酶的结合自由能。异位抑制剂体系中分子片段2-甲基环己醇结合在 exo 位,有利于抑制剂被束缚在活性位点附近。异位抑制剂体系中抑制10剂 tl-3 与蛋白酶的结合自由能为-85.78 kcal/mol,活性位抑制剂体系中为-79.45 kcal/mol。这 些结果有助于深入了解 hiv-1 pr 的动力学过程,为设计新型强效抑制剂提供了新见解。 关键词：计算化学;hiv-1 蛋白酶;异位抑制剂;mm-pb/gbsa;cuda;amber ff12sb 力场 中图分类号：o641.1215long-time scale molecular dynamics simulations of allosteric inhibitor system for hiv-1 proteasewang jialei, meng xianmei, zhang qinggang, zhang yici, zhang shaolong(college of physics and electronics, shandong normal university,jinan shandong 250014)20abstract: in this paper, we performed 100 nanosecond long time scale molecular dynamics simulation for the allosteric inhibitor system and active site inhibitor system of hiv-1 protease respectively,employed the recently developed ff12sb force field on nvidia cuda gpu,both systems bind free energy of inhibitor tl-3 with hiv-1 protease were calculated using mm-pb/gbsa method. fragment 2-methylcyclohexanol,binds in the exo site,stabilizes the inhibited conformation of25the target. the free energy of inhibitor tl-3 binding to protease is -85.78 kcal/mol in allosteric inhibitor system and -79.45 kcal/mol in active site inhibitor system. these results might help gain abetter insight into the dynamics of hiv-1 pr,which may provide important guidelines for the design ofnovel potent inhibitors.key words: computational chemistry;hiv-1 protease;allosteric inhibitor;mm-pb/gbsa;cuda;30amber ff12sb force field0引言蛋白酶(hiv-1 protease)是艾滋病病毒生命过程必需的一种酶,蛋白酶抑制剂药物作为高 效抗逆转录病毒疗法(haart,俗称鸡尾酒疗法)的重要组成部分已经在帮助艾滋病人存活方35面发挥了很大作用1。美国食品药物管理局(food and drug administration,fda)目前批准的九 种蛋白酶药物都是以活性位(active site)为靶标进行抑制的。但是,艾滋病毒已进化出对现有药 物具有耐药性的变异形式,美国多于 10%的新艾滋病病例是由病毒的耐药株引起的,因此设计 新机制的蛋白酶药物迫在眉睫2。a. perryman 等人近来发现了 hiv-1 蛋白酶上两个在活性基金项目：高等学校博士学科点专项科研基金(20093704110001);国家自然科学基金(11274206);山东省自然 科学基金(zr2011hm048)作者简介：王加磊(1986-),男,硕士研究生,主要从事分子动力学模拟研究通信联系人：张少龙(1963-),男,副教授,分子动力学. e-mail: - 10 -位之外的新结合位点(allosteric binding sites,异位结合位点):exo 位和 outside/top 位(图 1),以及40结合在这两个异位结合位点上的化合物,2-甲基环己醇和 6-吲哚甲酸,它们既能结合到新位点 又能稳定活性位抑制剂的构象,虽然尚未表现出抑制作用,但可以成为开发针对耐药性艾滋病 病毒的新药的出发点3。这一发现奠定了开发一类新型抗艾滋病毒药物的基础,沿此路线设计 的新药可以加强现有疗法的疗效,治疗这种疾病的耐药类型,减慢病毒的耐药性进化。hiv-1 蛋白酶的结构如图 1 所示,它由各包含 99 个氨基酸残基的两条肽链组成,活性位由45asp25/25-thr26/26-gly27/27残基组成,asp25/25在催化中起关键性的作用4,挡板区(flap)位 于活性位上面并将活性部位遮盖。flap 对抑制剂或底物的进出起很大作用。图 1 hiv-1 蛋白酶的结构fig.1 the crystal structure of hiv-1 protease50分子动力学模拟可以得到体系中各原子运动随时间变化的细节,是理解蛋白质-配体相 互作用的有力工具。势函数的不准确(力场问题)和模拟时间的有限是分子动力学模拟成功应 用受到的主要限制。amber ff12sb 力场是目前最新的改进的非极化力场,它在 ff99sb 基础上加入新的主链和55侧链扭转势,主链骨架扭转势改进后更精确,侧链的改进是通过拟合更精确的量子从头算数据 获得的,与实验数据符合得比较好5。gpu 计算是一种大规模并行计算的新方法,该方法通过芯片内数百个甚至上千个处理核 心的通信和协作,以最高可超过传统方法百倍的速度解决复杂的计算问题。本研究应用最新的 amber ff12sb 力场在 nvidia gtx670 gpu 上对 hiv-1 蛋白酶异位60抑制剂体系进行了长时间分子动力学模拟计算,得到了比以前更精确的结果。1理论与计算方法结合自由能是蛋白质与抑制剂结合的重要指标,结合自由能的正负反映了蛋白质与抑制 剂结合的可能性,结合自由能的大小反映了结合紧密程度6。mm-pb/gbsa(molecular mechanics poisson-boltzmann/generalized born surface area)是计算溶液环境中结合自由能的65流行方法。受体和配体的结合自由能为：gbinding= gcomplex gligand greceptormm-pb/gbsa 方法把各类的自由能(complex,ligand,receptor)分解为气相能量焓,溶解自 由能和熵项7：g = egas + esolv t s= ebat + evdw + ecoul + gsolv, polar + gsolv,nonpolar t s70其中 ebat 为成键的键、角和扭转角项能的和,evaw 和 ecoul 为非成键的范德华能和静电相 互作用能,gsolv,polar 为溶解自由能的极化项,gsolv,nonpolar 为溶解自由能的非极化项8。ebat、evaw、 ecoul 之和为完整的气相力场能,即 mm 部分。溶解自由能极化项 gsolv,polar 可通过 poisson-boltzmann(pb)方程求解,也可用广义波恩(gb)模型求解,溶解自由能非极化项 gsolv,nonpolar 通常用经验公式得到：75gsolv,nonpolar= sasa + const其中,sasa 是溶剂可接触表面积(solvent-accessible surface)。另外,溶解自由能非极化项gsolv,nonpolar 在更精确方法下可以分解为排斥项和吸引项9。2计算过程分子动力学模拟使用 amber12 和 ambertools12 软件包。初始构象取自蛋白质库(pdb id:803kf0),pdb 中包含 hiv-1 pr、3tl、4dx、dmso 和 bme, 其中二甲基亚砜(dimethyl sulfoxide,dmso)和-巯基乙醇(beta-mercaptoethano,bme)在 hiv-1 蛋白酶的表达、提纯和 结晶中作为有机溶剂起溶解和稀释作用3, 在本研究中予以忽略; 提取 pdb 中 hiv-1pr+3tl+4dx 作为异位抑制剂体系,提取 hiv-1 pr+3tl 作为活性位抑制剂体系;由于结晶水 在 hiv-1 蛋白酶与抑制剂的结合中起了非常重要的作用10，pdb 中的全部结晶水予以保留。85将得到的结构在 ambertools12 的 leap 模块中补全缺失的氢原子并添加 ff12sb 力场参数。 抑制剂 3tl 和分子片段 4dx 的部分电荷用 sqm 模块中的 am1-bcc 方法计算。两个体系周 围加有 9 个氯离子以保持体系中性。使用 shake 算法进行键约束计算11。将两个体系质心 置于截角八面体盒子中心,盒子壁到溶质原子的距离为 10 ,加入水分子,水分子采用 tip3p 模型。接着采用 amber12 的 pmemd 模块对两个体系进行能量最小化12,能量最小化先是采90用 25 000 步最陡下降法,继而是 25 000 步共轭梯度法,然后用 500 000 步将两个体系温度从 0k 升到 300 k,并在此温度下进行 2 000 000 步动力学平衡。最后,对两个体系分别进行 100 ns 的分子动力学模拟,运动方程积分步长(时间步长)2 fs,每隔 1 ps 记录一次坐标文件,用于后续 的分析。953结果和分析3.1动力学轨迹分析100105图 2 分子片段 4dx 和蛋白酶表面 exo 位结合示意图fig.2 space filling rendering of the exo site on the pr surface into which the fragment 4dx bindx 射线晶体衍射实验发现分子片段 2-甲基环己醇(4dx）结合在蛋白酶的非活性位点 exo 位,exo 是由一对环区(残基 14/1418/18和 38/3842/42)和一链(残基 59/5965/65)形成 如图 2 所示的蛋白酶表面凹槽3。2-甲基环己醇的羟基与残基 gly17 形成氢键,它的甲基环己 基(the methylcyclohexyl group)与 lys14 和 leu63 存在疏水性结合。exo 位的 38/3842/42 和 59/5965/65序列紧邻蛋白酶的挡板区(flap),flap 的变化活动与 exo 位为反相关关系,当 flap 打开时,exo 位凹槽会被压缩,反之 flap 关闭时 exo 位凹槽就会扩大,类似于杠杆系统3。图 3 exo 位的 lys14、leu63、gly17 残基示意图fig.3 the graph of lys14、leu63、gly17 residues in exo site110115图 4 异位抑制剂体系挡板区的运动fig.4 the graph of the motion of flap in inhibitor system图 5 分子动力学过程中骨架碳原子相对初始结构的均方根偏差随时间变化fig.5 rmsd vs. time graph of the backbone c atoms relative to their initial structures in md120动力学运动轨迹由不同时刻体系中所有原子的坐标组成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 位凹槽开口更大,同时比较两体系得异位抑制剂体系中挡板区活动方125130135向如图 4 所示,即挡板区关闭更牢固,这说明 flap 活动与 exo 位为反相关关系。系 统是否 达到动 态平衡 通常用 体系的 温度密 度总 能量和 均方根 偏差 (root-mean-square deviation,rmsd)来衡量14。两个体系相对初始结构的 rmsd 如图 5 所示, 可以看出在 100ns 分子动力学模拟过程中两个体系之间蛋白质骨架碳原子的均方根波动差 异明显,异位抑制剂体系的 rmsd 从开始到 10ns 缓慢上升至 1.0,并在随后的 90ns 内保持在1.0 上下浮动。活性位抑制剂体系的 rmsd 在开始时快速跃升到 1.2,并在后面模拟时间段 内保持在 1.2 上下浮动。两个体系在 30ns 以后达到平衡。整体上观察,异位抑制剂体系的 rmsd 比活性位抑制剂体系要小,这说明异位抑制剂体系更稳定。图 6 分子动力学过程中挡板区骨架碳原子相对初始结构的均方根偏差随时间变化fig.6 rmsd vs. time graph of the flap backbone c atoms relative to their initial structures in md图 6 为两个体系挡板区骨架碳原子相对于初始结构的 rmsd 随时间变化示意图。可以 看出异位抑制剂体系中的挡板区在整个模拟时间段无明显起伏,保持在 0.8 上下浮动,较为 稳定。活性位抑制剂体系中,挡板区在整个模拟时间段内起伏变化剧烈,比较不稳定。140145图 7 两体系抑制剂相对初始结构的 rmsd 随时间变化fig.7 rmsd vs. time graph of the inhibitor atoms relative to their initial structures图 7 为两个体系抑制剂 tl-3 相对于初始结构的 rmsd 随时间变化示意图。可以看出在 开始的 30ns 内,由于两体系未达到平衡状态,rmsd 变化剧烈.。达到平衡后异位抑制剂体系 的 rmsd 值整体低于活性位抑制剂体系,因此 4dx 与 exo 位的结合有利于抑制剂在活性位 点与蛋白质结合,异位抑制剂体系稳定性好。3.2蛋白酶残基 b 因子分析150155160165图 8 两体系残基的 b 因子fig.8 b-factor graph of individual residue of the two systemsb 因子(b-factor)是衡量蛋白质某个残基在分子动力学过程中偏离其平均结构的指标,反 映了蛋白质残基的灵活性和柔性度。b 因子越大,说明相应残基灵活性大柔性度高15。我们重点对 exo 位和挡板区氨基酸残基 b 因子进行分析,挡板区由序列 43/4358/58 组成的肽段构成,exo 位由序列 14/1418/18、38/3842/42和 59/5965/65构成。他们在 图 8 表示为：exo-1 位(1418、3842、5965); exo-2 位(113117、137141、158164); flap-1(4258);flap-2(141157)。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 位影响比较明显,这些残基活动范围明 显减少,刚性更强,更稳定,使活性位抑制剂不易脱落。1703.3结合自由能分析本研究用 mm-pb/gbsa 方法分别计算了两个体系中活性位抑制剂 tl-3 与 hiv-1 蛋白 酶的结合自由能,结果如表 1、2 所示。表中 gele 表示分子力场中静电能,gvdw 表示分子 力场中范德华能,ggas 表示气相能,gpbcal 和 gpbsur 表示 pb 计算方法中极性溶解能和非 极性溶解能,gpbsol 表示 pb 计算方法中极性和非极性的溶解能总和,ggb 和 ggbsur 表示175180185gb 计算方法中极性的溶解能和非极性的溶解能,ggbsol 表示 gb 计算方法中极性和非极性溶解能总和,gpbtot 和 ggbtot 分别表示 pb、gb 模型中最终的结合自由能12。其中,气相 能为分子力场中静电能、范德华能和内能三者之和。由于分子力场中内能在单轨迹近似方法 中总是为零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)tab.1 binding free energy from mm-pbsa(kcal/mol)gelegvdwggasgpbcalgpbsurgpbsolgpbtot异位抑制剂体系-99.67-94.64-194.31121.06-12.54108.53-85.78活性位抑制剂体系-83.20-89.27-172.47103.06-10.0493.02-79.45差值-16.47-5.37-21.8418.00-2.5015.51-6.33表 2 mm-gbsa 计算结合自由能(kcal/mol)tab.2 binding free energy from mm-gbsa(kcal/mol)gelegvdwggasggbggbsurggbsolggbtot异位抑制剂体系-99.67-94.64-194.31100.95-12.5488.42-105.89活性位抑制剂体系-83.20-89.27-172.4788.86-8.4180.44-92.03差值-16.47-5.37-21.8412.09-4.137.98-13.861901954结语本研究对 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 位凹槽,扩大了凹槽体积,使挡200板区关闭更牢固,同时,两体系挡板区 b 因子的对比表明异位抑制剂体系挡板区活动范围减少,刚性变强,更稳定。这些分子动力学结果证实异位抑制剂体系使活性位抑制剂 tl-3 与 hiv-1蛋白酶的结合更牢固,因此体系更稳定。这为开发艾滋病新型抑制剂提供了有价值的线索。参考文献 (references)2052102152202252301 jasklski 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 inhibitorj. biochemistry,1991,30(6): 1600-1609.2 watkins t,resch w,irlbeck d,swanstrom r. selection of high-level resistance to human immunodeficiencyvirus type 1 protease inhibitorsj. antimicrob agents chemother, 2003,47(2): 759-769.3 perryman al,zhang q,soutter hh,rosenfeld r,mcree de,olson aj,elder je,stout cd. fragment-based screen against hiv proteasej. chem biol drug des,2010 march;75(3):257-268.4 时术华,扈国栋,陈建中,等.分子动力学模拟研究质子化态在 hiv-1 protease-indinavir 复合物中的作用j.化学学报,2009,67(24):2791-2797.5 lindorff-larsen k,piana s,palmo k,maragakis p,klepeis jl,dror ro,shaw de.improved side-chain torsion potentials for the amber ff99sb protein force fieldj.2010;78:1950-1958.6 侯廷军, 徐筱杰. 基于分子动力学模拟和连续介质模型的自由能计算方法j. 化学进展,2004,16(2):153-158.7 bill r.miller,iii,t. dwight mcgee,jr.,jason m. swails,nadine homeyer,holger gohlke,adrian e. roitberg.mmpbsa.py:an efficient pro