A molecular dynamics perspective into estrogen receptor inhibition by selective flavonoids as alternative therapeutic options

Abstract

Zearalenone is an estrogenic mycotoxin which is a common food contaminant and has been implicated in increasing the incidence of carcinogenesis and other reproductive health ailments through the estrogen receptor alpha (ERα) pathway. Competitive ERα blockers such as 4-Hydroxytamoxifen (OHT), are synthetic FDA approved drugs which, albeit being an effective anticancer agent, induces life altering side effects. For this reason, there is an increased interest in the use of naturally occurring medicinal plant products such as flavonoids. This study aimed to identity flavonoid ERα inhibitors and provide insights into the mechanism of inhibition using computational techniques. The Molecular Mechanics/Generalized Born Surface Area calculations revealed that quercetrin, hesperidin, epigallocatechin 3-gallate and kaempferol 7-O-glucoside out of 14 flavonoids had higher binding affinity for ERα than OHT. The structural analysis revealed that the binding of the compounds to the receptor lead to dynamic alterations, which induced conformational shift in the structure and orientation of the receptor resulting in stabilised, compact and low energy systems. The results of this study provide imperative information that supports the use of flavonoids in the inhibition of ERα to prevent or ameliorate the consequential adverse effects associated with zearalenone exposure.

Communicated by Ramaswamy H. Sarma

Keywords:

• Zearalenone
• estrogen receptor α
• flavonoids
• inhibitor
• molecular dynamic simulations

Introduction

Mycotoxins are secondary metabolites produced by fungi that can be toxic to man and animals and are ubiquitously found in the environment. Common ones include aflatoxin, fumonisins, ochratoxin A, citrinin, ergot alkaloids, trichothecenes, patulin, and zearalenone (Andersen et al., 2018). Unlike the other mycotoxins, zearalenone exhibits estrogenicity and hence referred to as a mycoestrogen (Andersen et al., 2018).

Zearalenone (ZEA) is a non-steroidal mycoestrogen produced by several Fusarium fungi species especially F. graminearum and F. cuolmorum (Bennett & Klich, 2003; Kuiper-Goodman et al., 1987). These species of fungi are commonly found in human and animal foodstuffs such as cereal crops, dried fruits, spices, wine, beer and meat products. Contamination of cereal crops with ZEA generally occurs during pre-harvest and persists thereafter due to the mycotoxin’s high resistance to elevated food processing temperature and pressure conditions. This mechanism allows the bioaccumulation of ZEA in humans and animals that consume the contaminated produce (Andersen et al., 2018).

Chemically described as 6-(10-hydroxy-6-oxo-trans-1-undecenyl)-β-resocyclic acid lactone, the structural similarity of ZEA (Figure 1B) to the endogenous estrogens such as 17β-estradiol (Figure 1A) plus the high binding affinity to the estrogen receptors are the reasons for its harmful estrogenic activities. ZEA has been implicated in carcinogenesis (breast cancer, oesophageal cancer, cervical cancer, prostate cancer) through estrogen signaling which influences DNA synthesis and cell reproduction (Ahmed et al., 2017; Claeys et al., 2020; Kowalska et al., 2018; Luo et al., 1990). Animal studies have reported ZEA induced infertility, spontaneous abortions and still births. ZEA also modulates the synthesis of hormones (estradiol, cortisol, testosterone and progesterone) by binding to estrogen receptors (ERα and ERβ), leading to the transcriptional activation of estrogen responsive genes (Andersen et al., 2018).

Figure 1. Chemical structures of the ERα ligands. (A) Estradiol, the endogenous stimulatory ERα-ligand (ChemSpider ID: 5554). (B) ZEA, the stimulatory mycoestrogen (ChemSpider ID: 4444897). (C) 4-Hydroxytamoxifen, the synthetic selective estrogen-receptor modulator (ChemSpider ID: 395987).

Selective estrogen receptor modulators (SERMs) such as 4-Hydroxytamoxifen (used to treat breast cancer, Figure 1C), Clomiphene (used to treat infertility), and Raloxifene (used to treat osteoporosis), which acts as competitive blockers of the estrogen/ZEA binding to the ER, have been used successfully (Y. Kim et al., 2020; Neven et al., 2018; Qin et al., 2020). Although these drugs are efficacious, they are associated with drug resistance and life altering side effects such as uterine cancer, endometrial cancer, menopausal symptoms, edema, pulmonary embolism and so on (Ignatov & Ortmann, 2020; Jeon et al., 2020; Zhang et al., 2016). Hence, the worldwide interest in the use of naturally occurring medicinal plant compounds such as flavonoids to circumvent debilitating side effects.

Flavonoids are phytochemicals synthesised as secondary metabolites in several parts of plants where they function as growth factors and a defence mechanism from infection and other stressors (Santos et al., 2017). They have been associated with a broad spectrum of health-promoting effects and are an indispensable component in a variety of nutraceutical, pharmaceutical, medicinal and cosmetic applications (Panche et al., 2016). They have antioxidative, anti-inflammatory, anti-mutagenic and anti-carcinogenic properties coupled with their capacity to modulate key cellular enzyme functions and are potent SERMs (Panche et al., 2016). Flavonoids are polyphenolic compounds consisting of 15-carbon skeleton: two phenyl rings (A and B) and a heterocyclic ring (ring C) (see Figure 2). They are classified into 6 subgroups (Flavones, anthocyanidins, flavanones, flavanols, isoflavones, flavonols) based on the degree of oxidation and hydroxylation of the ring structures (Panche et al., 2016). The hydroxylation often occurs at positions 3,5,7,3´,4´ and 5´(Figure 2) and often these hydroxyl groups are sulphated, methylated or acetylated (Santos et al., 2017). Flavonoids naturally occur as aglycones or glycosides with the sugar (D-glucose, L-rhamnose, glucose rhamnose, arabinose or galactose) moiety conjugated at either position 3 or 7. (Santos et al., 2017). Under natural conditions, most of the flavonoids are present in glycoside form. The modification to their structures alters their stability, reactivity and solubility (Santos et al., 2017). The hydrophobic and polar components on the rings interact with amino acid residues within the ligand binding domain of Estrogen receptors.

Figure 2. Basic chemical structures of several different flavonoid classes (Panche et al., 2016).

The estrogen receptor α (ERα) is a ligand activated transcription factor (Yaşar et al., 2017). The ERα has 6 functional domains: N-terminal (A/B) domain, the DNA-binding domain (C), the hinge domain (D) and the ligand binding domain (LBD, E/F).

The LBD is responsible for encompassing the ligand, it consists of a co-regulator interaction site (activation function 2, AF2), dimerisation interface and a binding pocket (Kumar et al., 2011). This domain has 12 alpha helices (H1-H12) and two beta sheets (S1 and S2) (see Figure 3A) (Brzozowski et al., 1997). The binding of a ligand to the LBD induces dynamic conformational changes in the helical structures of the binding pocket (BP), H12 conformational changes are mainly responsible for overall receptor activation (Puranik et al., 2019).

Figure 3. X-ray crystallography structures of the ERα ligand binding domain (LBD). (A) The LBD encompasses helices: H1-H12. The ligand binding pocket is comprised of H3, H6, H8, H11 and H2, PDB ID: 3ert. (B) ERα LBD complexed with zearalenone in the active/closed conformation, PDB ID: 5krc. (C) ERα LBD complexed with 4-hydroxytamoxifen in the inactive/open conformation, PDB ID: 3ert. (D) The apo form of ERα LBD complexed with 17β-estradiol, PDB ID: IA52.

When 17β-estradiol (E2) binds to the BP, it rests within the hydrophobic core of the BP formed by H3, H6, H8, H11 and H12 (see Figure 3B). The hydrophobic core of E2 forms non-polar interactions with the hydrophobic residues of the BP, the phenolic hydroxyl group from the A ring forms hydrogen bonds with Glu353 in H3 and Arg394 in H6 while the 17β-hydroxyl group of the D ring forms a hydrogen bond with His524 (H11) (Brzozowski et al., 1997). A similar picture is seen upon ZEA binding to the BP, the hydroxyl group on position 14 forms a hydrogen bond with Glu353 while the keto/hydroxyl group in 6^th position is bonded with Hie524 (Ehrlich et al., 2015; Zhang et al., 2019). The combination of the hydrophobic and polar interactions of ZEA and E2 within the BP results in an active/closed conformation. This shift in conformation causes the inner hydrophobic surface of helix H12 to be positioned across H3 and H11 thus forming a lid over the cavity to accommodate co-regulator binding (Figure 3B), which ultimately leads to enhanced gene transcription (Brzozowski et al., 1997; Shtaiwi et al., 2019).

Conversely, ERα competitive blockers such as 4-hydroxytamoxifen (OHT) and raloxifene occupy the same hydrophobic BP as ZEA and E2. The phenolic hydroxyl group of OHT interacts with Glu353 and Arg394, while the dimethylaminoethoxy group side chain lies near Asp351 and forms a direct hydrogen bond. The bulky side chain protrudes from the BP between H3 and H11 preventing the hydrophobic interactions of H12 with H3/H11 thus attaining an inactive/open conformation where H12 is positioned directly opposite H3 and H5, occupying the co-regulator binding site (Figure 3C), which reduces the level of estrogen-regulated gene transcription (Pike et al., 1999; Shtaiwi et al., 2019). However, the unbound form of ERα exhibits a unique conformation in which H12 does not share hydrophobic interactions with the BP instead it is directed away from the LBD core, as depicted in Figure 3d (Shtaiwi et al., 2019).

Over the last few decades, computer-aided drug design (CADD) strategic methods have greatly contributed toward the drug discovery, design and development that act against a wide diversity of bacteria, viruses and diseases (Yu & MacKerell, 2017). Computational techniques are frequently used to rapidly facilitate the discovery process and intentionally direct focus toward drug compounds that express potential pharmacological activities, provide enhancements for promising candidates and intricately details enzyme-drug interactions (Das et al., 2017; Kore et al., 2012).

Virtual screening (VS) is a process that utilises a series of computational tools to search through compound databases/libraries for suitable small molecules that have similar structuring to the template scaffold. Structures identified through VS should bind to the amino acid residues within the target’s active site (Appiah-Kubi & Soliman, 2017; Berinyuy & Soliman, 2017; N Mhlongo & ES Soliman, 2016). Structure-based VS searches for potential ligand candidates primarily based on the 3 D-structure of a target enzyme’s active binding pocket. On the other hand, ligand-based VS searches for favourable ligand features with known bioactive properties for pharmacophore model generation (Cele et al., 2016). Pharmacophore/ligand-based VS involve intermolecular interactions, which include hydrogen bonding, electrostatic energies and charge. This approach also screens for potential drug candidates that form interactions with highly contributing residues based on the binding free energy (Islam & Pillay, 2015).

Herein, we implemented molecular dynamics (MD) simulations, and post dynamic analysis to identify potential ERα inhibitors. This study provides the fundamental basis that will facilitate the selection of natural products with better suited antiestrogenic properties to ameliorate ZEA-induced pathologies.

Computational methodology

System preparation

The X-ray crystallographic structure of human ERα LBD complexed with OHT was obtained from the protein data bank (PDB code: 3ERT) (Burley et al., 2019). The complex was then processed using the UCSF Chimera software (Pettersen et al., 2004). The list of flavonoids was obtained from the Biological Magnetic Resonance Data Bank (BMRB) (Romero et al., 2020). We obtained 30 flavonoids (Table 1: supplementary materials) from the BMRB and docked all of them against the estrogen receptor. Furthermore, all the docked complexes with dock scores more than 8 were selected for further analysis (Molecular Mechanics/GB Surface Area method-MM/GBSA).

Table 1. Flavonoid compounds obtained from the Biological Magnetic Resonance Data Bank.

Compound	Abbreviation	Pubchem compound ID	Classification	Source	References
4-Hydroxytamoxifen (Afimoxifen)	OHT	449459	Selective estrogen receptor modifier	Synthetic hormone.	(Ahmad, 2018; Neven et al., 2018; Sanchez-Spitman et al., 2019)
Quercetrin (Quercetin 3-rhamnoside)	QRM	5280459	Flavonol	star anise, grapes, cocoa, pepper.	(Dabeek & Marra, 2019; Ulusoy & Sanlier, 2020)
Hesperidin (Hesperitin-7-rhamnoglucoside)	HPN	53477767	Flavanone	Citrus fruits	(Binkowska, 2020; Ganeshpurkar & Saluja, 2019; Selmi et al., 2017)
Epigallocatechin 3-gallate	EGCG	65064	Flavan-3-ols	Green tea extract	Vázquez Cisneros et al., 2017)
Kaempferol 7-O-glucoside	K7G	10095180	flavonol	Delphinium, Witch-hazel, grapefruit.	(Cid-Ortega & Monroy-Rivera, 2018; Huynh et al., 2018)
Naringenin-7-O-glucoside (Prunin)	NGN	92794	flavanone	Apigenin metabolite.	(Karim et al., 2018; Venkateswara Rao et al., 2017)
Apigenin 7-O-glucoside	APN	5385553	flavone	Chamomile.	(Cvetanović et al., 2017; Wang et al., 2019)
Daidzein	DZN	5281708	Isoflavone	Soybeans	(Sumardi et al., 2017)
Phloretin	PRT	4788	Chalcone	Fruits	(Behzad et al., 2017; Mariadoss et al., 2019)
Baicalein 6-glucoside	BCN	44258404	Flavone	Fruits, vegetables, medicinal plants	(de Oliveira et al., 2015; Dinda et al., 2017)
Luteolin	LTN	5280445	flavone	Flowering food plants.	(Somerset & Johannot, 2008)
Fisetin	FST	5281614	flavonol	Berries, grapes.	(Grynkiewicz & Demchuk, 2019; Khan et al., 2013; Somerset & Johannot, 2008)
Myricetin	MCT	5281672	flavonol	Vegetables. fruits and beverages.	(Gupta et al., 2020)
Tangeretin	TGN	68077	Flavone	Apples, tangerines, Fortunella, japonica	(Ashrafizadeh et al., 2020)
Cyanidin	CYN	128861	Anthocyanidin	Fruits and vegetables.	(Somerset & Johannot, 2008)

All 2-D (Figure 4) and 3-D ligand structures were obtained from PubChem (S. Kim et al., 2019) and processed using the chimera software (Pettersen et al., 2004). All 14 ligands including OHT, which was used as a ERα inhibitor reference, were docked to the Erα-LBD.

Figure 4. 2 D structures of the flavonoids.

Molecular docking

The compounds were docked into the LBD of ERα. The docking protocol was validated by redocking the original substrate OHT to the processed LBD and superimposing the redocked complex to the original crystal structure. The RMSD obtained was 0.004 Å (Figure 5). In preparation for docking, non-polar hydrogen atoms were merged to carbon atoms and Gasteiger charges were added to all ligands docked into the binding pocket of the ERα LBD through UCSF Chimera (Pettersen et al., 2004), using the AutoDock Vina plugin (Pettersen et al., 2004), compounds were docked by means of the grid box function to express the binding site (spacing of 1Å was marked out at a size of 22.3554 × 22.7154 × 26.9884 in the x, y and z three-dimensional axis respectively. All the docked complexes were then subjected to molecular dynamic simulations.

Figure 5. Lateral view of ERα LBD showing the validation of the docking protocol. The original crystal ERα LBD structure (PDB: 3ert) in green was superimposed with the newly docked complex (in blue and OHT in brown). RMSD = 0.004A˚.

Molecular dynamic (MD) simulations

Molecular dynamics (MD) simulation can provide not only plentiful dynamical structural information on biomacromolecules but also a wealth of information concerning the energetics of protein and ligand interactions. Such information is very important to understanding the structure-function relationship of the target molecules, the essence of protein–ligand interactions and to guiding the drug discovery and design process. Thus, MD simulations have been applied widely and successfully in each step of modern drug discovery (Liu et al., 2018).

The MD simulations were performed using AMBER18 with GPU acceleration (Hollingsworth & Dror, 2018). The enzyme was dehydrogenated and parameterized using the AMBER FF18SB force field (Hollingsworth & Dror, 2018). Hydrogen atoms were added to the ligands and charged via Gasteiger charges. The antechamber module was used to generate partial atomic charges for the ligands by applying General Amber Force Field (GAFF) and restrained electrostatic potential (RESP) procedures (Hollingsworth & Dror, 2018). The LEAP module was applied to neutralise and solvate the systems (Le Grand et al., 2013) by adding hydrogen, sodium and chloride atoms, whilst suspending them in an orthorhombic box of TIP3P water molecules so that all atoms were within 8 Å of the box edges. The amino acid residues of the LBD were numbered 305-554.A 2000 step partial minimisation was performed, while applying a 1000 step restraint potential of 500 Kcal/mol, using the steepest descent method followed by 1000 steps of conjugate gradients. Furthermore, 1000 steps of full minimisation were executed using conjugate gradient algorithm without restraint. A heating MD simulation from 0 K to 300 K was gradually carried out for 50 ps, maintaining a fixed number of atoms and volume. A potential harmonic restraint of 10 kcal/mol and collision frequency of 1.0ps were applied on to all the solutes in the system and Langevin thermostat with collision frequency of 1.0 ps⁻¹ was applied to each system. After heating, about 500ps of equilibration was executed while keeping a constant operating temperature at 300 K. Features such as several atoms and pressure were also kept constant mimicking an isobaric-isothermal ensemble. The Berendsen barostat was used to preserve the system pressure at 1 bar and the SHAKE algorithm was utilized for hydrogen bond constraint (Gonnet, 2007; Le Grand et al., 2013). SPFP (Single precision with fixed point integer arithmatic) precision model was used and the step size of each simulation was 2fs (Le Grand et al., 2013).

Post simulation analysis

The coordinates of the systems were saved, and the trajectories were analyzed for every ps using CPPTRAJ and PTRAJ modules in the Amber18 package (Shunmugam & Soliman, 2018). The following post analysis were conducted root mean square deviation (RMSD), root means square fluctuation (RMSF) and radius of gyration (RoG). The Origin data analysis software (Seifert, 2014) was used to plot data and generate graphs.

RMSD

The root mean square deviation (RMSD) is a similarity measure widely used to estimate the deviation in the Cα (alpha carbon)-atoms of protein residues indicative of the degree of protein stability (Dixit et al., 2006). The Cα RMSD plots of the apo and bound systems were analysed to assess the effect of flavonoid binding on the structural stability of the ERα-LBD.

RoG

Radius of gyration (RoG) is defined as the mass weighted root mean square distance of a collection of atoms from their common centre of mass. It is an indicator of protein structural compactness (Baig et al., 2014). This method of analysis allowed the assessment of ERα-LBD compactness upon flavonoid binding along the 200 ns trajectory in the MD simulation.

RMSF

Root means square fluctuation (RMSF) is a measure of the displacement of a particular atom, or group of atoms, relative to the reference structure, averaged over the number of atoms. It is frequently used to discern whether a structure is stable in the timescale of the simulations or if it is diverging from the initial coordinates (Martínez, 2015). The RMSF analysis provides insights into the flexibility of the ERα-LBD amino acid residues upon binding to the flavonoids.

Binding free energy calculations

The MD allows free energy differences to be accepted, which regulate the underlying mechanism of all biological processes (Gapsys et al., 2015).Binding free energy (BFE) calculations are a key tool for studying the in-depth binding mechanism between a protein and ligand, which involves both entropic and enthalpic contributions (Schauperl et al., 2017). The BFE is used in combination with the Molecular Mechanics/Generalized Born Surface Area (MM/GBSA) technique to estimate the binding affinity within a docked system. In contrast to a full-scale MD free energy perturbation/thermodynamic simulation, this approach can be used to provide reproducible relative binding affinities of compounds with good precision and needs significantly less computational resources (Shunmugam & Soliman, 2018).

The free binding energy was calculated using the Molecular Mechanics/GB Surface Area method (MM/GBSA) to estimate and compare the binding affinity of the systems. This method was employed as previously described (Idowu et al., 2020). Binding free energy was averaged over 200 000 snapshots extracted from the 200 ns trajectory. The free binding energy (ΔG) computed by this method for each molecular species (complex, ligand, and receptor) can be represented as: (1) $$\Delta G_{\mathit{bind}}=\mathrm{ }\Delta G_{\mathit{complex}}-\mathrm{ }\Delta G_{\mathit{receptor}}-\mathrm{ }\Delta G_{\mathit{ligand}}$$(1) (2) $$\Delta G_{\mathit{bind}}=\mathrm{ }{E}_{\mathit{gas}}+\mathrm{ }{G}_{\mathit{sol}}–\mathrm{ }T\Delta S$$(2) (3) $$E_{\mathit{gas}}=\mathrm{ }{E}_{\mathit{int}}+\mathrm{ }{E}_{\mathit{vdw}}+\mathrm{ }{E}_{\mathit{ele}}$$(3) (4) $$G_{\mathit{sol}}+G_{GB}+\mathrm{ }{G}_{SA}$$(4) (5) $$G_{SA}+\mathrm{ }\gamma \mathit{SASA}$$(5)

The term E_gas denotes the gas-phase energy, which consists of the internal energy E_int; Coulomb energy E_ele and the van der Waals energies E_vdw. The E_gas was directly estimated from the FF14SB force field terms. Solvation free energy, G_sol, was estimated from the energy contribution from the polar states, GGB, and non-polar states, G. The non-polar solvation energy, SA. GSA, was determined from the solvent-accessible surface area (SASA), using a water probe radius of 1.4 Å. In contrast, the polar solvation, GGB, the contribution was estimated by solving the GB equation. S and T denote the total entropy of the solute and temperature, respectively.

Results and discussion

Thermodynamic binding free energy of the flavonoids to ERα

To determine the inhibitory potential of the flavonoids, the binding affinities of the flavonoids to the ERα was estimated by calculating their binding free energy (ΔG_bind) using the MM/GBSA method. This method was previously used to calculate ΔG_bind of ligands to receptors in several computational studies (Idowu et al., 2020; Mazurek et al., 2020; Shunmugam & Soliman, 2018). Following MD simulation of 200 ns, the binding energies of all the bound systems were calculated (Table 2).

Table 2. Thermodynamic binding free energy components (Kcal/mol) for all the flavonoids compared to OHT (reference drug).

	Energy components (kcal.mol^-1)
	ΔEvdw	ΔEelec	ΔGgas	ΔGsolv	ΔGbind
OHT	−50.744 ± 3.170	−170.383 ± 24.083	−221.127 ± 24.936	176.402 ± 22.526	−44.724 ± 4.306
HPN	−69.495 ± 4.825	−41.980 ± 15.293	−111.476 ± 16.009	55.410 ± 11.203	−56.066 ± 6.696
QRM	−40.415 ± 5.505	−74.587 ± 22.181	−115.002 ± 18.785	59.369 ± 10.954	−55.633 ± 9.099
EGCG	−39.509 ± 5.004	−77.166 ± 11.204	−116.675 ± 11.998	68.5883 ± 8.165	−48.087 ± 6.098
K7G	−45.012 ± 3.806	−54.586 ± 12.723	−99.598 ± 11.091	53.089 ± 8.014	−46.509 ± 4.428
NGN	−37.498 ± 2.903	−30.468 ± 4.055	−67.967 ± 3.623	25.023 ± 2.066	−42.943 ± 2.688
APN	−39.777 ± 5.352	−52.169 ± 11.325	−91.947 ± 12.127	49.196 ± 8.647	−42.750 ± 5.105
PRT	−39.837 ± 2.908	−29.380 ± 6.617	−69.217 ± 6.343	27.201 ± 3.352	−42.016 ± 4.224
DZN	−35.808 ± 2.811	−33.620 ± 3.661	−69.429 ± 3.213	27.440 ± 1.868	−41.989 ± 2.514
BCN	−50.919 ± 3.978	−19.380 ± 9.151	−70.300 ± 9.561	27.861 ± 6.837	−42.438 ± 4.751
LTN	−29.569 ± 3.948	−52.599 ± 8.135	−82.168 ± 7.714	41.990 ± 5.680	−40.178 ± 4.209
FST	−39.759 ± 3.657	−28.113 ± 7.935	−67.873 ± 6.668	28.618 ± 4.392	−39.254 ± 3.568
MCT	−30.390 ± 3.710	−57.558 ± 8.732	−87.948 ± 7.334	51.083 ± 6.281	−36.865 ± 3.158
TGN	−42.146 ± 6.023	−18.402 ± 7.392	−60.549 ± 10.457	25.905 ± 6.648	−34.644 ± 5.854
CYN	−38.289 ± 5.012	−156.49 ± 12.523	−194.78 ± 12.216	160.238 ± 9.477	−34.548 ± 4.567

ΔEelec: electrostatic, ΔEvdW: van der Waals, ΔGbind: calculated total free binding energy, ΔGgas: gas phase interaction, and ΔGsolv: solvation free energy.

The binding affinities of flavonoid bound systems were compared to the OHT complex. We postulated that compounds with the lowest ΔG_bind are considered to have better binding affinity and enhanced inhibitory potential. HPN (-56.07 Kcal/mol), QRM (-55.63 Kcal/mol), EGCG (-48.09 Kcal/mol) and K7G (-46.51 Kcal/mol) had lower ΔG_bind compared to OHT (-44.72 Kcal/mol). This finding may imply that these specific compounds express a much more favourable binding affinity for ERα than the reference synthetic drug, OHT, ultimately suggesting that HPN, QRM, EGCG and K7G are potentially better ERα Inhibitors.

Majority of previous studies based on the binding of flavonoids use docking scores to predict the binding affinity of flavonoids to ERα (Indira et al., 2020; Maruthanila et al., 2019; Pang et al., 2018; Powers & Setzer, 2015; Puranik et al., 2019; Widowati et al., 2019). To the best of our knowledge, this is one of a few studies to use the MM/GBSA method to estimate the binding affinity of flavonoids to ERα. As a consequence, to HPN, QRM, EGCG and K7G superior binding affinity, the compounds were selected for further structural and dynamical analysis.

ERα- ligand complex stability

The RMSD readings of alpha carbon (Cα) atoms throughout the 200 ns MD simulation of the apo and bound systems were analysed to assess the effect of the binding mechanisms of QRM, HPN, K7G and EGCG on the structural stability of the ER (Figure 6). The RMSD analysis serves as a measure of system convergence and stability, which influences overall receptor function (Idowu et al., 2020). These measurements are determined by ligand-receptor interactions that arise from ligand movements in the binding pocket (BP) to achieve optimum geometric orientation (Mazurek et al., 2020).

Figure 6. (A) The RMSD plots of Cα for the apo (black), OHT (red), HPN (cyan) and QRM (purple). (B) The RMSD plots of Cα apo (black), OHT (red), K7G (green) and QRM (blue).

The RMSD plots showed that OHT, HPN and QRM complexes each converged at approximately 20 ns (Figure 4a), while the apo reached stability at 50 ns and the remaining systems converged at 60 ns (EGCG and K7G) (Figure 4b).

The findings of this study are similar to those observed in previous studies conducted by Celik et al. (2007) and Shtaiwai et al. (2019) in which the Cα atom RMSD of a OHT complex converged at 20 ns and the apo system showed greater Cα atom deviations than the bound systems.

The results suggest that the binding of the ligands to the receptor stabilised the structural complexes following an initial period of high fluctuations in order to achieve the optimum geometric conformation (Mazurek et al., 2020).

To gain more insight into the stability and structural conformation of the flavonoid complexes, the radii of gyration (RoG) were analysed and graphically represented in Figure 7. RoG is a measure of protein folding, stability and structural compactness. According to relevant figure, it is seen that the RoG value for the apo and OHT, HPN, EGCG, QRM and K7G-bound systems are reasonably similar, with an average value of 18.8Å, 18.6Å, 18.5Å,18.6Å, 18.7Å and 18.7Å respectively. These findings suggest that compound binding to the LBD of ERα induced a conformational shift to a slightly more compact structure, despite there only being a minimal change in distance, [OHT (0.2Å), HPN (0.3Å), EGCG (0.2Å), QRM (0.1Å) and K7G (0.1Å)], resulting in a steady and lower energy system. This will ultimately form a relatively stable protein (Shunmugam and Soliman. 2018).

Figure 7. Radii of gyration (RoG)of apo and bound systems were measured over a 200 ns simulation. Data obtained were presented as a plot exhibiting the differences arising in radius deviation between the apo and bound systems. (A) The RoG plots for apo(black), OHT (red), HPN (cyan) and QRM (purple). (B) The RoG plots for apo (black), OHT (red), K7G (green) and QRM (blue).

ERα- ligand complex residue fluctuations

The RMSF values of the apo and bound systems were analysed to determine the binding effects of OHT, HPN, EGCG, QRM and K7G on the dynamic residue behaviour within the ERα-LBD. The RMSF with respect to the averaged MD simulation conformation was used to demonstrate the differences occurring in residue flexibility.

The RMSF of all residues were calculated for flexibility of the structure backbone and graphically depicted in Figure 8. The RMSF values of HPN, EGCG, QRM and K7G were compared to OHT.

Figure 8. (A) The RMSF plots for the apo system (black), OHT complex (red), HPN (cyan) complex and QRM (purple) complex. (B) The RMSF plots for the apo system (black), OHT (red) complex, K7G (green) complex and QRM (blue) complex.

All systems, free and bound, exhibited increased dynamic fluctuations of the amino acid residues (R) in four specific regions corresponding to: loops preceding H3 (R333-R341), H9 (R454-R471) and those surrounding the terminal helix H12 (R525-R535 and R545-R551). This finding is congruent with previous studies in literature (Fratev, 2015; McGee et al., 2008; Ng, 2016; Shtaiwi et al., 2019)

The most evident amino acid fluctuations observed in the apo system was identified as residues that surround H12 (Table 3). These findings are consistent with the RMSD and RoG data obtained. The flexibility of the Cα atoms as reflected in the RMSD plots and the loose structural configuration represented in the RoG graphs is attributed to the amino acid residue movements in these regions (Figure 8) (Ng, 2016). The residue variations in these regions could be due to the lack of ligand interactions, which lock to keep the helices in place for appropriate receptor function (Ye et al., 2018). This occurrence was also observed in a previous study whereby removal of the ligand from the ERα BP led to increased residue fluctuations specifically in the loop region located at the head of H12 (McGee et al., 2008).

Table 3. The average Cα RMSF values of specific regions of the ER LBC.

Helices	Residues	Apo	OHT	HPN	QRM	EGCG	K7G
H3	R333-R341	2.8	2.6	4.4	3.8	5.4	4.7
H9	R454-R471	2.2	3.3	3.6	2.6	4.7	3.3
H12	R525-R535	4.8	4.2	2.8	3.2	3.6	3.6
H12	R545-R551	8.9	6.0	2.7	4.8	3.3	6.0

The units of measurements are in Angstrom.

Similarly, the OHT, QRM and K7G bound system residues exhibited their highest movements in the terminal H12 regions while that of the EGCG complex were lowest in these regions. The flexibility seen in the OHT bound system is attributed to the conformation of H12 attained when the bulky side chain for OHT protrudes from the BP between H3 and H11, thus preventing hydrophobic interactions of H12 with the core of the binding pocket (Brzozowski et al., 1997). All the polar and hydrophobic interactions in the QRM and K7G complexes were localised in H3 and H5 (Figure 9). The lack of interactions with the residues of H11 and H12 accounts for increased movements in the loops surrounding these helices because they are devoid of the ligand interactions which would restrict their flexibility (Ye et al., 2018). However, EGCG formed 6 out of its 8 hydrogen bonds with residues of H11 and H12. The presence of the EGCG hydrogen bonds pulls the helices inwards, and locks them in place, hence reducing their movement.

Figure 9. The graphical representation of HPN, QRM, EGCG and K7G residue interactions within the LBD of ERα.

Contrastingly, the residue fluctuations of the HPN complex were highest in the H3 and H9 regions (Table 3). Six out of the seven hydrogen bonds in the HPN complex are in H3 and H6 (Table 4). The presence of these interactions are likely to reposition the helical residues towards HPN (Ye et al., 2018) and this indicates that during this phenomenon, the loops that surround the affected helices are more than likely to exhibit increased dynamic fluctuations (Brzozowski et al., 1997).

Table 4. The flavonoid- ERα interaction details. The number of interactions and the interacting residues.

Compound	Total interactions (Hydrogen bonds)	Types of hydrogen bonds	Interacting residues (hydrogen bonds)
OHT	27 (9)	Conventional:2 Carbon hydrogen: 7	HIE524, CYS530. THR 347, THR347, ASP351, GLY521, CYS530, CYS530, ASN532.
HPN	26 (7)	Conventional:2 Carbon hydrogen:5	GLU339, THR347 GLU339, LEU346, THR347, TRP383, CYS530
QRM	17 (9)	Conventional: 6 Carbon hydrogen:3	GLU323, GLU353, GLU353, HIE356, ARG394, HIE356 Arg394, PHE445, HIE356
EGCG	17 (8)	Conventional: 7 Carbon hydrogen: 1	ASP351, ASP351, MET522, ASN532, ASN532, VAL534, VAL534 TYR526
K7G	18 (6)	Conventional: 4 Carbon hydrogen: 2	LEU327, GLU353, GLU353, ILE386 GLY390, LYS444

ERα and ligand interactions

The structural properties of flavonoids determine the nature of their interactions with the amino acid residues of the Erα-LBD pocket. To assess the binding mechanisms of QRM, HPN, K7G and EGCG with the receptor’s amino acid residues, snapshots of each compound were captured at 200 ns of the MD simulation and subjected to analysis (Figure 9).

As observed in Table 4, HPN had the highest number of interactions (26) and carbon hydrogen: conventional hydrogen bond ratio (5:2), which is comparable to OHT with 26 interactions and a 7:2 carbon hydrogen: conventional hydrogen bond ratio. It is well established in literature that HPN has a double sugar side chain with multiple hydroxyl groups, which are responsible for most of the hydrogen bonds (Binkowska, 2020; Kuntić et al., 2014). QRM and EGCG have 17 interactions each and carbon hydrogen: conventional hydrogen bond ratio of 3:6 and 1:7, respectively. K7G has the lowest overall number of hydrogen bonds. QRM and K7G are glycosides while EGCG is aglycone (Panche et al., 2016). However, the most reactive hydroxyl groups (7-OH in flavones or the 3-OH in flavonols) in flavonoids are glycosylated (Panche et al., 2016). Hence, this explains the lower overall number of interactions, number of carbon hydrogen bonds and the binding affinities observed (Table 2).

The amino acids residues Glu353, Arg394, and Hie524 form the binding mode of the endogenous ERα stimulator 17β-estradiol (E2) and are therefore considered as estrogenic residues as they are important in the recognition of estrogenic ligands (Brzozowski et al., 1997; Shtaiwi et al., 2019). However, the estrogenic residues in addition to Asp351 (Brzozowski et al., 1997; Celik et al., 2007; Shtaiwi et al., 2019) Thr347 and Leu346 form the binding modes of most SERMs (such as OHT and raloxifene) and they are therefore considered as antiestrogenic residues which are important in recognising antiestrogenic ligands (ERα inhibitors) (Li et al., 2013).

The alkyl side chain of OHT forms hydrogen bonds with Thr347, Asp351, Gly521, Hie524, Cys530 and Asn532. Notably, OHT formed strong C-H covalent bonds with the antiestrogenic residues Thr347 and Asp351.The hydroxyl groups on the rhamnoglucoside chain of HPN formed hydrogen bonds with Glu339, Leu346, Thr347, Trp383, and Cys530; and a covalent bond with Thr347 and Leu346. These results are consistent with previous docking studies that reported OHT and HPN binding to Asp351 and Thr347 (Puranik et al., 2019). The glucoside chain of QRM formed direct hydrogen bonds with Glu323, Glu353, Hie356, Arg394, Phe445 similar to Puranik’s (2019) team who reported the binding of the aglycon version of QRM to Glu353 and Arg394. On the other hand, EGCG and K7G formed hydrogen bonds with Asp351, Met522, Asn532, Val534, Tyr526 and Leu327, Glu353, Ile386, Gly390, Lys444, respectively. EGCG also shared a conventional hydrogen bond with Asp351. Contrary to the findings of this study, Pang et al. (2018) reported a lack of binding interactions between EGCG and Asp351

Interestingly, QRM and K7G formed hydrophobic interactions with the same amino acid residues. These interactions may be due to similarities arising in the unconjugated skeleton structures of both compounds. The only difference observed is the OH group on C-3 of the B- ring in QRM, which is replaced with H in K7G (Dabeek & Marra, 2019). However, the sugar sidechains of both compounds formed double conventional hydrogen bonds with Glu353, but none were observed to be formed with the antiestrogenic residues. These findings contradict those expressed by Puranik et al. (2019) in which K7G and QRM interacted with Asp351 and Thr347, respectively (Puranik et al., 2019). This difference may be accredited to the use of aglycon versions of K7G and QRM in the present study. (Martínez, 2015)

Conclusions

Zearalenone is an estrogenic mycotoxin which is a common human and animal food contaminant. It has been implicated in carcinogenesis and other reproductive ailments through the estrogen receptor (ER) pathway. FDA approved synthetic selective estrogen receptor modulators (SERM) such as tamoxifen (OHT) which are used to treat these ailments induce life altering side effects albeit being effective. Hence the need for natural occurring ER inhibitors as medicinal alternatives such as flavonoids to circumvent side-effects. This study aimed to identify flavonoids ERα inhibitors.

Inhibitors of ERα extrapolated from a library of flavonoids were identified and the binding mechanism of the inhibition were investigated using molecular mechanics/GB surface area method (MM/GBSA) calculations and molecular dynamic simulations. The MM/GBSA binding free energy calculations revealed that HPN QRM, K7G and EGCG had a higher binding affinity for the receptor than reference drug, OHT. The ligand protein interactions revealed that HPN had the greatest total number of interactions and carbon-hydrogen bonds while QRM and K7G interacted with the same amino acid residues. Collectively the results revealed that the binding of the compounds to the receptor lead to dynamic alterations, which induced conformational shift in the structure and orientation of the receptor resulting in stabilised, compact and low energy systems.

Cumulatively, the results of this study support the use of flavonoids, especially HPN and QRM, as a naturally sourced alternative therapeutic agents against the pathological effects of zearalenone exposure in humans.

Abbreviations
ΔGbind	=	Binding free energy
BP	=	Binding pocket
Cα	=	Alpha carbon
EGCG	=	Epigallocatechin 3-gallate
ERα	=	Estrogen receptor alpha
HPN	=	Hesperidin
K7G	=	Kaempferol 7-O-glucoside
LBD	=	Ligand binding domain
MD	=	Molecular dynamics
MM/GBSA	=	Molecular Mechanics/Generalized Born Surface Area
OHT	=	4-Hydroxytamoxifen
QRM	=	Quercetrin
RMSD	=	Root mean square deviation
RMSF	=	Root mean square fluctuation
RoG	=	Radius of gyration
SERMs	=	Selective estrogen receptor modulators
ZEA	=	Zearalenone

Acknowledgements

The authors wish to acknowledge the College of Health Science at the University of KwaZulu-Natal and the National Research Foundation (NRF) for the financial support and The Centre of High-Performance Computing (CHPC, https://www.chpc.ac.za), Cape Town, RSA, for computational resources.

Disclosure statement

The authors declare no potential conflicts of interest.

Funding

This work did not receive direct funding however the student and supervisor involved received some financial support from the parties mentioned.