Interaction between β-lactoglobulin and EGCG under high-pressure by molecular dynamics simulation

The binding between β-lactoglobulin and epigallocatechin gallate (EGCG) under the pressure of 600 MPa was explored using molecular docking and molecular dynamics (MD) simulation. EGCG bound mainly in two regions with site 1 in internal cavity of the β-barrel and site 2 on the surface of protein. 150 ns MD was performed starting from the structure with the optimal binding energy at the two sites in molecular docking, respectively. It was found that the protein fluctuated greatly when small molecule bound to site 2 at 0.1 MPa, and the protein fluctuation and solvent accessible surface area became smaller under high-pressure. The binding of small molecules made the protein structure more stable with increasing of α-helix and β-sheet, while high-pressure destroyed α-helix of protein. The binding energy of small molecules at site 1was stronger than that at site 2 under 0.1 MPa, with stronger van der Waals and hydrophobic interaction at site 1 while more hydrogen bonds were present at site 2. The binding energy of both sites weakened under high-pressure, especially at site 1, causing the binding force to be weaker at site 1 than that at site 2 under high-pressure.


Introduction
Milk is an essential food for modern people. In addition to direct drinking, there are also many milk drinks, for example, tea was added into milk to make the product more diverse in form and taste and more comprehensive in nutrition [1]. Tea contains a lot of polyphenols, which have many physiological functions, such as antioxidation, reducing the incidence of cancer and cardiovascular [2]. Polyphenols are easy to form complex with milk protein, such as α-and β-caseins [3], β-lactoglobulin [4], thus affecting the quality of products and the biological activity of polyphenols. Many scholars have studied the binding of protein with polyphenols, such as the binding type, binding site and binding mechanism [5][6][7]. The binding between them is mainly non covalent, such as van der Waals force, hydrophobic interaction, electrostatic interaction and hydrogen bond, and covalent cross-linking also occurs under certain conditions [8,9].
At present, dairy products are still mainly sterilized by heat, thus having some negative effect on product quality. Ultra-high pressure is a kind of cold sterilization technology, which has good sterilization effect without destroying the flavor and nutrient of small molecules in food. Therefore, it is widely used in food, including dairy products [10,11]. The application of ultra-high pressure in milk will change the structure of macromolecular milk protein, and then affect the quality of dairy products. Some researchers reported the effect of high-pressure treatment on the structure of milk protein [12][13][14]. In terms of the effect of ultra-high pressure on the binding of polyphenols and proteins, Chen, Wang, Feng, Jiang, and Miao [15] reported that soybean protein and tea polyphenols were mainly crosslinked by hydrophobic interaction and hydrogen bond under 400 MPa, protein had protective effect on the biological activity of tea polyphenols, and the addition of tea polyphenols could also reduce the change of protein structure under high pressure. To our best understanding, there was no report on how high pressure affected the binding between polyphenols and proteins in milk drinks. Spectroscopy, calorimetry and atomic force microscopy are widely used to study protein structure, but only static structure is measured. Protein in organism is always in dynamic change. Now molecular dynamics simulation has become another means with enough small time and space scale to understand protein structure from molecular level after experimental and theoretical means [16]. The binding process and mechanism of milk protein with many small molecules were studied by molecular modelling, such as citrus flavonoids [17], curcumin [18], apigenin [19], capsaicin [20], and tea polyphenols [4]. Furthermore, a few scholars have simulated the structural changes of some proteins [21][22][23] or enzymes [16,24] under high pressure. However, the simulation of the interaction between small molecules and proteins under high pressure has not been reported so far.
β-lactoglobulin is rich in milk with a lot of genetic variants particularly A and B phenotypes [25], because of good function and nutritional characteristics, and good affinity for both hydrophobic and amphoteric ligands, it is widely used in food industry. Each monomer of βlactoglobulin contains 162 amino acids with a molecular weight of 18.4 kDa [26]. As shown in Fig 1, it has a calyx composed of 8-standard antiparallel β-sheets and an α-helix structure on the outside. EGCG, namely epigallocatechin gallate, is one of the most important components in tea polyphenols, which has very important health care function for human body [27,28]. Therefore, molecular dynamics simulation was used to study the changes of binding energy and binding mechanism between β-lactoglobulin and EGCG at 600 MPa in this study, so as to lay a theoretical foundation for further improving the quality of tea milk beverage and the application of high-pressure technology in milk beverage.

Preparation of ligands and protein
β-lactoglobulin was downloaded from RCSB website [29] with an ID of 3npo, which had no missing atom, no crystallized ligand and had a reasonable resolution. Then the crystal water of protein was removed with PyMOL software, and missing hydrogen atoms were added to protein by AutoDock Tools (version 1.5.6). EGCG was prepared by ChemDraw Professional 17.1 and Chem3D 17.1, and the structure was optimized by the self-contained MM2 force field. The coordinate files (PDBQT) of EGCG and β-lactoglobulin were prepared in AutoDock Tools environment respectively.

Molecular docking of β-lactoglobulin with EGCG
To investigate the binding sites of EGCG to β-lactoglobulin, a blind docking was performed using AutoDock Vina program [30] to search the top 20 optimal conformations. In the process of docking, the small molecule was set to be fully flexible and the center coordinate was the geometric center of the receptor, a box of 40 � 40 � 40 Å 3 was set with a grid spacing of 1 Å, and the docking parameters of energy_ range and exhaustiveness were set to 5 and 100, respectively.

Molecular dynamics simulation of EGCG binding to β-lactoglobulin under high pressure
Although molecular docking can explore the binding site and binding energy between small molecule and protein, they could not fully understand the dynamic changes in the binding process. Therefore, 150 ns of molecular dynamics simulation was conducted. The optimal conformations of site 1 and site 2 derived from docking simulations were selected as the initial structure to perform molecular dynamics simulation with GROMACS (2019.6) software package [31,32]. The topology parameters of protein and small molecule were generated by Gromacs program and Automated Topology Builder (ATB version 3.0) server [33], respectively. The GROMOS54a7 force field was applied to describe the interaction [34], and the complex was immersed into a periodic water box of cube shape using a single point charge (SPC) water model. The shortest distance between the protein and the edge of the box was 1 nm, and Na + was added to make the system electrically neutral. Energy minimization of the system was performed to relieve unfavorable interactions using the steepest descent method [35]. Simulation comprised of equilibration and production phases, firstly a 400 ps simulation in NVT ensemble was carried out to equilibrate the system, and temperature was maintained at 300 K using velocity-rescale method. Then another 400 ps equilibrium in NPT ensemble was followed using the Parrinello-Rahman barostat to maintain the pressure at 0.1 or 600 MPa. Finally, a 150 ns MD was performed with a time step of 2 fs. PME (Particle Mesh Ewald) was used to calculate long-range electrostatics, and the cutoff of van der Waals and Coulomb interaction was set to 1.0 nm [7,16,17,20]. The atomic coordinates were recorded to the trajectory file every 5 ps. At the end of the simulation, the periodic boundary conditions were removed, and then the corresponding commands of GROMACS analysis tools was used to evaluate the root mean square deviation (RMSD), root mean square fluctuation (RMSF), gyration radius (Rg), hydrogen bond between proteins, protein secondary structure and solvent accessible surface area. Origin 8.0 software was used for plotting.

Binding energy between small molecule and protein using MMPBSA
The 500 snapshot structures of each complex were extracted from the trajectory file of last 50 ns at an interval of 100 ps, and then MM-PBSA (molecular mechanics Poisson -Boltzmann surface area) method, which was implemented in g_mmpbsa tool, was used to calculate the binding energy between small molecule and protein [36]. Solvent accessible surface area (SASA) was used as a non-polar solvent model.

The 2-D plot and surface structural for interaction between EGCG and βlactoglobulin
The average PDB structure of the complex from 100ns to 150ns under different simulation conditions was extracted, and the hydrophobic and hydrogen bonding interactions between small molecule and protein at the binding sites were analyzed using LigPlot+(V 2.2) software [37]. The surface structure of the protein and its binding to small molecule was drawn by PyMOL.

The binding sites of EGCG to β-lactoglobulin
Through molecular docking, it was found that the top 20 preferential binding poses of EGCG in β-lactoglobulin were mainly distributed in two regions (The PDBQT file of EGCG, including the coordinates of 20 docking poses obtained by molecular docking was shown in S1 File), as shown in Fig 1, 14 of which were at the top of the hydrophobic cavity (region 1), and another 6 were located in the protein surface (region 2). The average docking energy in region 1 was greater than that in region 2. Under natural conditions β-lactoglobulin is in equilibrium between the monomeric and dimeric form, Abdollahi, Ince, Condict, Hung, and Kasapis found that the preferred binding site of ferulic acid lied at the interface of the two monomers when pH was at 7.3 [38], while the binding sites explored in this study were far from the crevice of dimer interface; therefore, the presence of dimer would not affect the binding of small molecules to protein in this study. There were many reports on the binding sites of polyphenols to β-lactoglobulin, most of which believed that small molecules can bind to the hydrophobic cavity and protein surface simultaneously, and the binding energy in hydrophobic cavity was stronger [39,40]. However, Gholami and Bordbar [7], Li et al [5], and Pu et al [41] reported that the small molecules were adsorbed only on the surface of β-lactoglobulin. These results may be related to the differences in the structure of small molecules, the initial structure of docking protein or the simulation conditions.

Molecular dynamics simulation of EGCG binding to β-lactoglobulin
Changes of RMSD in molecular dynamics simulation. RMSD (root mean square deviation) is a measure of the average deviation of protein structure from the original conformation at a given time and is an important indicator to evaluate the stability of the protein structure. Fig 2 showed the RMSD comparison of β-lactoglobulin under different treatment conditions. It could be seen that the RMSD value fluctuated greatly when EGCG bound to site 2 at 0.1 MPa, and reached equilibrium after 80 ns, while the other three systems achieved equilibrium at about 50 ns. The protein average RMSD after 100ns was 0.2898 nm when EGCG bound to site 2 at 0.1 MPa, which was significantly higher than that for site 1 at 0.1MPa (0.2459 nm), site 1 at 600 MPa (0.2464 nm) and site 2 at 600 MPa (0.2513 nm), indicating an instability of protein structure was observed when the small molecule was bound to site 2 under normal pressure, but this instability would disappear under high pressure. Vanaei, Parizi, Abdolhosseini, and Katouzian also reported an increased RMSD when they studied the binding of safranal to β-lactoglobulin [42].
Changes of RMSF in molecular dynamics simulation. RMSF (root mean square fluctuation) reflected the root mean square displacement of each amino acid residue in the protein compared with the average conformation. The greater the fluctuation of the residue, the greater the flexibility of the residue in the pressure treatment. The higher the average flexibility, and the more unstable the protein structure. It can be seen from Fig 3 that the protein residues fluctuated greatly when the small molecule bound to the protein at site 2 under 0.1 MPa, which was consistent with the RMSD analysis. In terms of the average fluctuation of all residues, the simulation under high pressure was smaller than that under normal pressure, indicating that high pressure could reduce the fluctuation of protein residues. These results were consistent with some other reports, among which Kurpiewska et al [21] found that the residue displacement of insulin decreased when it was treated at 200 or 500 MPa, and Hata, Nishiyama, and Kitao [43] also believed that the structural fluctuation of protein decreased under high pressure. In terms of specific amino acid residues, the ones at the entrance of hydrophobic cavity fluctuated greatly (Res110-115, Res85-87), while those near binding site 2 fluctuated slightly in the process of binding.
Changes of protein gyration radius in molecular dynamics simulation. Fig 4 showed the change of protein gyration radius (Rg) under different simulation conditions. The smaller the gyration radius was, the more compact the protein structure was. It could be seen that the Rg decreased gradually with the extension of simulation time and remained stable after 70ns.

PLOS ONE
No matter where the small molecule bound, the Rg of protein decreased at 600 MPa, showing that high pressure was conducive to the formation of compact structure of protein, which may be due to the deformation of hydrophobic cavity caused by water molecules entering into it under high pressure, resulting in the change of secondary and tertiary structure of protein, and then the decrease of protein volume [43].

PLOS ONE
Changes of hydrogen bonds between proteins in molecular dynamics simulation. It could be seen from Fig 5 that after the system was stabilized (after 100 ns), the binding of small molecules to protein at site 2 resulted in the decrease of intramolecular hydrogen bonds within the protein, which may be related to the formation of more hydrogen bonds between protein and small molecule. Hydrogen bond is one of the main non-covalent forces to stabilize protein structure, the drop of hydrogen bonds may result in the instability of protein. So, the increasing fluctuation of protein RMSD and RMSF when the small molecule bound to protein at site 2 may be related to the decrease of hydrogen bonds between proteins; however, the fluctuation decreased under 600 MPa, which may be due to the stabilizing effect of pressure. In addition, under high pressure, the hydrogen bonds between proteins were stable and the average hydrogen bonds were decreased only by a number of 1.5 compared to 0.1 MPa.
Changes of solvent accessible surface area of protein in molecular dynamics simulation. Changes of solvent accessible surface area of protein with simulation time from 100 to 150 ns were shown in S2 File, and the comparation under different conditions was summarized in Table 1. It could be seen that the solvent accessible surface area of protein decreased under 600 MPa, which may be due to the decrease of protein volume. The hydrophilic and hydrophobic areas were reduced simultaneously, but hydrophobic surface area of protein binding EGCG at site 1 decreased more at high pressure. The decrease of hydrophobic surface area may be due to the water molecules infiltrating into the protein molecules under high pressure, so that more hydrophobic areas are exposed to polar aqueous solution [44]. This study further showed that when the hydrophobic region bound with small molecules, it was easier to be exposed to polar water solution under high pressure, which may be due to the change of hydrophobic cavity structure caused by the combination of small molecules. Some researchers have also studied the changes of protein hydrophobicity under high pressure by experimental methods, but the results were not consistent, some thought that the hydrophobicity of protein weakened with the increase of pressure [15,43], while others reported that it increased under high pressure [45]. Changes of protein secondary structure in molecular dynamics simulation. The secondary structure changes of β-lactoglobulin under different treatment conditions were calculated with DSSP code [46]. Changes of secondary structure of protein with simulation time from 100 to 150 ns were shown in S3 File, and the comparation under different conditions was summarized in Table 2. It could be seen that under 0.1 MPa, the binding of small molecules to protein, no matter at site 1 or site 2, can stabilize the secondary structure of protein with the βsheet and α-helix increasing while the irregular coil decreasing. Which was in agreement with the result reported by Kanakis et al [4] when they studied the interaction between milk β-lactoglobulin and tea polyphenols. Under 600 MPa, the secondary structure of the protein was destroyed to some degree, especially the reduction of α-helix. Qiu [45] also found that the αhelix of myofibrillar protein decreased under high pressure and it was sensitive to pressure. This may be due to the fact that the hydrogen bond of protein was destroyed to a certain extent under high pressure, and the hydrogen bond was crucial to the stability of α-helix. Chen et al [15] found that the pressure of 400MPa could reduce the α-helix and increase the β-sheet of soybean protein. Chen [16] reported that the secondary structure of lipase changed little when the pressure was 400MPa, while some α-helix was transformed into irregular coil or β-sheet when the pressure was above 500MPa. These differences may be due to different proteins.

Binding energy between small molecule and protein using MMPBSA method
It could be seen from Table 3 that the binding between small molecule and protein mainly depended on van der Waals force, followed by electrostatic force and non-polar solvation energy (SASA energy), while the electrostatic interaction was shielded by polar solvation, the values of which were positive and played an unfavorable effect on the binding [20]. Under the pressure of 600 MPa, no matter which binding site, the binding free energies decreased significantly compared to 0.1 MPa, mainly resulting from the changes of van der Waals force. So, at

PLOS ONE
Molecular dynamics interaction between β-lactoglobulin and EGCG under high-pressure 600 MPa, the binding of small molecules with protein was not as stable as at normal pressure. From the comparison of two binding sites, the binding force of small molecules at site 1 was greater than that at site 2 under 0.1 MPa, however, it was just the opposite at 600 MPa, indicating that the change of protein structure at binding site 1 under high pressure was unfavorable to the binding of small molecules. This was related to a more dramatic decline of van der Waals force at site 1 under 600 MPa; furthermore, a reduction of 4kJ/mol for SASA energy at site 1 were also partially responsible for this change, which closely related to hydrophobic interaction. Chen [16] discussed the binding energy of lipase with ethyl acetate at 0.1, 200, 400 and 600 MPa and found that the binding energy increased at 200 MPa, but decreased at 400 and 600 MPa, the observed decline of binding energy at 600 MPa was consistent with this study.
The amino acid residues that played a greater contribution (<-0.5kJ/mol) or repulsion (>0.5kJ/mol) to the binding of protein with small molecules were shown in Fig 6. It could be seen that at 0.1 MPa, when the small molecule bound at site 1, Leu31, Pro38, Leu39, Val41, Ile71 and Ile84 were the crucial residues contributing to the binding energy, while residues Lys60, Lys69, Asp85 and Asn90 prevented the protein from binding to small molecule. Under the pressure of 600 MPa, residue Lys60 still prevented the binding of small molecule to protein strongly, and the hindrance of residues Asp85 and Asn90 was greatly reduced, while the hindrance of residue Lys69 was increased. In addition, the role of residues Glu62 and Asn109 changed from contribution to hindrance, and the contribution of residue Leu58 increased significantly under pressure. In Fig 6(A), there were 14 hydrophobic amino acids, which contributed -46.93kJ/mol binding energy in total, while in Fig 6(B), there were only 8 hydrophobic amino acids with a total energy contribution of -35.61kJ/mol, so it could be judged that the hydrophobic force between small molecule and protein decreased at 600 MPa, which was consistent with the result of SASA energy.
At atmospheric pressure, when the small molecule bound at site 2, residues Glu65, Cys160 and Ile162 contributed a lot to the binding between EGCG and protein, while residue Glu157 had the greatest role on the interaction under 600 MPa. At 0.1 MPa, there were 8 amino acid residues which had a great hindrance, while the number increased to 14 at 600 MPa, leading to the decrease of total binding energy. In particular, the contribution of residue Glu65 was significant at atmospheric pressure, but not significant at high pressure. The role of residues Glu45 and Gln59 changed from contribution to hindrance, and the contribution of residues Cys160 and Ile162 also decreased significantly at 600 MPa. In terms of hydrophobic interaction, there were only 4 and 5 hydrophobic amino acids in Fig 6(c) and 6(d), respectively, and their contribution to the total binding energy was very limited, which may be due to the fact that small molecules bound to the protein surface.

PLOS ONE
Molecular dynamics interaction between β-lactoglobulin and EGCG under high-pressure (b). It could be seen that at 0.1 MPa, the small molecule formed four hydrogen bonds with protein, two with Asp85, and one with Asn109 and Glu112 respectively. At 600 MPa, five hydrogen bonds were observed, including two with Glu62, and one with Leu58, Ala67 and Pro38 respectively. It could be seen that the hydrogen bonds formed under normal and high pressure were completely different, indicating that the binding of small molecules with protein changed greatly under high pressure. Under both pressures, the small molecule interacted hydrophobically with the 12 amino acids of the protein, but the hydrophobic environment changed a lot, at normal pressure, 10 of 12 amino acids belonged to hydrophobic amino acids, while at high pressure, only 6 hydrophobic amino acids. These results were also consistent with the conclusion above mentioned that the hydrophobic interaction was weakened under high pressure when small molecule bound at position 1 of protein. So, the binding law and environment of small molecule and protein changed greatly at 600 MPa, and thus resulting in the change of binding energy. Fig 8 showed the two-dimensional plot of the interaction between small molecule and protein at site 2 under 0.1MPa (a) and 600MPa (b). It could be seen that at 0.1 MPa, the EGCG formed six hydrogen bonds with protein, two with Glu44, and one with Gln59, Glu68, Glu158, and Cys160 respectively. At 600 MPa, seven hydrogen bonds were found, including two with Glu158, and one with Trp19, Glu44, Glu45, Gln159, and Gln159 respectively. There were three same hydrogen bonds at different pressure. More hydrogen bonds were observed between small molecule and protein at position 2 than at position 1, which may be related to the decrease of hydrogen bonds between proteins when EGCG binding to the surface of protein. Site 2 was located on the surface of protein, so the hydrophobic interaction between small molecule and amino acid residues was weak, and there were only three and four residues interacted hydrophobically with EGCG at 0.1 and 600 MPa respectively.
So, there existed van der Waals force, hydrophobic force, and hydrogen bond simultaneously when small molecule bound in hydrophobic cavity, while hydrophobic force was very small when EGCG bound on the surface of protein, mainly van der Waals force and hydrogen bond. At 600 MPa, the reduction of binding energy at site 1 was due to the decrease of hydrophobicity and van der Waals force, while the reduction of van der Waals force was the main reason for site 2.
In order to further understand the effect of high pressure on the surface structure of protein and its binding to small molecules from a macro perspective, PyMOL software was used to draw the surface structure diagram of complexes. Fig 9(A) and 9(B) showed the average surface structure of protein and the binding of polyphenols to protein at site 1 during 100-150 ns of simulation at 0.1 and 600 MPa, respectively. It could be seen that the binding site of the small molecule remained at the top of the hydrophobic cavity under 600 MPa, however, the binding pose of small molecule changed greatly, which was consistent with the results of the previous analysis on hydrogen bond and hydrophobic environment. In addition, the surface structure near the hydrophobic cavity changed obviously under high pressure, which may also affect the binding of small molecules to protein. The average surface structure of protein and the binding of polyphenols to protein at site 2 during 100-150 ns of simulation at 0.1 and 600 MPa were shown in Fig 9(C) and 9(D), respectively. It could be seen that under high-pressure, the small molecule had a large displacement on the protein surface, which may be one of the main reasons for the great changes in binding energy at 600 MPa.

Conclusion
EGCG interacted with β-lactoglobulin mainly in two sites. The binding at hydrophobic cavity (site 1) mainly depended on van der Waals force, hydrophobic interaction and hydrogen bond. Under pressure of 600 MPa, the binding energy significantly weakened mainly due to the van der Waals force and hydrophobic interaction. The interaction at protein surface (site 2) mainly depended on van der Waals force and hydrogen bond. Among which Van der Waals force was weaker while hydrogen bond interaction was a little stronger than those at site , in conjunction with the fact that the hydrophobic interaction was relatively weak at protein

PLOS ONE
surface, resulting in lower total binding force of small molecules at this site than at hydrophobic cavity under 0.1 MPa. Under the pressure of 600 MPa, the binding site of small molecules on the protein surface shifted significantly, and the van der Waals force decreased significantly, resulting in the binding energy also decreased significantly, but stronger than that of site 1 under 600 MPa. Therefore, the preferential binding site of EGCG was in the hydrophobic cavity at atmospheric pressure, while it lied in protein surface at 600 MPa. The findings in this study laid a good theoretical foundation for further improving the quality of tea milk beverage and the application of high-pressure technology in milk beverage.