Next Article in Journal
Network Pharmacology Study on Morus alba L. Leaves: Pivotal Functions of Bioactives on RAS Signaling Pathway and Its Associated Target Proteins against Gout
Next Article in Special Issue
Treatment of HT29 Human Colorectal Cancer Cell Line with Nanocarrier-Encapsulated Camptothecin Reveals Histone Modifier Genes in the Wnt Signaling Pathway as Important Molecular Cues for Colon Cancer Targeting
Previous Article in Journal
Auxin Metabolome Profiling in the Arabidopsis Endoplasmic Reticulum Using an Optimised Organelle Isolation Protocol
Previous Article in Special Issue
Structure-Activity Relationships of the Imidazolium Compounds as Antibacterials of Staphylococcus aureus and Pseudomonas aeruginosa
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Elucidation of Agonist and Antagonist Dynamic Binding Patterns in ER-α by Integration of Molecular Docking, Molecular Dynamics Simulations and Quantum Mechanical Calculations

Division of Bioinformatics and Biostatistics, National Center for Toxicological Research, U.S. Food and Drug Administration, 3900 NCTR Road, Jefferson, AR 72079, USA
*
Author to whom correspondence should be addressed.
Current affiliation: Department of Bioinformatics, Alagappa University, Karaikudi 630 003, Tamil Nadu, India.
Int. J. Mol. Sci. 2021, 22(17), 9371; https://doi.org/10.3390/ijms22179371
Submission received: 27 May 2021 / Revised: 25 August 2021 / Accepted: 27 August 2021 / Published: 29 August 2021
(This article belongs to the Special Issue Application of In Silico Techniques in Drug Design)

Abstract

:
Estrogen receptor alpha (ERα) is a ligand-dependent transcriptional factor in the nuclear receptor superfamily. Many structures of ERα bound with agonists and antagonists have been determined. However, the dynamic binding patterns of agonists and antagonists in the binding site of ERα remains unclear. Therefore, we performed molecular docking, molecular dynamics (MD) simulations, and quantum mechanical calculations to elucidate agonist and antagonist dynamic binding patterns in ERα. 17β-estradiol (E2) and 4-hydroxytamoxifen (OHT) were docked in the ligand binding pockets of the agonist and antagonist bound ERα. The best complex conformations from molecular docking were subjected to 100 nanosecond MD simulations. Hierarchical clustering was conducted to group the structures in the trajectory from MD simulations. The representative structure from each cluster was selected to calculate the binding interaction energy value for elucidation of the dynamic binding patterns of agonists and antagonists in the binding site of ERα. The binding interaction energy analysis revealed that OHT binds ERα more tightly in the antagonist conformer, while E2 prefers the agonist conformer. The results may help identify ERα antagonists as drug candidates and facilitate risk assessment of chemicals through ER-mediated responses.

1. Introduction

Estrogen receptor (ER) is one of the important targets of drugs and endocrine disrupting chemicals in the endocrine system [1]. It is a ligand-dependent transcriptional factor in the steroid type 1 nuclear receptor family [2]. ER plays a major role in various biological functions such as bone modeling, reproductive system, cardiovascular system, metabolism, and cell proliferation [3]. ER is an extensively studied target among the endocrine receptors. There are two major ER isoforms, ERα and ERβ. Like other nuclear receptors, ERα consists of three distinct domains: N-terminal domain (residue 1–180), DNA binding domain (residue 181–263), and C-terminal domain or ligand binding domain (LBD, residue 303–552) (Figure 1). The activation function domain 1 (AF1) is present in the N-terminal domain and plays a major role in the protein–protein interaction [4,5]. The mitogen-activated protein (MAP) kinase pathway regulates the activity of AF1 through the growth factors [6]. The LBD is composed of twelve helices and two antiparallel β-sheets which are arranged as a three-layer antiparallel α helical sandwich [5,7]. The first layer is formed by helices 1 to 4 and 7, the middle layer is made up of helices 5, 6, 9 and 10 and the final layer is composed of helices 8, and 11 [8,9,10]. The activation function domain 2 (AF2) in LBD is responsible for binding of cofactors. AF2 undergoes conformational change due to the binding of a compound in the ligand binding pocket (LBP). The conformational change of AF2 determines the types of binding cofactors which play a major role in activating or inhibiting the target genes of ER [11]. The hinge (residue 264–302) regions connect the DBD and LBD.
The H12 acts as a molecular switch that turns ER activity on and off depending on the binding chemicals [11,12,13]. The AF1 and AF2 play a major role in the transcriptional activation of ER [14]. The estrogenic compounds bind in the hydrophobic pocket of ER LBD. The hydrophobic pocket is composed of Met342 to Leu354 of H3, Trp383 to Arg394 of H6, Val418 to Leu428 from the preceding loop of H8, Met517 to Met528 of H11, Leu539 to His547 of H12, and Leu402 to Leu 410 of S1/S2 hairpin [7]. Binding of antiestrogenic compounds to ER induce a H12 conformational changes by placing H12 across the H3 and H11 and moving H12 away from the LBP. Due to the H12 conformation change, the AF2 in the LBD is distorted and not suitable for binding cofactors [15]. ER enhances and represses its function via various pathways [10,16,17]. Understanding the ERα dynamic binding patterns with agonists and antagonists is crucial for discovery of ERα agonists and antagonists. Dynamic binding pattern represents the forming and breaking of non-covalent interactions such as hydrogen bonding and Van der Waals interactions between a protein and a ligand, as well as conformational changes caused by the binding ligand throughout a molecular dynamics (MD) simulation. More than 350 3D structures of ER bound with various ligands are deposited in the Protein Data Bank (PDB). Those structures are useful to understand the structural changes due to agonist and antagonist binding in the ERα LBP. Various computational techniques such as molecular docking [18,19,20,21,22,23,24], MD simulations [25,26,27,28,29,30], predictive modeling [31,32,33,34,35,36,37,38,39,40,41], and in vitro studies were conducted to predict ER binders or non-binders [42,43] and agonists or antagonists [44,45].
Many ligand-based computational methods were used to predict ERα activity of chemicals based on chemical features, including ERα binders and nonbinders, and agonist and antagonist activities [31,46,47,48,49]. However, the dynamic binding patterns of ER agonists and antagonists are not clearly understood. Hence, in this study we applied QM-Polarized Ligand Docking (QPLD) and MD simulations to elucidate the dynamic binding patterns of ERα agonists and antagonists using 17β-estradiol (E2) and 4-hydroxytamoxifen (OHT). E2 is the natural steroid hormone that activates ER. The activated ER modulates gene expression in cells. E2 binds with ER in the nucleus and forms a dimer. Subsequently, the dimer interacts with the estrogen response element of ER and regulates transcription of the target gene. E2 has two hydroxyl groups, one at C3 and another at 17β (Figure 2). The hydroxyl group at C3 forms hydrogen bonds with Glu353 and Arg394 of ER. The 17β-OH group forms a hydrogen bond with His524. The planar part of A/B ring forms a sandwich between Ala350 and Leu387. The D ring forms a nonpolar contact with Ile424, Gly521 and Leu525 [7]. OHT is a selective estrogen receptor modulator and acts as an antagonist towards ER in specific tissues [6]. The hydroxyl group in OHT (Figure 2) has a high binding affinity towards ER [50]. Binding of OHT in the LBP of ER pushes the H12 away to occupy part of the AF2 site, blocking coactivators binding in AF2 [51]. In QPLD, the ab initio molecular charges were applied to obtain binding orientation of the two chemicals in the LBP of ERα [52]. MD simulations were used to elucidate the dynamic binding patterns of the agonist E2 and antagonist OHT.

2. Results

2.1. Molecular Docking

QPLD and Glide docking are two widely used docking methods to identify orientations of compounds in binding sites of proteins [52,53,54]. The EXtra-Precision (XP) Glide scores, QPLD scores and docking energy values of the four complexes with the best ligand orientations are shown in Table 1.
In both docking methods, the agonist E2 had a lower docking energy in the agonist conformation (ERα1) than in the antagonist conformation (ERα2), while the antagonist OHT had a higher docking energy in the agonist conformation (ERα1) than in the antagonist conformation (ERα2). The orientations of E2 and OHT in ERα1 and ERα2, as well as the ERα residues interacting with E2 and OHT in the four complexes, are depicted in Figure 3. E2 in the binding site of ERα1 forms interactions with Glu353, Arg394, and His524. E2 also forms interactions with Glu353 and Arg394 but fails to interact with His524 in the binding site of ERα2. OHT forms hydrogen bond interactions with both Arg344 and Glu353 in the binding site of ERα2. OHT interacts with Arg344 but fails to interact with Glu353 in ERα1. The interaction analysis revealed that E2 forms more hydrogen bond interactions in ERα1, while OHT forms more hydrogen bond interactions in ERα2. These four complexes were subjected to MD simulations to elucidate the dynamic binding patterns of agonist E2 and antagonist OHT in ERα.

2.2. MD Simulations

In the MD simulations for each of the four complex structures obtained from molecular docking, structures were recorded for every 4.8 picoseconds (ps) in the trajectory file. Thus, each trajectory file contains 20,835 structures. The details of the simulation systems are summarized in Table 2.
To understand the dynamics of ER binding with the agonist E2 and antagonist OHT, root mean square deviations (RMSD) were calculated between the 20,835 structures for ERα1_E2, ERα1_OHT, ERα2_E2, and ERα2_OHT using a MATLAB script. The obtained RMSD matrixes for the four complexes are shown in Figure 4. Examining the RMSD values from the MD simulations (shown in Figure 4) found that structural changes were not the same during the MD simulations and the RMSD matrixes formed patterns. Furthermore, the RMSD patterns are different among the four complexes, indicating the binding dynamics of the agonist and antagonist in ERα are different.
The interaction energy was calculated using prime MM/GBSA for whole trajectory files and the calculated energy values are provided in Table 3.
To identify distinct structural patterns, hierarchical clustering analysis was conducted based on the RMSD matrixes. The major clusters (with >500 structures) from the MD simulations for ERα1_E2, ERα1_OHT, ERα2_E2, and ERα2_OHT are summarized in Table 4.
To further examine conformational changes of ER caused by agonist (E2) and antagonist (OHT) binding in the simulations, we used the distance between H12 and the centroid of LBP of ER to measure the conformational changes of H12. The residues (only in the secondary structure, not from loop region) around 4 Å to the ligands in 1GWR and 3ERT were selected for calculating the centroid of LBP: Met343, Leu346, Thr347, Leu349, Ala350, Asp351, Glu353, Leu384, Leu387, Leu391, Arg394, Phe404, Met421, Leu428, Gly521, His524 and Leu525. The residues Asp535 to Leu549 were selected to represent H12. First, a distance was calculated between the centroid and each of the H12 residues for a structure in the simulations. The maximum of the H12 distances was used to measure the distance of H12 from the centroid. For each structure in the simulations, a relative distance was calculated by subtracting the distance of the initial structure. The resulting relative distances for simulations of ERα2_E2 and ERα1_OHT are shown in Figure 5. Most of the structures in the simulation of ERα1-OHT had a longer distance between H12 and the centroid of LBP than the initial structure (positive relative distances, top of Figure 5), indicating the antagonist OHT made conformational changes in H12 from the initial active form towards to inactive form. On the other hand, most of the frames in the trajectory file from simulation of ERα2_E2 had a shorter distance between H12 and the centroid of LBP than the initial structure (negative relative distances, bottom of Figure 5), indicating the agonist E2 caused conformational changes in H12 from the initial inactive form towards the active form. Then, based on the determination of the active and inactive forms using the relative H12 distances, the receptor activation energies were calculated using the standard thermodynamic relation ΔG = −RTln(Nactive/Ninactive) [55], resulting 1.48 Kcal/mol for ERα1_OHT and 3.51 Kcal/mol for ERα2-E2.

2.3. Energy Analysis

A molecular mechanics-generalized Born/surface area (MM-GB/SA) approach was used to calculate the binding free energies for ERα1_E2, ERα2_E2, ERα1_OHT, and ERα2_OHT complexes. The binding free energy for the representative structure from each cluster of ERα1_E2, ERα1_OHT, ERα2_E2, and ERα2_OHT complexes are depicted in Figure 6. The binding free energy analysis plot shows different patterns of energy travel for complexes.
For ERα1_E2 and ERα2_E2, the free energy values shifted from −86.70 to −97.43 Kcal/mol and from −92.87 to −96.03 Kcal/mol in the MD simulations, respectively (top of Figure 6). In the early simulation time, a higher free energy value −86.70 Kcal/mol was observed for the complex ERα1_E2 than for the complex ERα2_E2 which had a free energy value of −92.87 Kcal/mol. At the end of the simulations, the complex ERα1_E2 had a lower free energy value (−97.43 Kcal/mol) than the complex ERα2_E2 (−96.03 Kcal/mol). The dynamic patterns of free energy revealed that E2 can accommodate better in the LBP of the ER agonist conformer than the ER antagonist conformer.
For ERα1_OHT and ERα2_OHT, the free energy travels from −134.20 to −123.20 Kcal/mol and from −106.31 to −101.07 Kcal/mol, respectively (bottom of Figure 6). The increase in free energy of ERα2_OHT in the simulation was smaller compared to ERα1_OHT. The larger increase in free energy for ERα1_OHT in the simulation might be due to H12 changing from an agonist conformation to an antagonist conformation. The analysis of free energy pattern in the simulations indicated that OHT had stronger binding than E2 in the hydrophobic LBP of ERα.
The dynamic binding interaction pattern analysis revealed that E2 and OHT bind tightly in the LBP of ERα1 and ERα2. Throughout the simulation, 12 residues of ERα1 interacted with E2. Among the 12 residues, His524 and Glu353 had hydrogen bond interactions in more than 60% and less than 20% of the simulation time, respectively. Phe404 had hydrophobic interactions in more than 60% of the simulation time. Though Ala350, Leu384, Met388, Leu391, Ile424 and Leu525 had hydrophobic interactions, the interactions were observed in less than 20% of the simulation time (Figure 7). The hydrogen bond between E2 and His524 and the hydrophobic interaction between E2 and Phe404 were the most stable in ERα1_E2 complex, while the other observed interactions were transient. Thirty-three residues of ERα2 formed interactions with OHT. Among the 33 residues, Glu353 had hydrogen bond interactions with OHT in more than 70% of the simulation time, while Thr347, Arg394 and Phe404 had hydrogen bond interactions with OHT in less than 20% of the simulation time. Ala350 and Leu525 had hydrophobic interactions with OHT in more than 40% of the simulation time, while Met343, Leu346, Leu354, Trp383, Leu384, Leu387, Met388, Leu391, Phe404, Met421, Ile424, Leu428 and His525 had hydrophobic interactions with OHT in less than 30% of the simulation time (Figure 7). Therefore, the hydrogen bond between OHT and Glu353 and hydrophobic interactions between OHT and Ala350 and Leu525 were more stable than other observed interactions. The dynamic binding interactions pattern analysis revealed that OHT tightly binds ERα2.

3. Discussion

ERα is one of the well-studied targets in the endocrine system. It plays a major role in various diseases such as cancer, bone modeling, reproductive system, cardiovascular system, metabolism, and cell proliferation [3]. Until now, many estrogenic activity chemicals were identified and crystallized in the binding site of the ER. Various computational and experimental studies were carried out to predict estrogenic activity of chemicals such as agonist, antagonist, binder, or non-binder. However, the ERα dynamic binding patterns for agonist and antagonist remain unclear. Hence, in this study, various computational techniques were applied to gain insight to the dynamic binding patterns of agonist and antagonist in the binding site of ERα. The 3D structure of ER complexed with E2 and OHT were retrieved from the PDB. The ERα1 and ERα2 represents the ER conformation in presence of E2 and OHT, respectively. The ERα1 and ERα2 were used as the target proteins to dock the E2 and OHT using QPLD.
E2 and OHT were redocked in the binding sites of ERα1 and ERα2, respectively, to validate the docking procedure and parameters [56]. The complexes from redocking were superimposed with the X-ray crystal structures 1GWR and 3ERT, respectively, to calculate RMSD values for the ligands (Figure 8). The RMSD values for E2 and OHT were 0.309 Å and 0.439 Å, respectively. The small RMSD values indicate that the docking procedure used and parameters are reliable.
The XP scores of E2 in the LBP of ERα1 and ERα2 were −11.00 and −9.65 Kcal/mol, respectively. The ERα1 and ERα2 complex with OHT had XP scores −9.07 and −8.59 Kcal/mol, respectively. Based on the Glide XP scores, E2 and OHT can bind tightly in the binding site of ERα1 compared to ERα2. The analysis based on the QPLD scores showed a different result compared to the Glide XP scores. Based on the QPLD scores, E2 and OHT can tightly bind in the binding site of ERα1 and ERα2, respectively. Applying ab initio charges can significantly increase the predictive power of the orientation of the compounds in the binding site of a protein. Hence, the ERα complex with E2 and OHT were selected based on the QPLD scores for interaction analysis instead of the Glide XP docking scores. The selected four ER complexes (ERα1_E2, ERα1_OHT, ERα2_E2, and ERα2_OHT) were subjected to MD simulations.
The in-house MATLAB script was used to calculate the RMSD matrixes from MD simulation trajectory files of the ERα complexes. The heat maps were generated based on the RMSD matrixes to identify structure clusters in the trajectory files (Figure 4). The heat maps showed various patterns of structural changes in the MD simulations of the four complexes. Hence, the hierarchical clustering method was used to group similar conformations of the ER complex from each trajectory based on its RMSD values. The clustering analysis revealed 6, 4, 5, and 4 clusters for ERα1_E2, ERα1_OHT, ERα2_OHT, and ERα2_E2, respectively. The binding free energy value was calculated for the representative structures from each cluster. For ERα1_E2 and ERα2_E2, the free energy values shifted from −86.70 to −97.43 Kcal/mol and −92.87 to −96.03 Kcal/mol, respectively. For ERα1_OHT and ERα2_OHT, the free energy travels from −134.20 to −123.20 Kcal/mol and from −106.31 to −101.07 Kcal/mol, respectively. The free energy analysis revealed that OHT had tighter binding than E2 in the hydrophobic ligand pocket of ERα1 and ERα2. So, to dissociate OHT from ERα higher external energy is required than for E2. Thus, the antagonist stays longer in the LBP of ER and represses its function.

4. Materials and Methods

4.1. Study Design

The overall workflow of this study is depicted in Figure 9. The ER complexed with E2 (PDB ID:1GWR, agonist) and OHT (PDB ID:3ERT, antagonist) were downloaded from the PDB. ER in the agonist and antagonist conformations was named ERα1 and ERα2, respectively. The QPLD method was applied to re- and cross-dock E2 and OHT in the LBP of ERα1 and ERα2 conformation. Four ER complexes (ERα1_E2, ERα1_OHT, ERα2_E2, and ERα2_OHT) were obtained from the molecular docking. The four ER complexes were subjected to 100-nanosecond (ns) MD simulations using DESMOND (https://www.schrodinger.com/products/desmond, (accessed on 8 August 2021)-Maestro-Desmond v-44017 Interoperability Tools, Schrödinger, New York, NY, USA). A MATLAB script was written to generate a RMSD matrix from each trajectory file. The hierarchical clustering analysis was carried out based on RMSD values and a representative structure was selected from each cluster. The binding free energy value was calculated by prime MM-GB/SA for each representative structure of the ER complexes from clustering analysis.

4.2. Molecular Docking

4.2.1. Protein Preparation WIZARD

The E2 complex structures, 1GWR and 3ERT, were downloaded from the PDB. 1GWR is a dimer form of ER complexed with E2. In this study, only the chain A with E2 was selected. Protein Preparation Wizard from Maestro in Schrodinger suite (https://www.schrodinger.com/products/protein-preparation-wizard, (accessed on 8 August 2021) version 2015-4 Schrödinger, LLC, New York, NY, USA) was used to re-move the heteroatoms and water molecules which are not interacting with binding site residues from 1GWR_A (ERα1) and 3ERT (ERα2). The Prime module in Schrodinger suite (https://www.schrodinger.com/products/prime, (accessed on 8 August 2021) version 2015-4 Schrödinger, LLC, New York, NY, USA) was used to add proper hydrogen atoms and to build the missing atom or residues in ER. Subsequently, the structures were optimized and minimized by applying the OPLS-AA-2005 force field [57]. The minimized ERα1 and ERα2 structures were used as the receptors to dock E2 and OHT. A 2 Å grid box was generated around the E2 and OHT in the binding site of ERα1 and ERα2, respectively.

4.2.2. Ligand Preparation

E2 and OHT were extracted from the structures 1GWR and 3ERT, respectively. These two compounds were imported into the LigPrep module in Schrodinger suite (https://www.schrodinger.com/products/ligprep, (accessed on 8 August 2021) version 2015-4 Schrödinger, LLC, New York, NY, USA) to check the bond order, ionization states, steric isomers and search the tautomers. The conformers were generated using the ConfGen method in the LigPrep module with distance dependent dielectric solvation and OPLS-AA_2005 force field used for energy minimization. The prepared E2 and OHT were saved in the SDF format.

4.2.3. QPLD Molecular Docking

QPLD is one of the powerful docking methods to identify the best orientation of a chemical in the binding site of a protein. As illustrated in Figure 10, QPLD combines the Glide docking algorithm with QM/MM calculations by Q-site program which uses the Jaguar and Impact program for the ligand with binding site residues (QM) and the remaining regions of the protein (MM), respectively. Subsequently, the ab initio charge was applied for the chemicals calculated using the Q-site which uses the 6-31G**/LACVP* basis set, B3LYP and Ultrafine SCF accuracy level for the density function calculation.
Initially, the Glide molecular docking module was used to dock E2 and OHT in the LBP of ERα1 and ERα2 using XP mode. GRID files were prepared by selecting residues 5 Å around the ligands. QM-based charge generated by DFT method was incorporated. The poses generated from Glide docking for the four complex structures (ERα1_E2, ERα1_OHT, ERα2_E2, and ERα2_OHT) by XP mode were subjected to QPLD to redock E2 and OHT in the LBP of ER1α and ER2α. First, QPLD generates several unique ER complexes with E2 and OHT. The QM regions were assigned for the ligand and the residues within 5 Å to the ligands (Figure 11).
The generated ER complexes were automatically subjected to the Q-site to compute the single point energy for each ER complex. Once the charge was calculated, E2 and OHT were redocked in the LBP of ERα1 and ERα2 by Glide module. Finally, 20 different complexes with XP score, docking energy, and the QPLD score were calculated.

4.3. MD Simulations

MD simulations were used to optimize the conformation of the docked ER complexes by analyzing the movements of atoms. The best conformation of ERα1_E2, ERα1_OHT, ERα2_E2, and ERα2_OHT from QPLD were subjected to 100 ns MD simulations using DESMOND (https://deshawresearch.com, (accessed on 8 August 2021) Maestro-Desmond v-44017 Interoperability Tools, Schrödinger, New York, NY, USA). The simulation systems were prepared by applying the OPLS-AA-2005 force field to the ER and ligands (E2 and OHT) and an orthorhombic box was generated by 10 Å border from the ER complexes. The accuracy of MD simulation is impacted by the force field used for the components in a simulation system. A specific force field should be derived and used for a nonstandard molecule [58]. The OPLS parameters are optimized for a variety of structural features in small molecules such as E2 and OHT to reproduce the thermodynamic properties in the liquid state [59,60,61]. The orthorhombic box was filled with TIP3P water molecules. To neutralize the whole simulation system, 0.15 M Na+ and Cl were added based on the total charge of the ER complex. The ER simulation systems were subjected to two-step energy minimizations with and without the restraint applied on the solute. As a first step of energy minimization, the ER complex was subjected to 12 ps of NVT simulation carried out at 10 K with Berendsen barostat. A restraint was applied to all the heavy atoms of the ER complexes, followed by 12 ps of NPT simulations carried out at 1 atmospheric pressure and 10 K with Berendsen barostat. The temperature and pressure were kept constant at 310 K for 100 ns. Further, the relaxed systems were subjected to a 100-ns simulation with a time set of two femtoseconds. The final trajectory files were saved for every 4.8 ps and used for subsequent analysis.
The trajectory files were used to elucidate the ER molecular mechanism in the presence of E2 and OHT. Initially, the RMSD matrices were created using the MATLAB script. The structures in each trajectory file were clustered based on the RMSD values calculated between them. The clusters with more than 500 structures were selected. A representative structure was obtained from each cluster for energy analysis.

4.4. Energy Calculation

The representative structures of ERα1 and ERα2 complexes from the clustering analysis were used to calculate the binding free energy. The MM-GB/SA method and OPLS-AA_2005 force field were used to compute the electrostatic component of the solvation free energy. This approach combines the molecular mechanism and continuum solvent models to predict the protein–ligand binding free energy as illustrated by the thermodynamic cycle shown in in Figure 12.
First, the receptor (ERα1 or ERα2), ligand (E2 or OHT), and complex of a receptor bound by a ligand were optimized in the solvent environment. The optimized structures were then used to calculate energy terms such as coulomb energy, covalent binding energy, Van der Waals energy, lipophilic energy, generalized Born electrostatic solvation energy, hydrogen bonding correction, and pi-pi stacking correction. These energy terms for the receptor, ligand, and complex were summed up as free energy as shown in equations below.
E P r o t e i n   =   E P r o t e i n C o u l o m b + E P r o t e i n   C o v a l e n t + E P r o t e i n H b o n d + E P r o t e i n L i p o + E P r o t e i n S o l v _ G B + E P r o t e i n v d W + E P r o t e i n P a c k i n g
E L i g a n d   =   E L i g a n d C o u l o m b + E L i g a n d   C o v a l e n t + E L i g a n d H b o n d + E L i g a n d L i p o + E L i g a n d S o l v _ G B + E L i g a n d v d W + E P r o t e i n P a c k i n g
E C o m p l e x   =   E C o m p l e x C o u l o m b + E c o m p l e x   C o v a l e n t + E C o m p l e x H b o n d + E C o m p l e x L i p o + E C o m p l e x S o l v _ G B + E C o m p l e x v d W + E C o m p l e x P a c k i n g
ECoulomb, ECovalent, EHbond, ELipo, ESolv_GB, EvdW, and EPacking represent coulomb energy, covalent binding energy, hydrogen bonding correction, lipophilic energy, generalized born electrostatic solvation energy, Van der Waals energy, and pi-pi packing correction, respectively. The binding free energy is estimated from the free energies of the ligand, protein and the complex using the following equations.
ΔGbind = EcomplexEproteinEligand
ΔGbind is the estimated binding free energy. Ecomplex is the estimated free energy of the complex in solvent, EProtein is the free energy of the protein in solvent, and Eligand is the free energy of the ligand in solvent.

5. Conclusions

ER is one of the important endocrine targets in the endocrine system. Overexpression of ER leads to various diseases such as cancer. Hence, understanding the dynamic binding patterns of agonist and antagonist binding in the hydrophobic binding pocket of ER gives insight into designing and predicting effective estrogenic compounds. Here, we applied an approach combining QM/MM docking and MD simulations to determine the energy-based ER agonist and antagonist binding mechanisms. The molecular docking revealed that E2 and OHT tightly bind in the LBP of ERα1 and ERα2, respectively. QPLD more accurately predicted the binding orientation of the estrogenic compounds in the LBP of ER than Glide docking. The heat maps RMSD values revealed that different clusters formed for the structures during the MD simulations and the binding mechanisms of agonist and antagonist in ERα were different. The binding free energy analysis for the representative structures of the clusters revealed that OHT binds more tightly with ERα2, while E2 prefers to bind ERα1. The binding interaction analysis revealed that the hydrogen bond between OHT and Glu353 and the hydrophobic interaction between OHT and Ala350 and Leu525 are the most stable interactions in the ERα2_OHT complex. These interactions might be the reason for the tight binding of OHT rather than E2 in the LBP of ERα1 and ERα2. Our findings shed light on the structural basis of agonist and antagonist dynamic binding pattern. This insight into the binding pattern and free energy analysis may help to discover ERα agonists and antagonists as drug candidates and facilitate risk assessment of chemicals through ER-mediated responses.

Author Contributions

Conceptualization, H.H.; methodology, S.S., C.S. and H.H.; software, S.S., C.S., W.G. (Weigong Ge) and H.H.; validation, S.S., W.G. (Wenjing Guo) and J.L.; formal analysis, S.S., C.S. and H.H.; data curation, S.S. and C.S.; writing—original draft preparation, S.S. and H.H.; writing—review and editing, T.A.P.; supervision, H.H. All authors have read and agreed to the published version of the manuscript.

Funding

This research received no external funding.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Not applicable.

Acknowledgments

This research was supported in part by an appointment to the Research Participation Program at the National Center for Toxicological Research (Chandrabose Selvaraj) administered by the Oak Ridge Institute for Science and Education through an interagency agreement between the U.S. Department of Energy and the U.S. Food and Drug Administration.

Conflicts of Interest

The authors declare no conflict of interest.

Disclaimer

The views presented in this article do not necessarily reflect those of the US Food and Drug Administration.

References

  1. Carroll, J.S.; Brown, M. Estrogen Receptor Target Gene: An Evolving Concept. Mol. Endocrinol. 2006, 20, 1707–1714. [Google Scholar] [CrossRef]
  2. Germain, P.; Staels, B.; Dacquet, C.; Spedding, M.; Laudet, V. Overview of nomenclature of nuclear receptors. Pharmacol. Rev. 2006, 58, 685–704. [Google Scholar] [CrossRef]
  3. Schug, T.T.; Janesick, A.; Blumberg, B.; Heindel, J.J. Endocrine disrupting chemicals and disease susceptibility. J. Steroid Biochem. Mol. Biol. 2011, 127, 204–215. [Google Scholar] [CrossRef] [Green Version]
  4. Kampa, M.; Pelekanou, V.; Notas, G.; Stathopoulos, E.N.; Castanas, E. The estrogen receptor: Two or more molecules, multiple variants, diverse localizations, signaling and functions. Are we undergoing a paradigm-shift as regards their significance in breast cancer? Hormones 2013, 12, 69–85. [Google Scholar] [CrossRef] [PubMed]
  5. Ascenzi, P.; Bocedi, A.; Marino, M. Structure-function relationship of estrogen receptor alpha and beta: Impact on human health. Mol. Asp. Med. 2006, 27, 299–402. [Google Scholar] [CrossRef] [PubMed]
  6. Shiau, A.K.; Barstad, D.; Loria, P.M.; Cheng, L.; Kushner, P.J.; Agard, D.A.; Greene, G.L. The structural basis of estrogen receptor/coactivator recognition and the antagonism of this interaction by tamoxifen. Cell 1998, 95, 927–937. [Google Scholar] [CrossRef] [Green Version]
  7. Brzozowski, A.M.; Pike, A.C.W.; Dauter, Z.; Hubbard, R.E.; Bonn, T.; Engström, O.; Öhman, L.; Greene, G.L.; Gustafsson, J.-Å.; Carlquist, M. Molecular basis of agonism and antagonism in the oestrogen receptor. Nature 1997, 389, 753–758. [Google Scholar] [CrossRef] [PubMed]
  8. Rastinejad, F.; Huang, P.; Chandra, V.; Khorasanizadeh, S. Understanding nuclear receptor form and function using structural biology. J. Mol. Endocrinol. 2013, 51, T1–T21. [Google Scholar] [CrossRef] [Green Version]
  9. Huang, P.; Chandra, V.; Rastinejad, F. Structural overview of the nuclear receptor superfamily: Insights into physiology and therapeutics. Annu. Rev. Physiol. 2010, 72, 247–272. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  10. Yaşar, P.; Ayaz, G.; User, S.D.; Güpür, G.; Muyan, M. Molecular mechanism of estrogen-estrogen receptor signaling. Reprod. Med. Biol. 2016, 16, 4–20. [Google Scholar] [CrossRef] [PubMed]
  11. Ng, H.W.; Perkins, R.; Tong, W.; Hong, H. Versatility or promiscuity: The estrogen receptors, control of ligand selectivity and an update on subtype selective ligands. Int. J. Environ. Res. Public Health 2014, 11, 8709–8742. [Google Scholar] [CrossRef] [Green Version]
  12. Pike, A.C. Lessons learnt from structural studies of the oestrogen receptor. Best Pract. Res. Clin. Endocrinol. Metab. 2006, 20, 1–14. [Google Scholar] [CrossRef] [PubMed]
  13. Nettles, K.W.; Bruning, J.B.; Gil, G.; O’Neill, E.E.; Nowak, J.; Guo, Y.; Kim, Y.; DeSombre, E.R.; Dilis, R.; Hanson, R.N.; et al. Structural plasticity in the oestrogen receptor ligand-binding domain. EMBO Rep. 2007, 8, 563–568. [Google Scholar] [CrossRef] [PubMed]
  14. Kumar, R.; Zakharov, M.N.; Khan, S.H.; Miki, R.; Jang, H.; Toraldo, G.; Singh, R.; Bhasin, S.; Jasuja, R. The dynamic structure of the estrogen receptor. J. Amino Acids 2011, 2011, 812540. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  15. Bruning, J.B.; Parent, A.A.; Gil, G.; Zhao, M.; Nowak, J.; Pace, M.C.; Smith, C.L.; Afonine, P.V.; Adams, P.D.; Katzenellenbogen, J.A.; et al. Coupling of receptor conformation and ligand orientation determine graded activity. Nat. Chem. Biol. 2010, 6, 837–843. [Google Scholar] [CrossRef] [Green Version]
  16. Bjornstrom, L.; Sjoberg, M. Mechanisms of Estrogen Receptor Signaling: Convergence of Genomic and Nongenomic Actions on Target Genes. Mol. Endocrinol. 2005, 19, 833–842. [Google Scholar] [CrossRef] [Green Version]
  17. Pakdel, F. Molecular Pathways of Estrogen Receptor Action. Int. J. Mol. Sci. 2018, 19, 2591. [Google Scholar] [CrossRef] [Green Version]
  18. Ng, H.W.; Zhang, W.; Shu, M.; Luo, H.; Ge, W.; Perkins, R.; Tong, W.; Hong, H. Competitive molecular docking approach for predicting estrogen receptor subtype alpha agonists and antagonists. BMC Bioinform. 2014, 15 (Suppl. S11), S4. [Google Scholar] [CrossRef] [Green Version]
  19. Pang, X.; Fu, W.; Wang, J.; Kang, D.; Xu, L.; Zhao, Y.; Liu, A.L.; Du, G.H. Identification of Estrogen Receptor alpha Antagonists from Natural Products via In Vitro and In Silico Approaches. Oxidative Med. Cell. Longev. 2018, 2018, 6040149. [Google Scholar] [CrossRef] [Green Version]
  20. Zhang, L.; Sedykh, A.; Tripathi, A.; Zhu, H.; Afantitis, A.; Mouchlis, V.D.; Melagraki, G.; Rusyn, I.; Tropsha, A. Identification of putative estrogen receptor-mediated endocrine disrupting chemicals using QSAR- and structure-based virtual screening approaches. Toxicol. Appl. Pharmacol. 2013, 272, 67–76. [Google Scholar] [CrossRef] [Green Version]
  21. Muchtaridi, M.; Syahidah, H.N.; Subarnas, A.; Yusuf, M.; Bryant, S.D.; Langer, T. Molecular Docking and 3D-Pharmacophore Modeling to Study the Interactions of Chalcone Derivatives with Estrogen Receptor Alpha. Pharmaceuticals 2017, 10, 81. [Google Scholar] [CrossRef] [Green Version]
  22. Dutta, S.; Kharkar, P.S.; Sahu, N.U.; Khanna, A. Molecular docking prediction and in vitro studies elucidate anti-cancer activity of phytoestrogens. Life Sci. 2017, 185, 73–84. [Google Scholar] [CrossRef]
  23. Ksiazek, P.; Bryl, K. Molecular Docking Reveals Binding Features of Estrogen Receptor Beta Selective Ligands. Curr. Comput.-Aided Drug Des. 2015, 11, 137–151. [Google Scholar] [CrossRef]
  24. Powers, C.N.; Setzer, W.N. A molecular docking study of phytochemical estrogen mimics from dietary herbal supplements. Silico Pharmacol. 2015, 3, 4. [Google Scholar] [CrossRef] [Green Version]
  25. Fratev, F. Activation helix orientation of the estrogen receptor is mediated by receptor dimerization: Evidence from molecular dynamics simulations. Phys. Chem. Chem. Phys. 2015, 17, 13403–13420. [Google Scholar] [CrossRef] [PubMed]
  26. Hu, G.; Wang, J. Ligand selectivity of estrogen receptors by a molecular dynamics study. Eur. J. Med. Chem. 2014, 74, 726–735. [Google Scholar] [CrossRef]
  27. Lu, Q.; Cai, Z.; Fu, J.; Luo, S.; Liu, C.; Li, X.; Zhao, D. Molecular docking and molecular dynamics studies on the interactions of hydroxylated polybrominated diphenyl ethers to estrogen receptor alpha. Ecotoxicol. Environ. Saf. 2014, 101, 83–89. [Google Scholar] [CrossRef] [PubMed]
  28. Ng, H.L. Simulations reveal increased fluctuations in estrogen receptor-alpha conformation upon antagonist binding. J. Mol. Graph. Model. 2016, 69, 72–77. [Google Scholar] [CrossRef] [PubMed]
  29. Sakkiah, S.; Kusko, R.; Tong, W.; Hong, H. Applications of Molecular Dynamics Simulations in Computational Toxicology. In Advances in Computational Toxicology: Methodologies and Applications in Regulatory Science; Hong, H., Ed.; Springer International Publishing: Cham, Switzerland, 2019; pp. 181–212. [Google Scholar] [CrossRef]
  30. Selvaraj, C.; Sakkiah, S.; Tong, W.; Hong, H. Molecular dynamics simulations and applications in computational toxicology and nanotoxicology. Food Chem. Toxicol. 2018, 112, 495–506. [Google Scholar] [CrossRef] [PubMed]
  31. Mansouri, K.; Abdelaziz, A.; Rybacka, A.; Roncaglioni, A.; Tropsha, A.; Varnek, A.; Zakharov, A.; Worth, A.; Richard, A.M.; Grulke, C.M.; et al. CERAPP: Collaborative Estrogen Receptor Activity Prediction Project. Environ. Health Perspect. 2016, 124, 1023–1033. [Google Scholar] [CrossRef]
  32. Sakkiah, S.; Guo, W.; Pan, B.; Kusko, R.; Tong, W.; Hong, H. Computational prediction models for assessing endocrine disrupting potential of chemicals. J. Environ. Sci. Health C Environ. Carcinog Ecotoxicol. Rev. 2018, 36, 192–218. [Google Scholar] [CrossRef]
  33. Norinder, U.; Boyer, S. Conformal Prediction Classification of a Large Data Set of Environmental Chemicals from ToxCast and Tox21 Estrogen Receptor Assays. Chem. Res. Toxicol. 2016, 29, 1003–1010. [Google Scholar] [CrossRef] [PubMed]
  34. Ng, H.W.; Doughty, S.W.; Luo, H.; Ye, H.; Ge, W.; Tong, W.; Hong, H. Development and Validation of Decision Forest Model for Estrogen Receptor Binding Prediction of Chemicals Using Large Data Sets. Chem. Res. Toxicol. 2015, 28, 2343–2351. [Google Scholar] [CrossRef] [PubMed]
  35. He, J.; Peng, T.; Yang, X.; Liu, H. Development of QSAR models for predicting the binding affinity of endocrine disrupting chemicals to eight fish estrogen receptor. Ecotoxicol. Environ. Saf. 2018, 148, 211–219. [Google Scholar] [CrossRef]
  36. Ng, H.W.; Shu, M.; Luo, H.; Ye, H.; Ge, W.; Perkins, R.; Tong, W.; Hong, H. Estrogenic activity data extraction and in silico prediction show the endocrine disruption potential of bisphenol A replacement compounds. Chem. Res. Toxicol. 2015, 28, 1784–1795. [Google Scholar] [CrossRef] [PubMed]
  37. Bhhatarai, B.; Wilson, D.M.; Price, P.S.; Marty, S.; Parks, A.K.; Carney, E. Evaluation of OASIS QSAR Models Using ToxCast in Vitro Estrogen and Androgen Receptor Binding Data and Application in an Integrated Endocrine Screening Approach. Environ. Health Perspect. 2016, 124, 1453–1461. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  38. Niu, A.Q.; Xie, L.J.; Wang, H.; Zhu, B.; Wang, S.Q. Prediction of selective estrogen receptor beta agonist using open data and machine learning approach. Drug Des. Dev. Ther. 2016, 10, 2323–2331. [Google Scholar] [CrossRef] [Green Version]
  39. Ribay, K.; Kim, M.T.; Wang, W.; Pinolini, D.; Zhu, H. Predictive Modeling of Estrogen Receptor Binding Agents Using Advanced Cheminformatics Tools and Massive Public Data. Front. Environ. Sci. 2016, 4, 12. [Google Scholar] [CrossRef] [Green Version]
  40. Zhao, Q.; Lu, Y.; Zhao, Y.; Li, R.; Luan, F.; Cordeiro, M.N. Rational Design of Multi-Target Estrogen Receptors ERalpha and ERbeta by QSAR Approaches. Curr. Drug Targets 2017, 18, 576–591. [Google Scholar] [CrossRef] [PubMed]
  41. Browne, P.; Judson, R.S.; Casey, W.M.; Kleinstreuer, N.C.; Thomas, R.S. Screening Chemicals for Estrogen Receptor Bioactivity Using a Computational Model. Environ. Sci. Technol. 2015, 49, 8804–8814. [Google Scholar] [CrossRef]
  42. Sosnovcova, J.; Rucki, M.; Bendova, H. Estrogen Receptor Binding Affinity of Food Contact Material Components Estimated by QSAR. Cent. Eur. J. Public Health 2016, 24, 241–244. [Google Scholar] [CrossRef] [Green Version]
  43. Blair, R.M.; Fang, H.; Branham, W.S.; Hass, B.S.; Dial, S.L.; Moland, C.L.; Tong, W.; Shi, L.; Perkins, R.; Sheehan, D.M. The Estrogen Receptor Relative Binding Affinities of 188 Natural and Xenochemicals: Structural Diversity of Ligands. Toxicol. Sci. 2000, 54, 138–153. [Google Scholar] [CrossRef] [Green Version]
  44. Yang, R.; Li, N.; Rao, K.; Ma, M.; Wang, Z. Combined action of estrogen receptor agonists and antagonists in two-hybrid recombinant yeast in vitro. Ecotoxicol. Environ. Saf. 2015, 111, 228–235. [Google Scholar] [CrossRef]
  45. Puranik, N.V.; Srivastava, P.; Bhatt, G.; John Mary, D.J.S.; Limaye, A.M.; Sivaraman, J. Determination and analysis of agonist and antagonist potential of naturally occurring flavonoids for estrogen receptor (ERα) by various parameters and molecular modelling approach. Sci. Rep. 2019, 9, 7450. [Google Scholar] [CrossRef] [PubMed]
  46. Sakkiah, S.; Selvaraj, C.; Gong, P.; Zhang, C.; Tong, W.; Hong, H. Development of estrogen receptor beta binding prediction model using large sets of chemicals. Oncotarget 2017, 8, 92989–93000. [Google Scholar] [CrossRef] [PubMed]
  47. Shen, J.; Xu, L.; Fang, H.; Richard, A.M.; Bray, J.D.; Judson, R.S.; Zhou, G.; Colatsky, T.J.; Aungst, J.L.; Teng, C.; et al. EADB: An estrogenic activity database for assessing potential endocrine activity. Toxicol. Sci. 2013, 135, 277–291. [Google Scholar] [CrossRef] [Green Version]
  48. Shi, L.; Tong, W.; Fang, H.; Xie, Q.; Hong, H.; Perkins, R.; Wu, J.; Tu, M.; Blair, R.M.; Branham, W.S.; et al. An integrated “4-phase” approach for setting endocrine disruption screening priorities--phase I and II predictions of estrogen receptor binding affinity. SAR QSAR Environ. Res. 2002, 13, 69–88. [Google Scholar] [CrossRef] [PubMed]
  49. Hong, H.; Tong, W.; Fang, H.; Shi, L.; Xie, Q.; Wu, J.; Perkins, R.; Walker, J.D.; Branham, W.; Sheehan, D.M. Prediction of estrogen receptor binding for 58,000 chemicals using an integrated system of a tree-based model with structural alerts. Environ. Health Perspect. 2002, 110, 29–36. [Google Scholar] [CrossRef] [Green Version]
  50. Jordan, V.C.; Collins, M.M.; Rowsby, L.; Prestwich, G. A monohydroxylated metabolite of tamoxifen with potent antioestrogenic activity. J. Endocrinol. 1977, 75, 305–316. [Google Scholar] [CrossRef]
  51. Wang, Y.; Chirgadze, N.Y.; Briggs, S.L.; Khan, S.; Jensen, E.V.; Burris, T.P. A second binding site for hydroxytamoxifen within the coactivator-binding groove of estrogen receptor β. Proc. Natl. Acad. Sci. USA 2006, 103, 9908–9911. [Google Scholar] [CrossRef] [Green Version]
  52. Cho, A.E.; Guallar, V.; Berne, B.J.; Friesner, R. Importance of accurate charges in molecular docking: Quantum mechanical/molecular mechanical (QM/MM) approach. J. Comput. Chem. 2005, 26, 915–931. [Google Scholar] [CrossRef] [Green Version]
  53. Zhong, H.; Kirschner, K.N.; Lee, M.; Bowen, J.P. Binding free energy calculation for duocarmycin/DNA complex based on the QPLD-derived partial charge model. Bioorg. Med. Chem. Lett. 2008, 18, 542–545. [Google Scholar] [CrossRef] [PubMed]
  54. Reema Abu, K.; Hamada Abd, E.-A.; Dima, S.; Ghadeer, A.; Ghassan Abu, S. CETP Inhibitory Activity of Chlorobenzyl Benzamides: QPLD Docking, Pharmacophore Mapping and Synthesis. Lett. Drug Des. Discov. 2017, 14, 1391–1400. [Google Scholar] [CrossRef]
  55. Udommaneethanakit, T.; Rungrotmongkol, T.; Frecer, V.; Seneci, P.; Miertus, S.; Bren, U. Drugs against avian influenza a virus: Design of novel sulfonate inhibitors of neuraminidase N1. Curr. Pharm. Des. 2014, 20, 3478–3487. [Google Scholar] [CrossRef] [PubMed]
  56. Furlan, V.; Bren, U. Insight into Inhibitory Mechanism of PDE4D by Dietary Polyphenols Using Molecular Dynamics Simulations and Free Energy Calculations. Biomolecules 2021, 11, 479. [Google Scholar] [CrossRef] [PubMed]
  57. Robertson, M.J.; Qian, Y.; Robinson, M.C.; Tirado-Rives, J.; Jorgensen, W.L. Development and Testing of the OPLS-AA/M Force Field for RNA. J. Chem. Theory Comput. 2019, 15, 2734–2742. [Google Scholar] [CrossRef] [PubMed]
  58. Lešnik, S.; Hodošček, M.; Bren, U.; Stein, C.; Bondar, A.-N. Potential Energy Function for Fentanyl-Based Opioid Pain Killers. J. Chem. Inf. Modeling 2020, 60, 3566–3576. [Google Scholar] [CrossRef] [PubMed]
  59. Jorgensen, W.L.; Maxwell, D.S.; Tirado-Rives, J. Development and Testing of the OPLS All-Atom Force Field on Conformational Energetics and Properties of Organic Liquids. J. Am. Chem. Soc. 1996, 118, 11225–11236. [Google Scholar] [CrossRef]
  60. Storer, J.W.; Giesen, D.J.; Cramer, C.J.; Truhlar, D.G. Class IV charge models: A new semiempirical approach in quantum chemistry. J. Comput. Aided Mol. Des. 1995, 9, 87–110. [Google Scholar] [CrossRef]
  61. Shivakumar, D.; Harder, E.; Damm, W.; Friesner, R.A.; Sherman, W. Improving the Prediction of Absolute Solvation Free Energies Using the Next Generation OPLS Force Field. J. Chem. Theory Comput. 2012, 8, 2553–2558. [Google Scholar] [CrossRef]
Figure 1. The cartoon representation of the different domains of ER. LBP is shown in green; AF2 is shown in electrostatic representation; and H12 is shown in magenta and blue for active (agonist bound) and inactive (antagonist bound) forms, respectively.
Figure 1. The cartoon representation of the different domains of ER. LBP is shown in green; AF2 is shown in electrostatic representation; and H12 is shown in magenta and blue for active (agonist bound) and inactive (antagonist bound) forms, respectively.
Ijms 22 09371 g001
Figure 2. Two-dimensional structures of 17β-estradiol (E2) and 4-hydroxytamoxifen (OHT).
Figure 2. Two-dimensional structures of 17β-estradiol (E2) and 4-hydroxytamoxifen (OHT).
Ijms 22 09371 g002
Figure 3. Binding orientations of E2 and OHT in the hydrophobic binding pocket of ERα1 and ERα2. The dotted lines represent hydrogen bond interactions. Green circles represent the hydrophobic residues; cyan circles represent the polar residues; red and blue ovals represent the negative and positive charged residues, respectively.
Figure 3. Binding orientations of E2 and OHT in the hydrophobic binding pocket of ERα1 and ERα2. The dotted lines represent hydrogen bond interactions. Green circles represent the hydrophobic residues; cyan circles represent the polar residues; red and blue ovals represent the negative and positive charged residues, respectively.
Ijms 22 09371 g003
Figure 4. RMSD matrixes from MD simulations of the four complexes: (A) ERα1_E2, (B) ERα1_OHT, (C) ERα2_OHT, and (D) ERα2_E2. The complex is shown in the title of each panel. The time of a structure in the MD simulations is depicted by the axis. The RMSD values are color coded as shown in the subfigure color legends.
Figure 4. RMSD matrixes from MD simulations of the four complexes: (A) ERα1_E2, (B) ERα1_OHT, (C) ERα2_OHT, and (D) ERα2_E2. The complex is shown in the title of each panel. The time of a structure in the MD simulations is depicted by the axis. The RMSD values are color coded as shown in the subfigure color legends.
Ijms 22 09371 g004
Figure 5. Relative distances of H12 to the centroid of LBP of ER in the simulations of ERα1_OHT (A) and ERα2_E2 (B). The X-axis indicates simulation time in ns. The Y-axis shows the relative H12 distance compared to the initial structures.
Figure 5. Relative distances of H12 to the centroid of LBP of ER in the simulations of ERα1_OHT (A) and ERα2_E2 (B). The X-axis indicates simulation time in ns. The Y-axis shows the relative H12 distance compared to the initial structures.
Ijms 22 09371 g005
Figure 6. Binding free energy analysis on the average structure from each cluster of ER complexes ERα2_E2 (brown in A), ERα1_E2 (blue in A), ERα1_OHT (brown in B), and ERα2_OHT (blue in B) calculated using Prime MM-GB/SA.
Figure 6. Binding free energy analysis on the average structure from each cluster of ER complexes ERα2_E2 (brown in A), ERα1_E2 (blue in A), ERα1_OHT (brown in B), and ERα2_OHT (blue in B) calculated using Prime MM-GB/SA.
Ijms 22 09371 g006
Figure 7. Interactions between E2 and ERα1 (top) and between OHT and ERα2 (bottom). The X-axis represents amino acids and their numbering in ER. The Y-axis shows the fractions of simulation time for different interaction types which marked in different colors.
Figure 7. Interactions between E2 and ERα1 (top) and between OHT and ERα2 (bottom). The X-axis represents amino acids and their numbering in ER. The Y-axis shows the fractions of simulation time for different interaction types which marked in different colors.
Ijms 22 09371 g007
Figure 8. The superimposition of the redocked and X-ray crystal structure of ER complexes. (A) Docked ERα1_E2 complex and X-ray crystal structure (1GWR). (B) Docked ERα2_OHT complex and X-ray crystal structure (3ERT). The ER are shown in the stylized model. E2 and OHT are displayed in the stick model. X-ray structures are in cyan, docked ERα1_E2 in green and ERα2_OHT in magenta.
Figure 8. The superimposition of the redocked and X-ray crystal structure of ER complexes. (A) Docked ERα1_E2 complex and X-ray crystal structure (1GWR). (B) Docked ERα2_OHT complex and X-ray crystal structure (3ERT). The ER are shown in the stylized model. E2 and OHT are displayed in the stick model. X-ray structures are in cyan, docked ERα1_E2 in green and ERα2_OHT in magenta.
Ijms 22 09371 g008
Figure 9. Overview of the study design.
Figure 9. Overview of the study design.
Ijms 22 09371 g009
Figure 10. QPLD docking procedure.
Figure 10. QPLD docking procedure.
Ijms 22 09371 g010
Figure 11. Cartoon representation of the region assigned for QM and MM during the QPLD docking. The square represents ER, binding site is highlighted in red, the QM region is depicted in yellow, and the MM region is shown in orange.
Figure 11. Cartoon representation of the region assigned for QM and MM during the QPLD docking. The square represents ER, binding site is highlighted in red, the QM region is depicted in yellow, and the MM region is shown in orange.
Ijms 22 09371 g011
Figure 12. Thermodynamic cycle for binding free energy calculation.
Figure 12. Thermodynamic cycle for binding free energy calculation.
Ijms 22 09371 g012
Table 1. Docking scores and docking energy values for the four ER complexes.
Table 1. Docking scores and docking energy values for the four ER complexes.
ER ComplexGlideQPLD
XP Score Kcal/molDocking EnergyQPLD Score Kcal/molDocking Energy
ERα1_E2−11.00−39.74−11.45−38.78
ERα2_E2−9.65−32.54−9.82−33.41
ERα1_OHT−9.07−29.43−8.17−28.96
ERα2_OHT−8.59−35.81−10.94−38.23
Table 2. Summary of the MD simulation systems.
Table 2. Summary of the MD simulation systems.
ComplexAtoms in ComplexWatersIons
ERα1_E23980789030 Na+; 23 Cl
ERα1_OHT3994789029 Na+; 22 Cl
ERα2_E23946913836 Na+; 26 Cl
ERα2_OHT3960913835 Na+; 25 Cl
Table 3. Interaction energy values for the ER complexes.
Table 3. Interaction energy values for the ER complexes.
AcronymsΔG_BindΔG_Bind_CoulombΔG_Bind_vdWLigand EnergyComplex EnergyReceptor Energy
ERα1_E2−44.10−7.73−22.171.73−9402.71−9360.34
ERα1_OHT−44.68−22.03−17.4732.43−9352.97−9340.71
ERα2_E2−24.69−9.19−7.131.72−9325.79−9302.81
ERα2_OHT−42.24−39.65−17.0432.12−9342.13−9332.01
Table 4. Major clusters from hierarchical clustering analysis.
Table 4. Major clusters from hierarchical clustering analysis.
ER ComplexCluster NumberNumber of StructuresStart FrameEnd FrameRepresentative Structure
ERα1_E214314243282312
21012432953464545
31366536467876339
47027742314,51413,922
5193814,52316,49115,558
687416,73017,62917,171
ERα1_OHT1169621698979
2860245634222866
31328589472987154
413,528730620,83515,478
ERα2_E215247254512353
26915545312,9458216
3129415,18716,74016,064
4378616,73520,83518,646
ERα2_OHT13732537451813
21531455260865356
33497608796437771
45916964115,76913,467
5298815,77019,02118,075
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Sakkiah, S.; Selvaraj, C.; Guo, W.; Liu, J.; Ge, W.; Patterson, T.A.; Hong, H. Elucidation of Agonist and Antagonist Dynamic Binding Patterns in ER-α by Integration of Molecular Docking, Molecular Dynamics Simulations and Quantum Mechanical Calculations. Int. J. Mol. Sci. 2021, 22, 9371. https://doi.org/10.3390/ijms22179371

AMA Style

Sakkiah S, Selvaraj C, Guo W, Liu J, Ge W, Patterson TA, Hong H. Elucidation of Agonist and Antagonist Dynamic Binding Patterns in ER-α by Integration of Molecular Docking, Molecular Dynamics Simulations and Quantum Mechanical Calculations. International Journal of Molecular Sciences. 2021; 22(17):9371. https://doi.org/10.3390/ijms22179371

Chicago/Turabian Style

Sakkiah, Sugunadevi, Chandrabose Selvaraj, Wenjing Guo, Jie Liu, Weigong Ge, Tucker A. Patterson, and Huixiao Hong. 2021. "Elucidation of Agonist and Antagonist Dynamic Binding Patterns in ER-α by Integration of Molecular Docking, Molecular Dynamics Simulations and Quantum Mechanical Calculations" International Journal of Molecular Sciences 22, no. 17: 9371. https://doi.org/10.3390/ijms22179371

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop