β-乳球蛋白与EGCG结合规律的分子动力学探究

孔庆新， 黄业传， 徐伟平， 张径舟， 刘碧林， 李思阳

（1.重庆化工职业学院 制药工程学院，重庆 401228；2.荆楚理工学院 生物工程学院，湖北 荆门 448000；3.西南政法大学，重庆 401120；4.江苏食品药品职业技术学院 药学院，江苏 淮安 223001）

牛乳是一种营养丰富、老少皆宜的食物。近年来，以牛乳为主要原料，添加茶或水果开发出了奶茶、果奶等乳饮料产品，产品形式多样、风味更佳、营养全面。然而，茶、水果等植物原料中含有的多酚类物质易与牛乳中蛋白质作用形成复合物，影响多酚的活性和牛乳蛋白质的性质，进而影响产品的功能和营养特性。很多学者对蛋白质与多酚的结合进行了较深入的研究，如两者间的结合类型、结合位点、结合机理等[1-5]。小分子物质与蛋白质的交联主要靠非共价键[6]，如氢键、疏水相互作用、范德华力、静电作用，在一定条件下也可共价交联[7-9]。

外源植物成分的加入，导致牛乳pH发生变化，进而影响多酚与蛋白质的结合，从而影响产品品质。很多学者开展了pH对多酚与牛乳蛋白质结合影响的研究，如川陈皮素在酸性（pH 3.0）与中性（pH 7.4）条件下[10]、阿魏酸在酸性（pH 2.4）和中性条件（pH 7.3）下[11]、芹黄素在pH分别为2.6、6.2、7.1、8.2条件下[12]与β-乳球蛋白的交联，没食子酸和表没食子儿茶素没食子酸酯在pH分别为3.0和7.0条件下与天然乳清蛋白的交联[13]。刘婵[14]模拟了胃（pH 2.0）、肠（pH 7.5）环境下，茶多酚与乳蛋白、消化酶间的结合及三者的竞争情况。

几十年来，学者们将光谱技术、波谱技术、量热技术等用于研究蛋白质的结构，但测得的只是静态结构，而生物体内蛋白质一直处于动态变化中。近年来，分子模拟已成为研究蛋白质结构的一种有效手段，成为继实验和理论手段后，从分子水平认识蛋白质的第3种手段，是一种具有足够小的时间尺度和空间尺度的强大模拟技术[15-16]。用于蛋白质和小分子物质结合的主要分子模拟手段有分子对接和分子模拟。一些学者采用分子模拟对牛乳蛋白质与小分子物质的结合进行了研究，如柑橘黄酮[17]、姜黄素[18]、茶多酚[14]、芦丁[19]等。关于不同pH条件下多酚与乳蛋白的交联，一些学者也进行了分子模拟[10-12，14]，刘婵[14]的研究涉及茶多酚与β-乳球蛋白的交联，但该研究只使用了分子对接对两者的结合位点和结合能进行探索，分子对接虽能研究蛋白质与小分子物质间的结合模式、结合能，以及蛋白质上参与结合的残基，但无法获得结合的动力学过程，而分子模拟能更深入地研究其结合机制。

β-乳球蛋白是牛乳中的一种重要蛋白质，因其具有良好的功能和营养价值，广泛应用于食品中。β-乳球蛋白对疏水及两性配体均具有良好的亲和能力[20]。EGCG（表没食子儿茶素没食子酸酯）是茶多酚中主要成分之一，因此，本研究中以β-乳球蛋白和EGCG为研究对象，同时用分子对接和分子模拟技术探讨pH为2（酸性）、5（接近等电点）和8（微碱性）条件下，两者结合能和结合机制的差异，以期为进一步提高茶乳饮料的质量奠定理论基础。

1 材料与方法

1.1 β-乳球蛋白与茶多酚EGCG的分子对接

β-乳球蛋白从RCSB网站下载（ID号：3npo），下载后用PyMOL软件手动除去结晶水。小分子物质EGCG用ChemDraw Professional 17.1和Chem3D 17.1制备，并用其自带的MM2力场进行结构优化，通过AutoDock Tools（1.5.6版本）对β-乳球蛋白和EGCG进行处理，其中EGCG加氢、计算gasteiger电荷，β-乳球蛋白加氢、计算gasteiger电荷、合并非极性氢，并分别生成β-乳球蛋白和EGCG的pdbqt文件。对接采用全盲对接，EGCG设置为全柔性，盒子大小设定为40×40×40Å，中心坐标为默认值，设置间隔为1Å。对接时采用AutoDock Vina搜索前20个得分最高的构象，对接参数energy_range设置为5，exhaustiveness设置为100。

1.2 不同pH条件下EGCG与β-乳球蛋白结合的分子模拟

从上述1.1得到的能量最高的构象中，剥离出EGCG的pdb结构，并在ATB网站生成其top文件。采用PROPKA在线计算平台计算β-乳球蛋白在pH分别为2、5、8时的氨基酸残基质子化状态，根据计算结果来设定β-乳球蛋白的不同质子化状态以模拟相应的pH。采用Gromacs（2019.03）软件进行分子模拟，先将EGCG合并到β-乳球蛋白中，修改复合物的top文件；选用GROMOS54a7力场，将β-乳球蛋白置于立方体水盒子中，采用SPC水模型，将β-乳球蛋白离盒子边缘最短距离设置为1 nm并添加Na+或Cl-使体系达到电中性，使用最速下降法对体系进行能量最小化。在NVT和NPT系综下分别进行400 ps的平衡，达到预设的温度和压力（温度为300 K，压力为1 MPa）。最后进行150 ns的分子模拟，使用LINCS算法约束所有键，使用PME计算静电作用，范德华力使用截断半径为1.4 nm。选用Parrinello-Rahman压浴方式，Isotropic控压方式，V-rescale控温方式。模拟步长为2 fs，每10 ps贮存1次数据，每个处理平行模拟2次。

模拟结束后，使用gmx的rms命令分别分析150 ns模拟过程中各pH条件下EGCG和β-乳球蛋白的均方根误差（RMSD）；并用gmx的相应命令分析各pH条件下β-乳球蛋白结构稳定后（100~150 ns）蛋白质残基的波动（RMSF）、蛋白质二级结构、蛋白溶剂可及表面积。分析时每100 ps选取一个数据点，作图用Origin 8.0软件。

1.3 MM-PBSA计算EGCG与β-乳球蛋白的结合能

提取模拟过程中100~150 ns的轨迹，借助g_mmpbsa实 用工具，利用MM-PBSA（Molecular mechanics poisson–boltzmann surface area）方法计算EGCG与β-乳球蛋白之间的结合能[21-22]。

1.4 EGCG与β-乳球蛋白结合的二维平面图和表面结构图

提取各pH条件下EGCG与β-乳球蛋白复合物在100~150 ns的平均pdb，用LigPlot+（V 2.2）分析EGCG与β-乳球蛋白结合位点的疏水相互作用和形成的氢键，并在PyMOL软件中绘制β-乳球蛋白分子的表面结构图及与EGCG结合情况。

2 结果与分析

2.1 β-乳球蛋白对酚类的结合位点

如图1所示，β-乳球蛋白单体由162个氨基酸组成，由8条反向平行的β-折叠构成桶状结构，其外侧含有1个α-螺旋。

图1 EGCG在β-乳球蛋白的主要结合位点Fig.1 Main binding sites of EGCG to β-lactoglobulin

2.2 EGCG与β-乳球蛋白结合的分子动力学

通过AutoDock Vina对接，发现在前20个能量较高的结合姿势中，EGCG在β-乳球蛋白的结合区域主要有3个，如图1所示，其中最主要是1号区域，位于疏水腔的上部，其结合能也较高；另外两个是位于β-乳球蛋白表面的2号和3号位点，其中3号位点只有20个对接结果中的1个，且结合能较低。关于多酚在β-乳球蛋白结合位点的报道较多，其中但倩誉[10]的研究认为在pH为7.4时，小分子物质（川陈皮素）主要结合在疏水腔，而在酸性条件下（pH 3.0）时小分子物质结合在β-乳球蛋白表面；Abdulatif等[19]也认为芦丁既可以结合在疏水腔也可以结合在蛋白质表面，其中疏水腔的结合力更大。尽管其研究中小分子物质都是偏疏水性，本研究中茶多酚偏亲水性，但研究结果与本研究结果基本一致，因为本研究中下载的β-乳球蛋白为中性pH。而Gholami等[3]研究发现柚皮苷只吸附在β-乳球蛋白表面；Zhu等[12]研究阿魏酸与β-乳球蛋白的结合时，发现在酸性条件下（pH 2.4）阿魏酸主要结合在疏水腔，而在中性条件下（pH 7.3）阿魏酸则主要吸附在表面蛋白二聚体的结合处；Li等[23]研究4种不同黄酮与β-乳球蛋白的交联时，发现均在2号位点的结合能最大。以上结果的差异可能主要与小分子物质结构、初始蛋白质结构或模拟条件的差异有关。

2.2.1 RMSD（均方根误差）在模拟过程中的变化RMSD（均方根误差）是衡量特定时间蛋白质结构与原始构象的平均偏差，是评价研究体系是否稳定的重要指标，图2是不同质子化状态下体系中EGCG和β-乳球蛋白的RMSD对比。可以看出，EGCG在前60 ns波动较大，在60 ns后基本保持稳定，但pH为2时仍有一定幅度的波动，且就整体波动而言，pH为2时EGCG波动大于pH为5和8时，说明在酸性条件下EGCG与β-乳球蛋白结合过程中波动更大。就β-乳球蛋白的波动而言，3个pH条件下RMSD相当，且均在很短的时间内就达到了平衡（20 ns）。因此，总的来说，不同pH条件下EGCG与β-乳球蛋白的结合很平稳，EGCG只在有限的范围内存在一定的波动，不会对蛋白质结构造成大的影响。Sahihi等[17]模拟柑橘黄酮与β-乳球蛋白的交联时也发现结合很平稳，RMSD较小；Gholami等[3]研究柚皮苷和β-乳球蛋白的结合时发现小分子物质波动比β-乳球蛋白小；Abdulatif等[19]研究芦丁与β-乳球蛋白结合时，发现小分子物质与β-乳球蛋白的波动相当，这些都与本研究的结论基本一致。本研究中，pH为2时β-乳球蛋白与EGCG的波动相当，而pH为5或8时，EGCG的波动略小于β-乳球蛋白。

图2 不同pH条件下模拟过程中EGCG和β-乳球蛋白的RMSD对比Fig.2 RMSD comparison of small molecule and protein during different pH simulation

2.2.2 RMSF（均方根涨落）在模拟过程中的变化RMSF主要反映蛋白质中各氨基酸残基在不同pH条件下的波动情况，波动越大，说明该残基柔性越大，平均柔性越高蛋白质结构越不稳定。从图3中可以看出，各体系RMSF相当，没有明显差异，说明不同pH总体上对残基波动影响不大。就EGCG和β-乳球蛋白结合口袋附近而言，残基波动较大的有86号 位的Ala、87号 位 的Leu和88号 位 的Asn。RMSF的结果与Sahihi等[17]模拟柑橘黄酮与β-乳球蛋白的结合结果基本一致。

图3 分子模拟过程中不同pH条件下β-乳球蛋白的RMSF对比Fig.3 RMSF comparison of different pH during molecular dynamics simulation

2.2.3 不同pH对蛋白溶剂可及表面积的影响从表1可以看出，总溶剂可及表面积、疏水表面积及亲水表面积均是随pH的增加先减少再增加，亲水表面积大于疏水表面积，因此β-乳球蛋白在不同pH条件下总体上均为水溶性。pH从2增加到5时，总溶剂可及表面积的减少可能是由于β-乳球蛋白从单体变为二聚体的原因，而pH从5增加到8时总溶剂可及表面积的增加可能是由于Tanford的转变，即β-乳球蛋白上的EF环（残基85~90）形成了疏水腔的盖子，其在pH低于6.2时关闭，而在高于7.1时打开，盖子关闭时，部分溶剂可及表面积被掩埋[12]。从2.2.2的分析也可以看出，波动较大的86~88号残基正是位于EF环上，因此可以推测Tanford转变确有发生。

表1 不同pH对蛋白溶剂可及表面积的影响Table 1 Effect of different pH on the solvent accessible surface area of lactoglobulin

2.2.4 不同pH对β-乳球蛋白二级结构的影响从表2可以看出，随着pH的升高，α-螺旋和β-折叠均有所降低，pH为5时降至最低；伴随着α-螺旋和β-折叠的降低，无规则卷曲结构有所增加；而转角和弯曲的变化是相反的，pH为2时β-乳球蛋白中转角显著高于另两个pH，而弯曲的变化则相反。二级结构的变化，意味着蛋白质与多酚在结合时结合位点和结合能可能也会发生变化，从而导致结合的稳定性有所不同。Abdollahi等[11]研究阿魏酸与β-乳球蛋白结合时，也发现pH升高时，α-螺旋和β-折叠的比例均降低，而无规则卷曲增加，这与本研究中的结果一致。而Zhu等[12]则发现在不同pH条件下β-乳球蛋白的二级结构变化不大。

表2 不同pH对β-乳球蛋白二级结构的影响Table 2 Effect of different pH on the secondary structure of β-lactoglobulin单位：个

2.3 MM-PBSA分析EGCG与β-乳球蛋白的结合能

从表3可以看出，pH为5时β-乳球蛋白与EGCG的结合能最大，其次是pH为2时，而pH为8时两者的结合能最低。从具体的能量组成来看，在3个pH条件下非极性溶剂化能相当，均较低；范德华力为β-乳球蛋白与EGCG结合的主要作用力，其在pH为5时最高，pH为2时最低；静电作用力在pH为5时最低，pH为8时最高，pH不同会使部分残基的侧链带有不同电荷而引起静电作用力的差异，pH为5时最接近β-乳球蛋白的等电点，所带净电荷最少，从而静电作用力也最低；极性溶剂化能阻碍β-乳球蛋白与EGCG的结合，且随pH的增加阻碍作用也增大。Zhan等[24]研究β-乳球蛋白与辣椒素结合时，用MM-PBSA计算得到的结合能是-142.744 kJ/mol，范德华力为主要结合力，与本研究结果基本一致。Abdulatif等[19]研究芦丁和β-乳球蛋白的结合时，用MM-GBSA计算出两者在蛋白质表面的结合能为-285.647 kJ/mol，大于本研究的结果，这可能与小分子物质结构不同有关，也可能与计算方法不同有关。Zhu等[12]研究了β-乳球蛋白和芹黄素的结合，发现pH为7.1和pH为8.2时芹黄素主要靠疏水相互作用吸附在疏水腔，而pH为2.6和pH为6.2时主要是靠氢键和范德华力吸附在β-乳球蛋白表面，不同pH条件下吸附能（E）为E（pH 6.2）>E（pH 8.2）>E（pH 7.1）>E（pH 2.6），与本研究中的结论有所差异，因为该研究中先在各pH条件下分别进行对接，并用各自的最佳对接构象来进行分子模拟和结合能计算，而本研究中起始构象均在疏水腔，因此模拟的起始构象不一样。一些学者也研究了小分子风味物质，如醛、酮、酯等在不同pH条件下吸附于β-乳球蛋白的情况，发现在pH为3~9时，pH越高吸附作用越强[25-26]，原因可能是小分子风味物质相对分子质量较小，与β-乳球蛋白结合时空间位阻较小，可结合位点多，pH增加时β-乳球蛋白的结构变化更有利于这些小分子风味物质的结合或增加了结合位点。

表3 不同pH条件下β-乳球蛋白与EGCG结合的MM-PBSA结合能变化Table 3 Effect of different pH on the MM-PBSA binding affinity between protein and EGCG

从不同残基的能量贡献值来看，pH为2时排名 前10的 分 别 是Leu31、Pro38、Pro113、Asn109、Lys91、Ala37、Glu108、Ala118、Glu89、Leu1；pH为5时排 名 前10的 分 别 是Met107、Leu58、Ile84、Ile71、Val92、Val43、Ile56、Val41、Phe105、Leu31；pH为8时排 名 前10的 分 别 是Ile84、Ile71、Leu39、Val41、Glu62、Met107、Leu58、Val92、Pro38、Glu108。可以看出，pH为5时残基均为疏水氨基酸，这10个残基共贡献了-50.75 kJ/mol的结合能；pH为8和pH为2时分别有8个和6个疏水氨基酸，这8个和6个残基分别贡献了-40.22 kJ/mol和-24.42 kJ/mol的结合能，说明β-乳球蛋白与EGCG在pH为5时疏水相互作用最强，pH为2时最弱。另外，排名前10的残基中，pH为5和pH为8时有6个残基一样，而pH为2时与另外两个pH相同的残基很少，这说明酸性条件下β-乳球蛋白与EGCG的相互作用可能发生了较大的变化。

2.4 β-乳球蛋白与EGCG可能的结合机制分析

从图4可以看出，不同pH条件下，EGCG所处的疏水环境基本一致，均与周围的Leu39、Val41、Lys69、Ile71、Met107、Glu108存在疏水相互作用，氨基酸环境与Cheng等[27]报道的夭车菊素-3-O-葡萄糖苷在β-乳球蛋白的结合位点基本一致，也与Kanakis等[20]报道的EGCG与β-乳球蛋白结合后EGCG所处的氨基酸环境基本一致。从2.3的分析可以看出，虽然疏水环境相似，但不同pH时β-乳球蛋白疏水残基对结合能的贡献存在明显差异。不同pH条件下EGCG与β-乳球蛋白形成的氢键数目不同，pH为2时为5个氢键，分别是88、90、109的Asn和116号位的Ser，其中EGCG与109号位的Asn形成2个氢键；pH为5时形成了6个氢键，与pH为2时相比，Asn88与EGCG形成 了2个氢键；而pH为8时 只 形 成 了4个 氢 键，EGCG与Asn88和Ser116各形成1个氢键，与Asn109形成2个氢键，而Asn90则与EGCG不存在氢键。因此可以推断，不同pH条件下，EGCG与β-乳球蛋白结合力的差异也部分与氢键有关。从表3可以看出，pH为5时静电作用最弱，因此除了氢键以外，EGCG与β-乳球蛋白的其他弱静电作用也很重要[28]。

图4 不同pH条件下EGCG与β-乳球蛋白相互作用的二维平面图Fig.4 2D plot for interaction between EGCG and βlactoglobulin at pH 2，pH 5 and pH 8

Abdulatif等[19]研究发现芦丁与β-乳球蛋白的结合时形成7个氢键，但最佳结合位点与本研究有所差异。Cheng等[27]报道夭车菊素-3-O-葡萄糖苷也吸附在β-乳球蛋白的疏水腔，主要结合力为氢键和疏水相互作用，且形成的氢键有两个分别与Asn109和Ser116结合，与本研究中一致。贾晶晶[29]研究EGCG与β-乳球蛋白结合时发现结合力以疏水相互作用为主，同时也有范德华力和氢键。本研究中，EGCG与β-乳球蛋白之间的结合同时存在氢键、疏水相互作用、静电作用和范德华力。综合分析，pH为5时氢键、疏水相互作用、范德华力均为最大，虽然静电作用力最小，但其贡献有限，因此pH为5时EGCG与β-乳球蛋白的结合能最高。

为进一步探索高压对β-乳球蛋白表面结构及与EGCG结合的影响，用PyMOL软件绘制了不同pH质子化状态下100~150 ns模拟过程中β-乳球蛋白的表面结构及多酚与β-乳球蛋白的结合情况（见图5）。由图5可以看出，在不同pH条件下β-乳球蛋白的表面结构均有差异，特别是在位于疏水腔上部结合口袋处的结构差异较大，结构的差异也必然会导致β-乳球蛋白与EGCG的结合产生差异，且EGCG的结合姿势明显不同。在pH为5时EGCG更加伸展，也更加靠近疏水腔的内部。因此，pH为5时β-乳球蛋白与EGCG的结合能最高可能也与此有关，EGCG结构越伸展越能与更多的氨基酸侧链发生相互作用，增强氢键和疏水作用力。

图5 不同pH条件下蛋白质分子表面结构及与多酚结合情况Fig.5 Surface structure and polyphenol binding of protein at pH 2，pH 5 and pH 8

3 结语

茶多酚EGCG主要结合在β-乳球蛋白的3个位点，其中疏水腔的结合最稳定。选取疏水腔作为结合位点，进行EGCG与β-乳球蛋白结合的分子模拟，发现结合过程中EGCG在模拟前期的波动并不会影响β-乳球蛋白结构的稳定。pH为8时蛋白疏水表面积最大，pH为2时β-乳球蛋白的α-螺旋和β-折叠含量最高，pH为5时这三者均为最低。不同pH条件下β-乳球蛋白表面的结构差异较大，特别是位于疏水腔入口的结合位点处，这可能是影响EGCG与β-乳球蛋白结合能和结合姿势不同的重要原因之一。EGCG与β-乳球蛋白的结合能在pH为5时最高，其次是pH为2时，而pH为8时最低；pH为5时，尽管静电作用力低于另外两个pH，但由于其疏水作用力、氢键和范德华力最强，因此总的结合能较高。本研究中进一步明确了不同pH条件下蛋白质与多酚等小分子物质的结合机理，进而为 提高植物蛋白饮料的质量奠定了一定的理论基础。