分子对接和分子动力学模拟方法研究芬太尼类化合物与阿片受体的相互作用

分子对接和分子动力学模拟方法研究芬太尼类化合物与阿片受体的相互作用

apollo营养元素是apollo酶o的主要亚型之一。现在，这三种亚状体的基因已被克隆。根据研究，选择性激动剂激活脊柱和脊柱受体，可能会产生一些生理效应，尤其是疼痛、付款效应、心率减慢、呼吸抑制、肠蠕动抑制、甚至瘫痪和中毒。经典的apo三种受体激动剂和椰子醇是它们的抗剂。apo受体属于g蛋白偶联体（gcl）。apsr具有相同的基本结构，即一个细胞外氨基端区域、七个半膜区域和一个细胞内羟基端端末年。芬太尼(fentanyl)学名为1-苯乙基-4-(N-丙酰苯胺)哌啶,为阿片类镇痛药,它的镇痛效力约为吗啡的100-180倍,哌替啶的550-1000倍,并具有毒性低、对循环影响小、时效短(15-30min)、容易控制、术后自主呼吸恢复快等优点,故近年来越来越受到人们的重视.芬太尼、舒芬太尼和阿芬太尼等芬太尼类配体药效增强、起效迅速、作用消失快,静脉滴注容易控制止痛剂量、安全可靠.关于芬太尼类激动剂与阿片μ受体选择性结合的分子作用机制的详细研究报道不多,目前对芬太尼活性构象的报道中发现,芬太尼结构中的N-苯乙基和N-苯基基团的取向对活性构象贡献较大.本文采用芬太尼及本课题组设计的小分子配体与阿片μ受体蛋白进行分子对接和分子动力学模拟,以研究配体与蛋白间的疏水作用、氢键作用等关系,探讨其分子作用机制,对辅助设计新型高效镇痛剂和其他特殊用途的化合物有一定的意义.1模型和计算方法1.1结构1:gpcr受体分子动力学模拟我们前期的研究结果与文献报道基本一致,故阿片μ受体蛋白的三维结构采用Zhang等用同源模建方法搭建的结构.该结构以具有高解析度的视紫红激酶的X射线晶体结构为模板(PDB代码为:1F88).一般的GPCR受体经过同源建模后需要进行分子动力学模拟优化,目前,GPCR受体的MD模拟方法有两种,一种是先构建蛋白与脂水膜体系,然后用MD研究整个蛋白-膜-水体系.另一种是构建蛋白-水体系,主要用来研究配体结合性.本文采用后一种方法.1.2受体的遗传算法在IBMIntelliStationPOWER工作站上,采用AutoDock4.0程序对12个芬太尼受体激动剂小分子与阿片μ受体进行分子对接,对接参数:受体大分子的格点盒子大小为8nm×8nm×8nm,格点间距为0.0375nm,盒子中心位于受体中心,运用Lamarckian遗传算法,将局部能量搜索与遗传算法相结合,以半经验势函数作为能量打分函数,对小分子构象和位置进行全局搜索,对每个配体进行128次独立的对接实验.遗传算法程序关键参数为ga_pop_size300ga_num_evals2500000,ga_run=100,最后依据最低对接能和成簇分析的情况,选取合理的受体配体结合模式作为初步的复合物结构.1.3结构的提取和分子对接操作将MD模拟后得到的受体再次与配体进行对接,复合物的结构取MD达到稳态后的平均构象然后从复合物结构中提取出受体结构和配体结构再以新得到的受体按1.2节中的方法进行分子对接操作.1.4激动剂-溶剂-水gc-ms的系统优化方法用GROMACS程序包分别在水溶液体系中对12个芬太尼受体激动剂和阿片μ受体蛋白受体复合物进行了分子动力学模拟研究.12个复合物结构为1.2节所述对接的结果.每个复合物的质心位于立方体盒子中心,选择溶质原子到盒子壁的距离为0.7nm,然后加入水分子,水分子选用SPC模型采用GROMACS全原子力场,用2fs积分步长,用最陡下降法进行了1000步的能量优化,然后在300K下进行了10000步的限制性动力学优化,最后采用2fs的步长,使用热浴耦合使系统温度保持在300K,用Maxwell分布产生.每隔100步记录一次轨迹在MD模拟中,用PME(Particle-MeshEwaldelectrostatics)方法处理静电作用,范德华相互作用截断半径取0.9nm,用SUNWAY计算机集群并行计算模拟1200ps,直到能量达到稳态为止.分别对蛋白环境和溶剂环境下激动剂与环境之间相互作用能进行系综平均,从而计算出受体-配体结合自由能(△G)和结合常数(Ki).1.5剂能量项采用了MM-PBSA(molecularmechanicsPoissonBoltzmannsurfacearea)方法计算复合物结合自由能,此方法采用分子力学和连续介质模型估算复合物的结合自由能,并采用修正后的RobertYang的Perl脚本处理构象,修正后的脚本排除了程序中读取错误.MM-PBSA方法的计算方程式:式(1)可以分为气相能量项和溶剂能量项,气相能量项包含内能项(EMM)和熵部分(TSMM)项;溶剂能量项Gsolv可以分为极性(Gpolar,solv)和非极性(Gnon-polar,solv)项,见式(2).用MD模拟后提取能量文件中的蛋白-蛋白电性和vanderWaals相互作用而得到EMM项.参考了文献中的方法考虑了在对接研究中不同复合物的熵相同而忽略了熵部分TSMM项.本文根据分子动力学结果,计算复合物加权协方差矩阵再计算其熵部分数值.极性和非极性项的计算则采用APBS软件包,其电性项的参数grid-spacing取0.01nm.使用GROMOS9643a1力场参数设置原子电荷和半径,探针半径0.14nm,复合物介电常数设为1,溶剂的介电常数设为80.非极性项Gnon-polar,solv采用溶剂可及表面(SSASA)方法计算:其中γ=2.2kJ·mol-1·nm-2,β=3.84kJ·mol-1.每个复合物采用21个结合构象,选取方法如下:在全部1200ps模拟中,选取最后200ps为平衡状态,即从1000ps开始(包含1000ps)取样,间隔为10ps,选取一个结构,至1200ps结束,共21个构象,采用这21个构象的极性和非极性项的平均值作为计算值.2结果与讨论2.1阿片受体三维结构的结构人类阿片μ受体蛋白有七个跨膜(TM)区段,属于GPCR超家族.该受体由400个氨基酸残基组成,TM2、TM3、TM7区和第2胞内环区(IL2)具有较高的同源性,而TM1、TM4、TM5区则同源性较低.第2胞内环和第3胞内环区(IL2和IL3)是蛋白结合区域,这些区域在三种阿片受体中的高度相似性表明有与蛋白相互作用的可能性.阿片μ受体与κ、δ两种阿片受体结构的最大差异部分在氨基端、羧基端及第2、3胞外环区(EL2和EL3),这些区域可能是阿片受体配体结合区,是不同配体选择性的结构基础.该推测已经通过直接点突变得到证实.阿片μ受体三维结构见图1.2.2md模拟结果通过对初步的分子对接结果分析,发现芬太尼配体分子间距离0.5nm内的参与配体相互作用比较重要的氨基酸有TM1区段上的Asp114、Ala117;TM3区段上的氨基酸Ile144、Ser145、Asp147、Tyr148、Asn150、Met151;TM5区段上的氨基酸His226、Ile230、Tyr234;TM6上Trp293、Asn296、His297、Tyr300、Arg303;TM7区段上的氨基酸Cys321、Ile322、Gly325、Tyr326等,所有参加作用的氨基酸编号见表1.12个芬太尼类衍生物(包括已知活性在内的4个芬太尼类衍生物)和阿片μ受体蛋白的对接结果有多种取向和构象,按照对接能量排序,我们发现低能构象多集中于配体质子化的N原子正电中心,并与TM3上的Asp147电负性中心发生相互作用(如图2所示),此外还有几个主要的作用残基,分别为TM1的Asp114,TM3的Asp147,和TM6的His297,其中Asp147与芬太尼类衍生物的哌啶季氨正电荷有电性作用,这些研究结果与现有文献的报道是一致的.通过MD模拟研究,我们还发现复合物蛋白结构与无配体受体结构有较大差异(图2).因此,我们将MD模拟后得到的受体结构再次与各自配体进行分子对接研究,AutoDock4.0再次对接后的计算结合能结果见表2.表2计算结果显示,结合能除羟甲芬太尼(ohmefentanyl)外,其它与实验活性排序一致.我们分析认为,羟甲芬太尼结构中额外的羟基通过水分子介导与Tyr148残基形成氢键(见图3),而在对接过程中这个溶剂分子引入的额外的氢键未给予考虑,从而羟甲芬太尼对接评分排序与实验结果是不一致的.采用MM-PBSA方法,计算出的12个芬太尼类衍生物的结合自由能,以及与它们对接计算的结果对比见表3.用MM-PBSA方法计算得到的结合自由能,不但能够对抑制剂的结合强弱进行正确的排序,而且计算得到的数值和实验数值能够较好地符合.其中以芬太尼的计算结果最为接近,与本文前述AutoDock4.0计算结果相比较,MM-PBSA方法计算结果更为接近实验值,但计算数据与实验值仍略有偏差,除羟甲芬太尼外,计算值整体稍微偏高.2.3载铅样品的md模拟2.3.1平均构象的变化空载受体蛋白在MD模拟过程中均方根偏差(RMSD)值是衡量体系是否稳定的重要依据,α碳(Cα)和骨架的RMSD值在初始的100ps内变化大,随后基本稳定.蛋白结构经过2000ps的优化达到稳态.本课题中以空载受体蛋白在1200ps平衡后的平均构象作为分子对接的基础构象,研究了其内能的变化,如图4所示.结果表明,体系的内能变化不大,在短期模拟后即达到稳定状态,这可能是因为采用的初始结构已经进行了能量优化的原因,上述模拟结果表明我们所得系统已经达到了稳定状态.MD模拟中的Cα空间位置均方根(RMS)涨落是反映分子内部运动特征及柔性的一个重要参数,RMS值越高,柔性越大,否则相反.从图5所示的空载受体RMS图可以看出,全部序列中有7个跨膜区段(分别以TM1-TM7标示)柔性相对胞内(标示为IL)比胞外区(标示为EL)低,其TM4、TM5跨膜区氨基酸在平衡态时仍有一定的柔性.受体蛋白的柔性可用Cα序列位置(δ)的RMS变化来衡量,即观察蛋白的δ在模拟后与模拟前结构的变化,用位置偏离的RMS均方根偏差来衡量,如果δ的RMS变化较大,说明这些区段的柔性较大.如图6所示.2.3.2不同配比的受体结构对活性构象的影响如图7所示,以芬太尼受体复合物为例的所有复合物分子的模拟,总能量和势能都很快达到稳态状态.在开始的500ps内,RMSD变化稍大,500ps后通过RMSD值变化可知,已基本达到平衡状态(图7(b)).这说明复合物体系已经稳定.分析复合物的RMS数据(图8)发现,受体柔性较大的片段为TM1、TM4、TM5跨膜区段,而对活性口袋具有主要贡献的TM2、TM3、TM6、TM7跨膜区段,结构柔性较低,比较稳定.将芬太尼复合物蛋白部分在1200ps内的Cα与空载配体受体结构Cα比较得知,复合物跨膜区段IL2、IL3(胞内区段)和TM4区段构象变化较大.分析复合物的RMS变化(图8、图9),发现结合配体后,受体TM1、TM3、TM6、TM7区段柔性降低(图8),而IL2、IL3(胞内区段)EL3(胞外区段)和TM4的部分区段柔性增加,这说明配体能够稳定受体TM1、TM3、TM6、TM7的构象,通过改变增加部分区段构象和柔性来激活受体,配体对受体结构的影响不仅仅在活性区,可能还通过一系列构象和柔性变化影响到整个受体的功能,从而使蛋白转变为活性构象而发挥作用.我们前期通过构效关系和受体分子药理学研究表明,受体蛋白活性中心的变化影响生理活性的大小,最可能仅对刚性受体而言是正确的.本研究表明,对大多数柔性受体来说,非活性中心构象的变化也将影响整个受体的活性构象及其与药物分子的相互作用,影响到药物分子的生物活性.不同配体的复合物结构RMS变化有差异.Carfentanil24结合受体后,受体的TM6、TM7区段RMS变化小,结构更加稳定,而IL2、IL3、EL3、TM4部分区段的RMS较高(如图10所示),实验值较好的化合物在稳定TM6、TM7区段跨膜蛋白结构的同时,还增加了其他蛋白区段的柔性,这些区段柔性和构象的改变可能与蛋白功能有关.与受体结合能较低的配体2、21、22与受体结合后,受体各个区段的RMS变化都呈下降趋势,柔性降低(图11).这种变化说明配体2、21、22与受体结合后,使得受体结构柔性整体降低,但是由于TM4、IL2、IL3柔性可能与活性构象有关,柔性整体降低不利于蛋白向活性构象转变从而降低了活性.2.3.3风险物质的添加突变数据和分子对接的结果都表明TM2,TM3,TM7跨膜区段都参与了同配体的结合,因此有必要重点分析构成活性口袋的螺旋区结构变化.如图12所示,分析结合配体后受体螺旋区的RMSD的变化,发现结合配体后受体TM4、TM5的RMSD波动较其他的几个螺旋显著.我们还分析了经历同样的MD模拟后,复合物和空载受体结构中氢键数量的变化,发现结合配体后,受体螺旋区的氢键除TM5区段外,其他区段的氢键数量稍有增加,特别是在TM4片段,比空载受体增加了两对氢键的相互作用,这与RMS变化分析中结合配体后引发TM4区段的变化相一致,而结合配体后整个体系氢键的相互作用增加,复合物与空载受体氢键数量比为63比55,见表4.在空载受体TM1和TM7、TM2和TM7、TM2和TM3、TM5和TM6与TM6和TM7之间都有氢键的形成.经过分析,这些螺旋之间形成氢键的现象在G蛋白偶联受体家族中是普遍存在的.在结合配体之后,TM1和TM7、TM2和TM7、TM2和TM3、TM5和TM6跨膜区段之间的氢键作用消失,TM1和TM2、TM2和TM4跨膜区段之间的氢键作用增加,TM2、TM3、TM7是结合口袋所在的区域,我们