Interaction of Coumarin Phytoestrogens with ERα and ERβ: A Molecular Dynamics Simulation Study.

Coumarin phytoestrogens, as one of the important classes of phytoestrogens, have been proved to play an important role in various fields of human life. In this study, molecular simulation method including molecular docking and molecular dynamics methods were performed to explore the various effects between four classical coumarin phytoestrogens (coumestrol, 4-methoxycoumestrol, psoralen and isopsoralen), and estrogen receptors (ERα, ERβ), respectively. The calculated results not only proved that the four coumarin phytoestrogens have weaker affinity than 17β-estradiol to both ERα, and ERβ, but also pointed out that the selective affinity for ERβ is greater than ERα. In addition, the binding mode indicated that the formation of hydrogen bond and hydrophobic interaction have an important effect on the stability of the complexes. Further, the calculation and decomposition of binding free energy explored the main contribution interactions to the total free energy.


Introduction
Phytoestrogens are a class of non-steroidal compounds presence in plants that bind to estrogen receptors (ER) in mammals and humans [1,2]. For the chemical structures, most of the phytoestrogens with heterocyclic polyphenols are similar to estrogens (such as 17β-estradiol) in mammals. In life process, phytoestrogens and ER combine to have a dual regulation effect, which can exert both estrogen-like effects and anti-estrogen effects. Therefore, phytoestrogens play an important role in the prevention and treatment of menopausal symptoms, breast cancer, osteoporosis and cardiovascular diseases. Due to the special mechanism of phytoestrogens, they now have been used as a natural hormone replacement to treat the estrogen-related diseases [3][4][5].
So far, among the discovered phytoestrogens, they can be divided into several categories according to their chemical structures [6]: isoflavones, coumarins, lignans, stilbenes, triterpenoids, sterols and others [7][8][9]. Phytoestrogens are widely found in various types of food diet, such as grains, fruits, and vegetables. For instance, isoflavones are mainly found in Legumes; lignans are found in grains, such as flaxseed, rye, etc., and coumarins are mainly distributed in plants including Leguminosae, Compositae, Apiaceae and Rutaceae. As we known, there are many studies that have focused on the isoflavone phytoestrogens in various fields. Some researchers also described the use of isoflavone phytoestrogens in breast cancer, osteoporosis and cardiovascular diseases as well as practice in scientific research [10][11][12][13][14][15]. Similarly, coumarin phytoestrogens are another kind of important phytoestrogens, and the like [16,17].
The above-mentioned compounds have also been proved to have a preventive or therapeutic effect on some diseases. Wu et al. [18] showed that coumestrol exerts a chemotherapeutic effect on gynecological tumors through the PI3K and ERK1/2-MAPK pathways, and is a potential novel treatment that prevents against ovarian cancer development. Zhao et al. [19] found that psoralen has the estrogen-like effect through cell experiments, but its affinity for ER is weaker than that of estradiol. It is further found that psoralen has a significant up-regulation effect on the expression level of ERα protein. Rajesh et al. [20] performed a 4D-QSAR study on coumarin derivatives and performed geometric optimization and molecular dynamics simulations for each compound. The results indicated that van der Waals interactions are important for increasing the activity of these compounds. For this, it is important to give further studies of the mechanism of coumarins at the molecular level.
In this study, four typical coumarin phytoestrogens were selected for computational studies, namely coumestrol, 4-methoxycoumestrol, psoralen and isopsoralen (See Figure 1). For the estrogenrelated compounds, they may have interaction or selection to the two estrogen receptor subtypes, ERα and ERβ. As is well known, the structures of ERα and ERβ are similar, the amino acid residues therein differ, resulting in different effects of the two subtypes on the ligands. Therefore, the four coumarin phytoestrogens were docked with the two subtypes of ER (ERα and ERβ) respectively. Then the complexes were explored by molecular dynamics simulation. The binding affinity of the ligands to the estrogen receptors and the interactions between the ligands and the receptors were investigated. In this study, 17β-estradiol, as a typical steroidal estrogen having a higher binding capacity than phytoestrogens [21,22] was used as a reference to compare with the selected coumarin phytoestrogens, and the differences between them indicate the strength of combining ability of the coumarin phytoestrogens.  Therefore, the four coumarin phytoestrogens were docked with the two subtypes of ER (ER α and ER β ) respectively. Then the complexes were explored by molecular dynamics simulation. The binding affinity of the ligands to the estrogen receptors and the interactions between the ligands and the receptors were investigated. In this study, 17β-estradiol, as a typical steroidal estrogen having a higher binding capacity than phytoestrogens [21,22] was used as a reference to compare with the selected coumarin phytoestrogens, and the differences between them indicate the strength of combining ability of the coumarin phytoestrogens.

Docking Analysis
In the present study, each of the four coumarin phytoestrogens acted as a ligand to bind to the estrogen receptor (ER α and ER β ). The docking results (docking scores) were listed in Table 1. The smaller the docking score, the more stable the complex formed. The results of docking showed that the affinity of coumarin phytoestrogens for ER α and ER β is generally lower than 17β-estradiol for both ER α and ER β [23,24]. From the docking results of coumestrol, it can be seen that the docking score of coumestrol with ER β is significantly higher than that of ER α , which proves that coumestrol has a stronger affinity for ER β than ER α . This result is consistent with the results of Zafar et al. [25]. For the corresponding system of ER α and ER β , the affinity of coumestrol and 4-methoxycoumestrol is higher than that of psoralen and isopsoralen. The results of the binding affinity of isopsoralen and psoralen for ER α are consistent with the results of the reference [26], that is, the docking score of isopsoralen is slightly lower than that of psoralen. In this study, the docking results do not make a major conclusion, but served as a reference for the subsequent research.

Molecular Dynamics Simulation Analysis
In order to evaluate the structural stability of the complexes during the simulation process, the Root Mean Square Deviation (RMSD) values of the ER α -ligands system and ER β -ligands system were separately analyzed using the Ptraj program in Amber 14. Figure 2 shows the RMSD values of the complexes in the ER α -ligands (a) and ER β -ligands (b) system. It can be seen from the Figure that the conformation of all the complexes reached a steady state after simulation of 100 ns. In the ER α system (Figure 2a), the mean of the RMSD of the stable conformations fluctuates around 0.16 nm, and the mean of the RMSD of the ER β system fluctuates around 0.14 nm (Figure 2b), indicating that the ER β system has less fluctuation than the ER α system during the MD simulation.
Similarly, in order to evaluate the changes of the ligands themselves during the MD simulation, the RMSD values of the ligands were also calculated, which were shown in Figure 3. The results showed that the RMSD values of all the ligands fluctuate within a small range, that is, the values were all below 0.1 nm, indicating that the ligands were in a stable state during the simulation process, and contribute positively to the formation of the complexes. Similarly, in order to evaluate the changes of the ligands themselves during the MD simulation, the RMSD values of the ligands were also calculated, which were shown in Figure 3. The results showed that the RMSD values of all the ligands fluctuate within a small range, that is, the values were all below 0.1nm, indicating that the ligands were in a stable state during the simulation process, and contribute positively to the formation of the complexes.  The root means square fluctuation (RMSF) values can be calculated to evaluate the stability of each amino acid in the complexes throughout the simulation. The results are shown in Figure 4. According to the curve trend in the figure, it can be seen that the fluctuation trend of the complexes in the same system are similar to the trend of the complex of ER-estradiol, indicating that the formation process of all the complexes are similar. Among them, the RMSF curve of the complexes containing coumestrol and 4-methoxycoumestrol showed even less fluctuation, compared to the RMSF curve of the complexes containing psoralen and isopsoralen, which are consistent with the docking results. The complex ERαestradiol has lower RMSF values around PRO20, GLU45, LEU76, LEU120, LEU145, ARG166， HIE213 and GLU231. The other complexes of the four coumarin phytoestrogens in the ERα system also have lower RMSF values at these similar positions, indicating that these amino acids play a positive role in the stability of the complexes. The amino acids in the ERβ system are ARG21, LEU38, MET78, ILE110, LEU142, VAL165, HIE207, and LEU223, which also performed the above-mentioned functions. The root means square fluctuation (RMSF) values can be calculated to evaluate the stability of each amino acid in the complexes throughout the simulation. The results are shown in Figure 4. According to the curve trend in the figure, it can be seen that the fluctuation trend of the complexes in the same system are similar to the trend of the complex of ER-estradiol, indicating that the formation process of all the complexes are similar. Among them, the RMSF curve of the complexes containing coumestrol and 4-methoxycoumestrol showed even less fluctuation, compared to the RMSF curve of the complexes containing psoralen and isopsoralen, which are consistent with the docking results. The complex ER α -estradiol has lower RMSF values around PRO20, GLU45, LEU76, LEU120, LEU145, ARG166, HIE213 and GLU231. The other complexes of the four coumarin phytoestrogens in the ER α system also have lower RMSF values at these similar positions, indicating that these amino acids play a positive role in the stability of the complexes. The amino acids in the ER β system are ARG21, LEU38, MET78, ILE110, LEU142, VAL165, HIE207, and LEU223, which also performed the above-mentioned functions.   Figure 5 shows the interaction of proteins with ligands. In the figure, the arrow indicates a hydrogen bond formed between ligands and amino acids, the arrow points to the donor, and the value between hydrogen bonds indicates the bond length. The results showed that 17β-estradiol forms three hydrogen bonds with the two receptors ERα, and ERβ, respectively. The hydrogen bond definition standard for this process is that the distance between the donor and the acceptor is less than 0.35nm, and the bond angle is greater than 120°. The coumestrol formed three hydrogen bonds with the ERα, and ERβ, respectively; while the 4-methoxycoumestrol formed two hydrogen bonds with them. Psoralen only formed one hydrogen bond with ERα, and no hydrogen bond was formed with ERβ. Isopsoralen did not form any hydrogen bonds with either ERα or ERβ. According to the above results, the number of hydrogen bonds formed between coumestrol and ERs is the highest, followed by 4-methoxycoumestrol, psoralen and finally isopsoralen. Among them, the phenolic   Figure 5 shows the interaction of proteins with ligands. In the figure, the arrow indicates a hydrogen bond formed between ligands and amino acids, the arrow points to the donor, and the value between hydrogen bonds indicates the bond length. The results showed that 17β-estradiol forms three hydrogen bonds with the two receptors ER α , and ER β , respectively. The hydrogen bond definition standard for this process is that the distance between the donor and the acceptor is less than 0.35 nm, and the bond angle is greater than 120 • . The coumestrol formed three hydrogen bonds with the ER α , and ER β , respectively; while the 4-methoxycoumestrol formed two hydrogen bonds with them. Psoralen only formed one hydrogen bond with ER α , and no hydrogen bond was formed with ER β . Isopsoralen did not form any hydrogen bonds with either ER α or ER β . According to the above results, the number of hydrogen bonds formed between coumestrol and ERs is the highest, followed by 4-methoxycoumestrol, psoralen and finally isopsoralen. Among them, the phenolic hydroxyl group of coumestrol forms hydrogen bonds with the residues at the two ends of the binding cavity, and the hydrogen bond interactions are greater than other ligands, which helps to improve the binding stability of coumestrol to ER. According to the references [27][28][29][30] and the comparison of the docking results, it can be seen that the hydrogen bond interactions play an important stabilizing role in the ligand-receptor binding process.

Binding Patterns Analysis
Molecules 2020, 25, x FOR PEER REVIEW 7 of 18 hydroxyl group of coumestrol forms hydrogen bonds with the residues at the two ends of the binding cavity, and the hydrogen bond interactions are greater than other ligands, which helps to improve the binding stability of coumestrol to ER. According to the references [27][28][29][30] and the comparison of the docking results, it can be seen that the hydrogen bond interactions play an important stabilizing role in the ligand-receptor binding process.
(a) To further study the binding pattern between the ligand and the receptor, the hydrophobic interaction analysis of the complexes in the last 1ns of the MD simulation was performed using Chimera1.13 software(National Institutes of Health: Bethesda, MD, USA). The results are shown in Figure 6. The amino acid residues at the pocket of the ligand-protein binding are mostly hydrophobic in nature, creating hydrophobic interactions with receptors, and further maintaining the stability of the complexes. It can be seen from the figure that all ligands are in contact with the receptor proteins, indicating the interaction between them. The hydrophobic interactions of the complexes in Figure 6 are similar, that is, the hydrophobic amino acids at the left and bottom of the binding cavity are To further study the binding pattern between the ligand and the receptor, the hydrophobic interaction analysis of the complexes in the last 1ns of the MD simulation was performed using Chimera1.13 software (National Institutes of Health: Bethesda, MD, USA). The results are shown in Figure 6. The amino acid residues at the pocket of the ligand-protein binding are mostly hydrophobic in nature, creating hydrophobic interactions with receptors, and further maintaining the stability of the complexes. It can be seen from the figure that all ligands are in contact with the receptor proteins, indicating the interaction between them. The hydrophobic interactions of the complexes in Figure 6 are similar, that is, the hydrophobic amino acids at the left and bottom of the binding cavity are greater Molecules 2020, 25, 1165 9 of 18 than other locations. This result is consistent with the reference [29], that hydrophobic interactions are the dominant force for stabilizing the binding of ligands to receptors.
Molecules 2020, 25, x FOR PEER REVIEW 9 of 18 greater than other locations. This result is consistent with the reference [29], that hydrophobic interactions are the dominant force for stabilizing the binding of ligands to receptors.

Binding Free Energy Analysis
The calculation results of the binding free energy of the ERα-ligands and ERβ-ligands system were shown in Tables 2 and 3. By comparing the calculation results of the two systems, it can be seen that the binding free energy values of the complexes in the ERβ system are generally lower than the binding free energy of the ERα system. The lower the binding free energy, the more stable the complexes formed, thus, implying that the complexes formed in the ERβ system are more stable. Furthermore, it is confirmed that the phytoestrogens have greater affinity for ERβ than ERα [23,24]. In the ERα system, according to the data in Table 2, the binding free energy value of the complex ERαestradiol is the lowest, indicating that the complex is the most stable. This result is consistent with the above conclusions that estradiol has a higher affinity for ER than phytoestrogens. For the complexes ERα-4-methoxycoumestrol and ERα-coumestrol, their binding free energy values are higher than ERα-estradiol. The structure of 4-methoxycoumestrol and coumestrol differ only in the substituents on the benzene ring, that is, one is a methoxy and the other is a hydroxyl. The free energy value of the complex ERα-4-methoxycoumestrol is slightly lower than that of ERα-coumestrol, indicating that the structure of the ligands may also have a certain effect on the stability of the complexes. The free energy values of the other two complexes, ERα-psoralen and ERα-isopsoralen are similar and somewhat higher, indicating that their stability is not as good as the former two. The complexes in the ERβ system have the same trend as the ERα system, indicating that among the four Figure 6. Hydrophobic interaction of complexes, the color from red to blue indicates a decrease in hydrophobic potential.

Binding Free Energy Analysis
The calculation results of the binding free energy of the ER α -ligands and ER β -ligands system were shown in Tables 2 and 3. By comparing the calculation results of the two systems, it can be seen that the binding free energy values of the complexes in the ER β system are generally lower than the binding free energy of the ER α system. The lower the binding free energy, the more stable the complexes formed, thus, implying that the complexes formed in the ER β system are more stable. Furthermore, it is confirmed that the phytoestrogens have greater affinity for ER β than ER α [23,24]. In the ER α system, according to the data in Table 2, the binding free energy value of the complex ER α -estradiol is the lowest, indicating that the complex is the most stable. This result is consistent with the above conclusions that estradiol has a higher affinity for ER than phytoestrogens. For the complexes ER α -4-methoxycoumestrol and ER α -coumestrol, their binding free energy values are higher than ER α -estradiol. The structure of 4-methoxycoumestrol and coumestrol differ only in the substituents on the benzene ring, that is, one is a methoxy and the other is a hydroxyl. The free energy value of the complex ER α -4-methoxycoumestrol is slightly lower than that of ER α -coumestrol, indicating that the structure of the ligands may also have a certain effect on the stability of the complexes. The free energy values of the other two complexes, ER α -psoralen and ER α -isopsoralen are similar and somewhat higher, indicating that their stability is not as good as the former two. The complexes in the ER β system have the same trend as the ER α system, indicating that among the four phytoestrogens, 4-methoxycoumestrol has the highest affinity for both ER α and ER β , followed by coumestrol, psoralen and finally isopsoralen.  According to the energy contribution data of van der Waals energy, electrostatic interaction, polar solvation energy and non-polar solvation energy, one can see that the van der Waals energy has the largest contribution (the energy value is the smallest). In addition, during the formation of the complexes, the van der Waals energy, electrostatic interaction energy and non-polar solvation energy in the ER α -ligands and ER β -ligands systems are negative, indicating that these three effects are beneficial to the combination of ligands and acceptors; the value of the polar solvation energy is positive, indicating that it hinders the formation of the complexes.
In order to compare our results with the experimental results, we searched the clinical assays reported in ChEMBL (https://www.ebi.ac.uk/chembl/) for these compounds. However, there are only two compounds that can be found, that is, the IC 50 values of ER α -coumestrol and ER α -estradiol as well as ER β -estradiol. The results are shown in Table 4. As can be seen in Table 4, the IC 50 values of estradiol binding to ER α and ER β are both 2 nM, which are lower than the IC 50 values of other complexes, indicating that the affinity of estradiol is higher than that of coumarin phytoestrogens. In addition, the IC 50 value of ER α -coumestrol is higher than ER β -coumestrol, which is also consistent with our calculation results of binding energy.

Binding Free Energy Decomposition
All amino acids of the protein were analyzed, in order to further understand the contribution of each amino acid to the binding free energy during complexes formation. Amino acids that contribute greatly to the binding free energy (more than −1 kcal/mol) were defined as key amino acids, and energy analysis was performed on these key amino acids. The results of the analysis were shown in Figures 7 and 8. In the ER α system, the number of key amino acids in the complexes is 9, 9, 5, 4, and 10, respectively (Residues name and sequence were shown in Figure 7). The number of key amino acids in the ER β system is 9, 9, 5, 2, and 10 (Residues name and sequence were shown in Figure 8). The larger the number of key amino acids, the more stable the complexes formed. The above results indicate that the complexes ER α -estradiol and ER β -estradiol are the most stable of all complexes (The two compounds have the largest number of key amino acids, 10). The number of key amino acids in complexes ER α -4-methoxycoumestrol (9), ER β -4-methoxycoumestrol (9), ER α -coumestrol (9) and ER β -coumestrol (9) are same and greater than the number of key amino acids in the complexes ER α -psoralen (5), ER β -psoralen (5), ER α -isopsoralen (4), and ER β -isopsoralen (2), which imply that the former are more stable than the latter. This is consistent with the above analysis. Almost all of these key amino acids are located around the binding cavity of the ligands and the receptor proteins, and the fluctuation is small during the MD simulation process, and the contribution to the binding free energy is large.
In the energy analysis process of the key amino acids, van der Waals energy and non-polar solvation energy were defined as hydrophobic interactions (black part in the figure), and electrostatic energy and polar solvation energy were defined as electrostatic interactions (grey part in the Figure). According to the contribution of key amino acids in Figures 7 and 8, the hydrophobic interactions are mostly to favor the formation of complexes (the energy is negative), while the electrostatic interactions are more likely to be detrimental to the formation of complexes (the energy is positive). In summary, the hydrophobic interactions promote the binding of ligands to receptors.   In the energy analysis process of the key amino acids, van der Waals energy and non-polar solvation energy were defined as hydrophobic interactions (black part in the figure), and electrostatic energy and polar solvation energy were defined as electrostatic interactions (grey part in the Figure). According to the contribution of key amino acids in Figures 7 and 8, the hydrophobic interactions are mostly to favor the formation of complexes (the energy is negative), while the electrostatic interactions are more likely to be detrimental to the formation of complexes (the energy is positive). In summary, the hydrophobic interactions promote the binding of ligands to receptors.

Protein and Ligand Preparation
The initial structures of the estrogen receptors ERα (PDB ID: 1GWR) and ERβ (PDB ID: 3OLS) were obtained from the RCSB Protein Data Bank [33]. The ligand binding domain of ERα and ERβ were extracted from the initial structures for the subsequent studies. Preparation of the receptor proteins (ERα and ERβ) were performed by the Protein Preparation Wizard module of the Schrodinger

Protein and Ligand Preparation
The initial structures of the estrogen receptors ER α (PDB ID: 1GWR) and ER β (PDB ID: 3OLS) were obtained from the RCSB Protein Data Bank [33]. The ligand binding domain of ER α and ER β were extracted from the initial structures for the subsequent studies. Preparation of the receptor proteins (ER α and ER β ) were performed by the Protein Preparation Wizard module of the Schrodinger 2015-2 software (Schrodinger, Inc., NY, USA) [34]. The process of preparations included the removal of water molecules, the addition of hydrogen atoms, and the optimization of the conformations under the Schrodinger's OPLS_2005 force field [35].
The structure files in SDF format for all the ligands were obtained from the Pubchem database [36]. The ligands were ionized and optimized using the Ligprep Wizard module of the Schrodinger software (Schrodinger, Inc., NY, USA) [37], and the preparation process considered all possible conformations. Finally, the most suitable ligand conformation was selected for molecular docking.

Molecular Docking Studies
First, a small molecule ligand was used as a center to prepare a 3D space grid with the residues surrounding the binding site. Then, the prepared ligands and the receptor proteins were docked using the Glide 6.7 docking procedure in Schrodinger 2015-2 (Schrodinger, Inc., NY, USA) [38][39][40]. Due to the existence of certain deviations in various docking methods, the docking results were only used for reference and did not make the key conclusions. Finally, the complexes formed by docking were used for molecular dynamics simulation studies.

Molecular Dynamics Simulations
The molecular dynamics (MD) simulation process was performed using the Amber14 program (University of California: San Francisco, CA, USA) [41,42]. The parameters of the system were generated by Antechamber module in Amber program. The restrained electrostatic potential (RESP) was used to describe the partial atomic charges. The FF03 force field [43] in Amber was used to describe the molecular parameters of the protein, and the molecular parameters of the ligand were described using the GAFF force field [44]. The complexes were hydro-treated using the tleap module. The complexes were dissolved in a TIP3P water molecule box with a distance of at least 12Å from the edge of the water box. For the ER α system, five sodium cations were added to neutralize the negative charge [45], making the system electrically neutral. Two sodium cations were also added to the ER β system for the same purpose.
The purpose of energy minimization was to make the system more balanced. In the optimization process of each step, the steepest descent method of 3000 steps was firstly performed, and then the conjugate gradient method of 2000 steps was optimized. The first optimization was performed to minimize the solvent; then to constrain the counter ions for energy minimization. Finally, the entire system was minimized without restrictions. The interception distance of a long-range Coulombic interaction was set to 1 nm using the Particle Mesh Ewald (PME) method [46]. The SHAKE algorithm was used to limit the key length, and the integration step in the MD simulation process was set to 2 fs. Subsequently, the temperature rise dynamics of the system was run. The temperature was raised from 0 K to 300 K in 50 ps, and the collision frequency was 2 ps −1 . Then, the system was subjected to an unrestricted equilibrium dynamics of 500 ps. Finally, all systems continue to run for 100 ns in a normal temperature, constant pressure NPT system, and the trajectory calculation interval was 1 ps.
The Root Mean Square Deviation (RMSD) and Root Mean Square Fluctuation (RMSF) were used to describe the structural features, which were calculated by the Ptraj module [47] of Amber14 (University of California: San Francisco, CA, USA).

Binding Free Energy Calculations and Energy Decomposition
The calculation of the binding free energy of all complexes was obtained by the MMPBSA.py script [48] in Amber14 software (University of California: San Francisco, CA, USA). This process extracted a snapshot from the trajectory file after the MD simulation and then performed the calculation. The binding free energy of the complexes was calculated as follows: ∆G MM = ∆G int + ∆G ele + ∆G vdw (2) Among them, the binding free energy (∆G bind ) was composed of the binding free energy of the complex (∆G comp ), the receptor protein (∆G pro ) and the ligand (∆G lig ). The gas phase binding energy (∆G MM ) was composed of internal energy (∆G int ), electrostatic interaction energy (∆G ele ) and van der Waals energy (∆G vdw ); the solvation free energy (∆G solv ) was calculated from the polar (∆G pol ) and non-polar (∆G nopol ) solvation free energy. For entropy change (T∆S), the change of entropy causes by conformational change after the ligand binds to the receptor. The purpose of this paper is only to compare the binding energy between several systems, and the calculation process of entropy change is more complicated, so the contribution of entropy change to binding free energy is neglected.
To form a complex by the ligand binding to the protein, each amino acid in the system contributed differently to the binding free energy. For the purpose of this analysis, the free energy decomposition process was used by the MM/GBSA free energy decomposition module in Amber14 software (University of California: San Francisco, CA, USA) to decompose the free energy and then explore the contribution of each amino acid. Among them, the contribution of three kinds of energy was investigated, that is, the van der Waals energy, the electrostatic interaction energy as well as the solvation energy.

Conclusions
In this paper, four coumarin phytoestrogens were used as ligands to dock with ER α and ER β , and the binding mode along with the binding free energy between ligands and receptors were explored by molecular dynamics simulation, respectively. The results showed that the affinity of the complex 4-methoxycoumestrol was slightly stronger than that of coumestrol, and was significantly higher than that of psoralen and isopsoralen. The four coumarin phytoestrogens had greater affinity and selectivity to ER β than ER α . In addition, the calculation and decomposition of the binding free energy pointed out that van der Waals energy and non-polar solvation energy play a beneficial role in the formation of complexes. In general, we hope that the results of this study can provide some help for further exploration of the weak interaction between the phytoestrogens and estrogen receptors.