The most representative (-)-Epicatechin metabolites were studied by DFT to determine the physicochemical properties of each one. The structures obtained after energy minimization are shown in Fig. 1. (-)-Epicatechin is a flavonoid, and as such, it shares a chemical core of benzo-γ-pyrone that consist in two benzene rings and a dihydropyran heterocycle (named A, B and C rings). In specific, (-)-Epicatechin is a flavan-3-ol characterized by lacking the keto group at the position 4 (Fig. 1A). (-)-Epicatechin has a biotransformation process after oral intake in humans, generating different classes of metabolites; some conserve the benzo-γ-pyrone chemical core with O-glucuronide group addition, like (-)-Epicatechin-3′-O-glucuronide (E3′G) (Fig. 1B), or a sulfate group addition in (-)-Epicatechin-3′sulfate (E3′S) (Fig. 1C). Also, the structure accepts more than one group like 3′-O-metyl-(-)-Epicatechin-5-sulfate (3′ME5S), and 3′-O-metyl-(-)-Epicatechin-7-sulfate (3′ME7S) (Figs. 1D and 1E). In addition, the benzo-γ-pyrone chemical core could be modified to C5 side chain ring fission and form 5-(4′-hidroxyphenyl)-γ-verolactone-3′-O-glucuronide (γVL3′G), and 5-(4′-hidroxyphenyl)-γ-verolactone-3′-sulfate (γVL3′S) (Figs. 1F and 1G).
A molecular potential analysis was applied to (-)-Epicatechin and metabolites in order to describe the possible electron density distribution (Fig. 2). Molecular electrostatic potential (MEP) is a descriptor that allow the prediction of molecular sites that may participate on noncovalent interactions [42]. The molecular sites with disposition of electrophilic attacks represented by blue surface are positioned in hydrogen atoms; in contrast, red surface represent molecular sites that could be part of nucleophilic attack and is located on oxygen atoms. In specific, sulfate group on E3′S, 3′ME5S, 3′ME7S, and γVL3′S could induce a more characteristic reactivity to these metabolites whereas a glucuronide group on E3′G, and γVL3′G suggest that these metabolites have more -OH reactive groups that could participate in more H-bonds formation.
Frontier molecular orbital (FMO) was calculated for metabolites and (-)-Epicatechin in order to determine the electron occupation throughout the molecule, and to predict the molecular sites that could participate donating or accepting electrons on a reaction [43]. FMO is configured by the highest occupied molecular orbital (HOMO) that revels the tendency to donate electrons and the lowest unoccupied molecular orbital (LUMO), with the ability to accept electrons [43]. HOMO orbitals (Fig. 3) were located with more frequency on A and B rings in the benzo-γ-pyrone chemical core of the (-)-Epicatechin and metabolites. HOMO energy values were between − 134 kcal/mol to -142 kcal/mol for (-)-Epicatechin and metabolites, except for γVL3′S, where the homo value was − 161.897 kcal/mol and the HOMO orbital was located on sulphate group (Fig. 5A). This shows that the C5 side chain ring fission of γVL3′S, induces significant alteration in physicochemical properties with respect to (-)- Epicatechin.
LUMO molecular orbitals representations are presented in Fig. 4 and values are plotted on Fig. 5. LUMO orbitals were located on the C- ring of the benzo-γ-pyrone chemical core of the flavanols. Higher values of LUMO are predicted for metabolites with sulfate group such as E3′S, 3′ME5S, 3′ME7S, and γVL3′S, suggesting that sulfate group induce on metabolites a predisposition to accept electrons. The energy difference between the HOMO and LUMO (GAP energy of HOMO-LUMO (ELUMO−HOMO)) can be used to determine kinetic stability and chemical reactivity of a molecule. Hence, a short gap suggests low kinetic stability and high chemical reactivity [44, 45]. The lowest GAP energy (Fig. 5A) was calculated for E3′S (114.39 kcal/mol) and the highest for γVL3′S (132.21 kcal/mol), also lower GAP energies were calculated for 3′ME5S and 3′ME7S compared to (-)-Epicatechin, suggesting that sulfate group increase the reactivity of the metabolites that could be favorable for their interaction with GPER. In contrast, the glucuronide group on E3′G and γVL3′G have higher GAP energy values, suggesting that glucuronide increase the stability of the metabolites. Then, this higher kinetic stability is associated with a more stable interaction when the molecule associates with specific proteins.
Then, reactivity indexes were calculated for (-)-Epicatechin and metabolites (Fig. 5 and Table 1). We observed that reactivity calculated by GAP energy has a correlation with the physicochemical properties calculated, specifically, Global Hardness which measures the degree that a certain molecule is unable to perform a charge transfer. Lower values of global hardness were calculated for those metabolites with a sulfate group such as E3′S, 3′ME5S, and 3′ME7S and were compared to (-)-Epicatechin, showing that these metabolites have higher ability to transfer charge, suggesting a higher chemical reactivity and lowest kinetic stability.
The chemical potential serves as an indicator of how likely electrons in a molecule are to transfer to other molecules[46–48]. A higher chemical potential implies greater reactivity. This descriptor was favorable for E3′S and showed acceptable values for other metabolites (see Fig. 5-B). In contrast, γVL3′S had the highest GAP energy and chemical potential among all the metabolites, along with significantly different values for other descriptors. This suggests that the chemical similarity between (-)-Epicatechin and γVL3′S is low, suggesting a probably reduced ability to interact with GPER.
Another set of physicochemical descriptors that can hint at a molecule's ability to interact in a complex are electron donor and acceptor properties. In this context, E3′S exhibited the highest values for both electron donation and acceptance (as shown in Fig. 5-B). Additionally, 3′ME5S is noteworthy because it has a lower electron donor index but a high electron acceptor index when compared to (-)-Epicatechin and 3′ME7S (Fig. 5-B). This observation, along with lower GAP and Hardness values and higher chemical potential, suggests that 3′ME5S possesses characteristics that could enhance its interactions with GPER.
Regarding E3′G, its chemical properties align closely with those of (-)-Epicatechin in terms of chemical potential, electron donor, and acceptor indexes. However, it exhibits an increased GAP energy (Fig. 5), implying that E3′G might have a similar interaction with GPER but potentially enhances the stability of the interaction.
A PCA analysis of metabolites using molecular descriptor information suggest that the E3′G metabolite exhibit the closest physicochemical descriptors to (-)-Epicatechin, in contrast γVL3′S has the lowest chemical similarity to (-)-Epicatechin (Fig. 5C and Table 1).
Table 1
Reactivity indexes were calculated for (-)-Epicatechin and metabolites. All the values are in Kcal/mol.
Molecular descriptor / Molecule | (-)-Epicatechin | E3′G | E3′S | 3′ME5S | 3′ME7S | γVL3′G | γVL3′S |
HOMO energy | -137.081 | -137.625 | -140.496 | -136.123 | -134.664 | -141.277 | -161.897 |
LUMO energy | -17.724 | -17.623 | -26.105 | -20.535 | -18.091 | -16.502 | -29.681 |
GAP | 119.357 | 120.001 | 114.391 | 115.587 | 116.573 | 124.775 | 132.216 |
Ionization potential | 137.081 | 137.625 | 140.496 | 136.123 | 134.664 | 141.277 | 161.897 |
Electroaffinity | 17.724 | 17.623 | 26.105 | 20.535 | 18.091 | 16.502 | 29.681 |
Chemical potential | -77.403 | -77.624 | -83.301 | -78.329 | -76.377 | -78.889 | -95.789 |
Global Hardness | 59.678 | 60.001 | 57.195 | 57.794 | 58.286 | 62.387 | 66.108 |
Global Softness | 3298.988 | 3281.268 | 3442.212 | 3406.573 | 3377.768 | 3155.741 | 2978.126 |
Electronegativity | 77.403 | 77.624 | 83.301 | 78.329 | 76.377 | 78.889 | 95.789 |
Electrophilicity Index (w) | 50.196 | 50.212 | 60.661 | 53.080 | 50.042 | 49.878 | 69.398 |
w + Electron Aceptor | 18.954 | 18.900 | 26.160 | 21.140 | 19.139 | 18.232 | 29.767 |
w- Electron Donator | 96.357 | 96.524 | 109.460 | 21.140 | 95.516 | 97.121 | 125.556 |
Electrophilicity Net | 115.311 | 115.424 | 135.620 | 120.609 | 114.655 | 115.353 | 155.323 |
After DFT characterization of (-)-Epicatechin metabolites, we evaluated its binding affinity to GPER protein. To address this, we first employed rigid molecular docking (Fig. 6-A). (-)-Epicatechin and metabolites exhibited favorable binding energies between − 7.0 to -9.0 kcal/mol, in fact, (-)-Epicatechin metabolites showed better binding affinity prediction to GPER that (-)-Epicatechin, E3′G exhibit the best affinity to GPER in contrast to γVL3′S, that was predicted to have less binding affinity. Then we performed a flex molecular docking (Fig. 6-B), and the results showed binding affinity energies between − 7.5 kcal/mol to -9.5 kcal/mol. Here γVL3′G and E3′G have the most favorable binding affinities. This suggests that glucuronide group induce better affinity because it has better LUMO energies (Fig. 5-A) associated to greater number of OH groups and a grater polarity. We observed variations in how strongly 3′ME5S and 3′ME7S bind. These molecules have the same atoms, but their sulfate group is in different positions, resulting in distinct physicochemical properties based on DFT calculations (e.g., GAP energies and electron donator index). These properties favor 3′ME5S, leading to a better binding affinity prediction for 3′ME5S compared to 3′ME7S. Furthermore, γVL3′S, which differs significantly in physicochemical characteristics from (-)-Epicatechin, demonstrates a poorer binding affinity to GPER compared to other metabolites that have a closer chemical resemblance to (-)-Epicatechin.
To further study the probable interaction of (-)-Epicatechin and its metabolites with GPER, we performed a 100ns molecular dynamics simulation of GPER in a membranal context. First, we evaluate the distance between the atoms, considering of the structure obtained on the time of simulation by root mean square deviation (RMSD). GPER protein RMSD in all complexes (Figure S 1-A) exhibited stabilization around 2.5nm. In contrast, the RMSDs for of (-)-Epicatechin metabolites take different configurations during the simulation, compared with the stability of (-)-of Epicatechin (Figure S 1-B). Figure 7-A shows the RMSD value considering the alpha-carbon of the protein and the specific ligand (named complex), here the complex stability exhibit differences depending on the metabolite, GPER/(-)-Epicatechin show lowest RMSD values and GPER/γVL3′S the highest; GPER/γVL3′G remains stable during the simulation but the last 20 ns turn to unstable interactions, the other complex metabolites (E3′G, E3′S, 3′ME5S, and 3′ME7S) with GPER maintains stable during the simulation around 0.3nm. RMSD analysis suggest that in contrast to (-)-Epicatechin, the metabolites could turn into different binding poses when interacts with GPER.
The molecular dynamics simulation also was evaluated in terms of radius of gyrate (Rg). Here we evaluated the fluctuation of the atoms distance to the center mass during the simulation; Rg values give information about the compactness and stability of the structure during the simulation. Rg values of GPER protein (Figure S 1-C) show lower fluctuation (0.017 nm average of standard deviation Rg values) during the simulation and that compactness was replicated on Rg values calculated for GPER complexes (0.019 nm average of standard deviation Rg values) (Fig. 7-B).
The average of GPER / metabolites Rg values are 2.322 nm with average standard deviation of 0.019 nm. The lowest values were obtained for GPER/E3′G (2.29611 ± 0.0223 nm) and the more stable values for GPER/(-)-Epicatechin. In contrast to GPER Rg (0.0189nm standard deviation of average Rg values) and complex Rg values (0.0158nm standard deviation of average Rg values), the calculation of Rg ligand show differences between metabolites (0.0395nm standard deviation of average Rg values), here the more stable values are associated to (-)-Epicatechin (0.372 ± 0.0029 nm) and 3′ME5S (0.430 ± 0.0065 nm), and the highest values for E3′G (0.482 ± 0.0156 nm). Rg values of metabolites suggest different binding poses associated to each metabolite.
Root means square fluctuation (RMSF) was calculated during the molecular dynamic simulation (Fig. 7-C). By doing this, it is possible to determine the flexibility associated to specific residue by the determination of its movement during the simulation. The flexibility increased on the intracellular domain for GPER/3′ME5S and GPER/(-)-Epicatechin. This region is reported to be important for GPER activation [30]. Higher flexibility was calculated for GPER/γVL3′G, and GPER/γVL3′S on H2, and GPER/E3′S on H7 are also regions reported to be important for GPER activation [30].
Molecular dynamics simulations were used to analyze GPER and (-)-Epicatechin metabolites structures. Based on RMSD and Rg values, it was observed that (-)-Epicatechin metabolites have multiple binding sites with the GPER protein, suggesting that these metabolites bind slightly differently compared to the predicted binding of (-)-Epicatechin.
The complexes underwent a clustering analysis to ascertain the most stable structure for each complex. The results indicate that the complexes exhibit a deviation ranging from 0.1 to 0.6 nm (Figure S2). The most populated GPER/metabolite cluster was studied using NCI software to identify and describe non-covalent interactions responsible for stabilizing the binding. To characterize the non-covalent interaction between metabolites and GPER protein a sign(λ2)ρ determination were used. These interactions encompass both favorable and unfavorable interactions. To distinguish between them, the sign of the second density Hessian eigenvalue (λ2) multiplied by the density is employed (sign(λ2)ρ). This value effectively characterizes the strength of interactions based on density and their nature determined by the sign of the second eigenvalue. Attractive and repulsive interactions are recognized in regions where λ2 < 0 and λ2 > 0, respectively. Weak van der Waals interactions are indicated by λ2 ≈ 0. To visualize these regions, a color code is often used based on sign(λ2)ρ: blue for strong attractive interactions (λ2 < 0), green for weak van der Waals interactions (λ ≈ 0), and red for strong repulsive interactions (λ2 > 0)[37]. Additionally, LigPlot software was used to create 2D representations (see Figs. 8 and 9).
In the case of the GPER/(-)-Epicatechin complex, it was found to be stabilized by four hydrogen bonds with Glu115, Ser62, His307, and Pro303, as well as hydrophobic interactions involving Asn310, Leu119, Gly306, Phe278, and Phe206.
The complex GPER/ E3′G is bound by 4 H-Bonds interactions with Thr40, Asn44, Glu203 and Glu48, and also by hydrophobic interactions with Phe206, Phe278, Val309, Thr305, Gly306, Phe208, His302, Arg299, Tyr55, Ala43 and Leu59. The H-bonds with Asn44 and Glu203 were facilitated by a glucuronide group, where the LUMO determination suggests the ability to this molecule to accept electrons. Also, interactions with phenylalanines 206, 208, and 278 are recognized interactions[49] that promote π-π interaction that is key on the aromatic cavity of GPER[29]. The complex of GPER/E3′S is stabilized by 3 H-bonds with Asp125, Cys205 and Cys207, also by hydrophobic interactions with Cys130, Tyr123, Ala204, Leu129, Met133, Phe206. The sulfate group in contrast with glucuronide group in E3′G only form the only H-bond with Asp125, but this binding allows the formation of two H-bonds with cysteines 207 and 205. These residues, Cys 130 and 207, are reported to generate a disulphide bond and this bond participate on structural changes that induce a rearrangement in the intracellular domain associated to GPER agonist properties[49].
Moreover, the interaction between 3′ME5S and GPER is stabilized through three hydrogen bonds. Two of these bonds facilitate the sulfate group's interaction with Glu203, while the third forms a hydrogen bond with Ala41. Additionally, the complex is reinforced by hydrophobic interactions involving Tyr55, His307, Leu59, Gly306, Phe278, His302, Phe206, Ala43, Thr201, and Leu42. The dual hydrogen bond with Glu203 contributes significantly to the overall stability of the complex and concurrently promotes π-π interactions with Phe206 and Phe278. The complex of GPER/3′ME7S like 3′ME5S maintains the double H-bond with Glu203 and sulfate group, but the H-bond with Ala41 is not reported. This metabolite is also stabilized by interactions with Phe206, Val309, Gly306, Ile279, Glu275, Phe278, His302, Pro303, and Ala43. The γVL3′G metabolite stabilizes its binding with GPER through 2 H-bonds by the glucuronide group with Cys207 and Glu275; and also by hydrophobic interactions with Phe208, Phe206, Thr40, His302, Ala43, His307, Asn310, Phe278 and Val309. The glucuronide group participates on important interactions with GPER by the presence of C5 side chain ring fission on the characteristic benzo-γ-pyrone chemical core. Finally, the complex of GPER/γVL3′S is stabilized by hydrophobic interactions with Asn118, Tyr123, Ile126, Leu129, Cys130, and Ile114; and by an H-bond with Glu115.
We also evaluated the GPER binding properties of (-)-Epicatechin metabolites by MMBGSA analysis. All the (-)-Epicatechin metabolites exhibit better binding affinity to GPER in comparation with (-)-Epicatechin (Fig. 10). Among the metabolites considered, 3′ME5S stands out with the highest binding affinity at -37.61 kcal/mol. This superior affinity is correlated with a more extensive array of non-covalent interactions, as evident from the sign(λ2)ρ distribution illustrated in Fig. 9-A. Specifically, the distribution pattern is linked to the formation of a double hydrogen bond in the sulfate group. Furthermore, the non-covalent interactions surrounding van der Waals energies are evenly dispersed around the molecule, as depicted in the isosurface of Fig. 9-A. In contrast, E3′G, E3′S, and 3′ME7S display binding energies of -34.02, -28.95, and − 24.06 kcal/mol, respectively. The sign(λ2)ρ plot for these metabolites indicates a reduction in density (depicted by red dots) of repulsive interactions with positive values of sign(λ2)ρ (refer to plots in Fig. 8-B, -C, and 9-B), as compared to the (-)-Epicatechin plot. In that case the glucuronide and sulfate groups enhance the formation of non-covalent interactions in order to van der Waals energies (see isosurfaces plotted on Fig. 8-B, -C, and 9-B) where the interaction as enhanced showing by the isosurfaces distribution around the metabolites.
In γVL3′G, the H-Bonds interactions on glucuronide group participates on attractive interactions and enhance van der Waals interactions (see Fig. 9-C) with GPER and in comparation to (-)-Epicatechin, that show a binding affinity energy of -26.98 kcal/mol. Finally, in γVL3′S complex exhibit binding affinity of -21.17 kcal/mol, that probably associated to attractive interactions where the H-bond formation induce a sign(λ2)ρ units around − 0.06 if it was compare with other metabolites.
To gain inside of relevant molecular descriptors to improve the GPER interaction based on (-)-Epicatechin or metabolites interactions a QSAR methodology was applied. Within (-)-Epicatechin and metabolites as a set, there were 431 3D-descriptors and 1,444 2D-descriptors calculated with PaDEL, while 1,665 descriptors with Mordred. For the QSAR analysis, we selected the regression model with the best adjusted R2 for each molecular descriptor dataset to ensure robustness and accuracy. Remarkably all the models showed high adjusted R2, indicating a high level of predictive accuracy. The optimal regression model for Mordred had an adjusted R2 of 0.9998 (Figure S3A):
Model 1: MMGBSA: 416.81467 + (-205.30923) * RNCG + (-54.65308) * MWC05 + (2.71864) * ATSC6p
This model integrates the Relative Negative Charge (RNCG), Walk Count of Leg-5 (MWC05), and Centered Moreau-Broto Autocorrelation of Lag 6 weighted by Pauling EN (ATSC6p) as pivotal molecular descriptors. Notably, RNCG stands out as the two least favourable MMGBSA values were associated with the two highest RNCG values ((-)-E and γVL3′G). As RNCG is calculated by taking the charge of the most negative atom and dividing it by the total negative charge, this could suggest that a wider distribution of negative charge could enhance the binding of molecules within the system[50]. PaDEL 3D showcased an adjusted R2 of 0.9970 in its most effective model (Figure S3B):
Model 2: MMGBSA: -19.18640 + (-0.00494) * TDB10m + (-7.04049) * RDF105p + (1.64234) * RDF115e
This model encapsulates 3D topological distance-based autocorrelation - lag 10 / weighted by mass (TDB10m), Radial distribution function − 105 / weighted by relative polarizabilities (RDF105p), and by relative Sanderson electronegativities (RDF115e) as essential descriptors. This model relies on the spatial arrangement of the ligand’s features. PaDEL 2D, demonstrated the highest adjusted R2 of 0.9999 for its optimal model (Figure S3C):
Model 3: MMGBSA: -10.90223 + (-0.04025) * ATS7s + (1.52915) * ATSC6p + (2.64078) * SssCH2
Within this 2D model, the identified descriptors include Broto-Moreau autocorrelation - lag 7 / weighted by I-state (ATS7s), Centered Broto-Moreau autocorrelation - lag 6 / weighted by polarizabilities (ATSC6p), and the Sum of atom-type E-State: -CH2- (SssCH2). This model highlights the importance of a specific spatial distribution within this metabolite. At the same time, SssCH2 stands out as the only descriptor directly related to chemical groups. Moreover, as 3′ME5S, the molecule with the highest affinity, possess the smallest positive value. In this way, as SssCH2 the sum of the E-States for carbon atoms with two hydrogens attached and two single bonds to a nonhydrogen atoms, a lower positive E-states of the CH2 could enhance the interactions within the system[51]. Therefore, Metabolites interaction in the system are influenced by a particular spatial distribution, a wider negative charge repetition, and a lower positive E-State.
In summary, for Mordred, the optimal model highlighted the importance of descriptors like RNCG, MWC05, and ATSC6p, suggesting a wider distribution of negative charge enhances binding. PaDEL 3D emphasized the significance of spatial arrangement with descriptors TDB10m, RDF105p, and RDF115e. PaDEL 2D revealed the importance of specific spatial distribution, wider negative charge repetition, and lower positive E-State for metabolite interactions with GPER. In conclusion, our findings provide insights into molecular features crucial for GPER interactions, offering a foundation for further optimization and understanding of the underlying mechanisms.