木脂素类植物雌激素对ERα及ERβ结合亲和力的分子动力学研究

王 婷,王云飞,栾 锋,赵春燕

(1.烟台大学化学化工学院,山东 烟台 264005;2.兰州大学药学院,甘肃 兰州 730000)

植物雌激素(Phytoestrogen,PE)是一类存在于植物中,能结合并激活哺乳动物或人体内雌激素受体(Estrogen Receptor,ER)从而发挥雌激素样和抗雌激素双重作用的植物成分。它们分为四大类:木脂素类、异黄酮类、香豆素类以及二苯乙烯类[1]。其中,木脂素类主要是通过P-羟基苯乙烯单体氧化耦合产生的一类小分子量次生代谢产物,其主要以游离状态存在于如夹竹桃科、爵床科、马兜铃科等植物中[2],其广泛分布于植物的根、根状茎、茎、叶、花、果实、种子以及木质部和树脂等部位[3]。根据木脂素的来源[4]可以将其分为植物木脂素和动物木脂素。植物木脂素中含量最多的两种木脂素是开环异落叶松脂酚(secoisolarlciresinol,SECO)和罗汉松脂酚(matairesinol,MAT)。动物木脂素是植物木脂素在哺乳动物的肠道菌群[5]作用下产生的代谢产物,它们分布于动物的血清、血浆、尿液和粪便中,主要含肠二醇(enteroldiol,END)和肠内脂(enteroctone,ENL)。

不论是植物还是动物木脂素的结构都与雌激素(如雌二醇)的结构相似,它们对预防及治疗与雌激素有关的疾病如骨组织疏松症、心血管疾病、前列腺疾病以及冠心病等[6-9]有重要作用。据文献调研,针对木脂素类植物雌激素与雌激素受体之间起作用的结合亲和力研究未见报道,因此,本课题选择分子动力学模拟方法探究两种常见的植物木脂素(罗汉松脂酚和开环异落叶松脂酚)和两种动物木脂素(肠内脂和肠二醇)分别与雌激素受体之间的结合亲和力以及产生的相互作用,以期为研究植物雌激素的类激素作用提供帮助。配体结构式如图1所示。

1 实验部分

1.1 受体蛋白与配体的准备

雌激素受体ERα(PDB ID∶1GWR)和ERβ(PDB ID∶3OLS)的晶体结构是从RCSB Protein Data Bank[10]中获得的,并从该结构中得到ERα和ERβ的配体结合域。受体蛋白的准备是在Schrodinger 2015-2软件中的Protein Preparation Wizard模块[11]中,用OPLS-2005力场[12]对结构进行预处理所得。图1中四种木脂素类植物雌激素作为小分子配体,其初始结构从Pubchem数据库[13]中得到。配体的准备过程是在Schrodinger 2015-2软件的LigPrep模块中[14],对配体进行能量优化处理。无论是蛋白还是配体的初始结构都将被用于后续的分子对接过程中。

图1 配体结构式

1.2 分子对接

配体与蛋白的对接过程如下:首先,以配体为中心,与周围氨基酸残基形成一个立体空间;其次,将准备好的配体与受体蛋白进行对接,并对对接结果进行打分评估;最后,将对接产生的蛋白-配体复合物以PDB的格式输出,用于下一步的分子动力学模拟过程。该分子对接过程采用Schrodinger 2015-2软件中的Glide6.7对接程序[15-16]完成。

1.3 分子动力学模拟过程

分子动力学模拟过程采用软件Amber14[17]进行。受体蛋白和配体的分子参数分别采用FF03立场和GAFF立场[18-19]。复合物需要进行加氢处理,并且在体系周围加上1 nm的TIP3P型的立方体水盒子进行溶剂化;分别在ERα体系添加5个钠离子,ERβ体系添加2个钠离子,来中和负电荷使体系呈电中性。其中参数截断距离[20]设为1 nm,积分步程设为2 fs。

分子动力学模拟过程:首先对体系进行溶剂能量最小化、约束抗衡离子能量最小化以及无限制整体能量最小化,每次能量最小化都先进行3000步的最陡下降法优化,然后再进行2000步的共轭梯度法优化;其次运行升温动力学模拟,在50 ps内将体系温度从0 K升至300 K;最后在300 K及恒压的NPT系统中进行100 ns的平衡分子动力学模拟。

1.4 结合自由能的计算及分解

结合自由能的计算是使用Amber14软件从MD模拟后的复合物的动力学轨迹的最后10 ns中按固定间隔提取出1000个瞬时结构进行计算的。其计算公式如下:

ΔGbind=ΔGcomp-ΔGpro-ΔGlig=

ΔGMM+ΔGsol-TΔS,

ΔGMM=ΔGint+ΔGele+ΔGvdw,

ΔGsol=ΔGpol+ΔGnopol。

其中,ΔGbind是结合自由能,由复合物(ΔGcomp)、受体蛋白(ΔGpro)和配体(ΔGlig)的能量组成。ΔGMM是由内能(ΔGint)、静电相互作用能(ΔGele)和范德华能(ΔGvdw)组成的气相结合能,其中由于模拟中使用单轨迹方案,配体与受体结合的构象重组能部分被忽略,即ΔGint=0,ΔGsol是由极性(ΔGpol)和非极性(ΔGnopol)溶剂化能组成的溶剂化自由能。由于熵变的计算过程复杂,在本研究中忽略熵变对结合自由能的影响。

结合自由能分解是应用Amber14软件中的MM/GBSA自由能分解模块进行的[21]。该过程主要是计算每个氨基酸残基在配体与受体结合过程中对总结合自由能的贡献,并通过对自由能的分解来考察如范德华能、静电相互作用能以及溶剂化能等各部分能量对结合自由能的贡献大小。

2 结果和讨论

2.1 动力学轨迹分析

针对小分子配体与受体蛋白对接后形成的复合物进行了100 ns的分子动力学模拟,并计算了各体系的均方根偏差(Root Mean Square Deviation, RMSD),结果如图2。图2中所有体系在30 ns的模拟之后RMSD值基本稳定(ERα体系RMSD=0.11～0.22 nm;ERβ体系RMSD=0.14～0.20 nm)。所有体系的RMSD波动值均在0.30 nm以下,处于合理范围内,这说明所有体系经过模拟后均达到平衡,可以进一步分析模拟结果。

图2 ERα和ERβ体系中复合物在100 ns模拟过程中的RMSD值

针对体系中各个配体在模拟过程中的状态变化,分别计算了各体系中配体的RMSD值,结果如图3。图3中所有配体的RMSD值均在0.25 nm以下,处于平衡状态,说明配体在复合物形成过程中起到一定的稳定作用。

图3 ERα和ERβ体系中配体在100 ns模拟过程中的RMSD值

接着计算了复合物中每个氨基酸的均方根波动值(Root Mean Square Fluctuation, RMSF),以分析各个体系中受体蛋白的构象变化,如图4所示。图中ERα及ERβ两个体系中的RMSF曲线波动趋势基本相似,即氨基酸在结合过程中的波动相似,说明同一体系中的复合物的结合过程相似。ERα体系中的复合物在氨基酸PRO19、GLU41、LEU75、PRO94、ILE139、ASP172、HIE212和LEU229的周围均有相似的波动。这些氨基酸的RMSF值较小,说明其波动较小,对复合物的稳定起到一定作用。在ERβ体系中也有相同的体现,即复合物在氨基酸PRO15、LEU36、LEU77、PHE94、PHE115、ILE142、LEU158和HIE203的周围波动较小。

图4 ERα和ERβ体系中复合物在100 ns模拟过程中的RMSF值

2.2 氢键相互作用分析

针对对接后复合物的构象进行氢键相互作用分析,四种木脂素类植物雌激素与ERα及ERβ之间产生氢键作用如图5。氢键相互作用是稳定配体与受体结合的主要力量之一,就ERα与ERβ两个体系而言,ERβ体系中产生的氢键作用大于ERα体系(ERβ体系中氢键个数为2、3、3、2,ERα体系中为3、2、2、1),表明木脂素类植物雌激素对ERβ的结合亲和力强于对ERα的。此外,形成的氢键位于配体两端的结构更有利于复合物的稳定,因此,植物木脂素(罗汉松脂酚和开环异落叶松脂酚)对ERα与ERβ的结合亲和力都强于两种动物木脂素(肠内脂和肠二醇)。

图5 配体与受体蛋白之间的氢键相互作用

2.3 疏水相互作用分析

针对复合物在模拟过程中最后1 ns的轨迹构象进行疏水相互作用分析,结果如图6。图中显示结合口袋处多为疏水性氨基酸,故产生的相互作用为疏水相互作用,且根据图中结合空腔的左侧和底部的疏水势明显强于其他位置,说明疏水相互作用对复合物的稳定起重要作用。

图6 复合物的疏水相互作用(颜色从红色变为蓝色表示疏水势降低)

2.4 结合自由能分析

为深入地研究四种配体与雌激素受体之间亲和力,用Amber14中的MM/PBSA分别计算了四种木脂素类植物雌激素与ERα及ERβ之间经过100 ns的MD模拟之后所形成的两个体系的结合自由能,结果列于表1和2。从表中数据可知,结合自由能主要是由范德华能、静电能、溶剂化能等几种能量构成。其中,范德华能在几种相互作用能的贡献中占主要部分,且其和静电能、非极性溶剂化作用能的数值均为负值,说明其对配体-蛋白的结合起促进作用。相反的是极性溶剂化作用能对结合自由能贡献为正值,这说明其不利于配体与受体蛋白的结合。

表1 ERα体系中复合物的结合自由能

表2 ERβ体系中复合物的结合自由能

此外,根据表中数据,ERβ体系中四种复合物的结合自由能数值均小于ERα体系中复合物的自由能数值,结合自由能数值越低,形成的复合物越稳定,可见ERβ体系的复合物较ERα体系稳定。复合物ERα-罗汉松脂酚、ERβ-罗汉松脂酚、ERα-开环异落叶松脂酚和ERβ-开环异落叶松脂酚的结合自由能数值低于ERα-肠内脂、ERβ-肠内脂、ERα-肠二醇和ERβ-肠二醇的自由能。结合上述氢键相互作用的分析结果,说明罗汉松脂酚和开环异落叶松树脂酚对ERα和ERβ的结合亲和力强于肠内酯和肠二醇。

2.5 结合自由能的分解

应用Amber14中的自由能分解模块进行分解来分析各氨基酸对总自由能的贡献,结果如图7。图7中标注的氨基酸是对结合自由能贡献较大(超过±4.184 kJ/mol)的关键氨基酸。关键氨基酸位于配体与受体形成的结合空腔中,关键氨基酸的数目越多,配体与受体之间的相互作用越强,可据此判断关键氨基酸在结合过程中起到的稳定作用大小。复合物ERα-罗汉松脂酚、ERα-开环异落叶松脂酚、ERβ-罗汉松脂酚和ERβ-开环异落叶松脂酚的关键氨基酸的数目(依次为13、13、12、12)多于复合物ERα-肠内脂、ERα-肠二醇、ERβ-肠内脂和ERβ-肠二醇中关键氨基酸的数目(依次为9、8、11、10)。因此,前四种复合物比后四种复合物更稳定。

图7 氨基酸对结合自由能的贡献

对关键氨基酸进一步能量分解的结果见图8,在此过程中将自由能分解为极性部分(图中条纹部分)和非极性部分(图中黑色部分)。其中LEU38、ALA42、LEU79、MET80、LEU83、PHE96和LEU214等非极性氨基酸对结合自由能的贡献大于极性氨基酸如GLU45、ARG86等的贡献。并且非极性部分主要是范德华能贡献,能量为负值,进一步证实疏水相互作用是配体与受体之间结合产生的主要相互作用。然而极性部分中的极性溶剂化能的贡献值为正值,导致极性部分对总结合自由能的贡献减小,发挥的作用减弱。因此,非极性作用更有利于复合物的形成。

图8 关键氨基酸的能量分析

3 结 论

本文用分子动力学方法研究了四种木脂素类植物雌激素与ERα及ERβ的结合亲和力。氢键、疏水相互作用分析以及结合自由能的分析结果表明，四种香豆素类植物雌激素对ERβ的结合亲和力均高于对ERα的，且植物木脂素(罗汉松脂酚和开环异落叶松脂酚)对两种ER受体的结合亲和力均高于动物木脂素(肠内脂和肠二醇)对ER的结合亲和力。自由能分解结果说明非极性相互作用更有利于配体与受体的结合。