Abstract
This study aimed to identify the docking and molecular mechanics-generalized born surface area (MM-GBSA) re-scoring parameters which can correlate the binding affinity and selectivity of the ligands towards oestrogen receptor β (ERβ). Three different series of ERβ ligands were used as dataset and the compounds were docked against ERβ (protein data bank (PDB) ID: 1QKM) using Glide and ArgusLab. Glide docking showed superior results when compared with ArgusLab. Docked poses were then rescored using Prime-MM-GBSA to calculate free energy binding. Correlations were made between observed activities of ERβ ligands with computationally predicted values from docking, binding energy parameters. ERβ ligands experimental binding affinity/selectivity did not correlate well with Glide and ArgusLab score. Whereas calculated Glide energy (coulomb-van der Waal interaction energy) correlated significantly with binding affinity of ERβ ligands (r2 = 0.66). MM-GBSA re-scoring showed correlation of r2 = 0.74 with selectivity of ERβ ligands. These results will aid the discovery of novel ERβ ligands with isoform selectivity.
Introduction
Oestrogen receptors α (ERα) and β (ERβ) are involved in regulating cell growth, proliferation, and differentiation in normal and cancerous tissues. ERα and ERβ are encoded by different genes and bind with endogenous oestrogen, 17β-estradiol (E2), with similar affinity, though they differ in size and share 60% ligand binding domain identity. In addition, ERα and ERβ receptors exhibit a similar, but not identical, binding affinity profile to different ligands. In addition, studies of knockout mice have shown two subtypes of ERs with distinct functions and are expressed differentially in certain tissues. These differences have stimulated the search for subtype-specific ligands that can elicit tissue or cell-specific ER activity. In particular, ERα receptor is associated with the proliferative effects of prostate, whereas the current evidence implies that ERβ has growth-suppressive activities in prostateCitation1. Taken together, these data establish ERβ as a potential therapeutic target for cancer. However, the high similarities in the ligand binding domain of the two ERs constitute a serious impediment to design and develop compounds that show high ERβ selectivity. The availability of three-dimensional structures of ERβ co-crystallized with various agonists provides a good basis for using structure-based techniques for the discovery of new ERβ agonists.
In silico screening has become a routine component of drug discovery. To accelerate the identification of a potential candidate, very large number of ligands are now passed through various types of rapid theoretical screens. Scoring functions and/or property range filters are used to identify potential ligands for further experimental analysis. The quality of the selected ligands depends upon the accuracy of their predicted binding modes and affinities. Due to the current advancement in docking methodology, the prediction of the binding modes of ligands with reasonable accuracy is achievable. However, accurate prediction of binding affinities remains challenging. Currently, molecular dynamics (MD) and Monte Carlo (MC) simulations coupled with free energy perturbation (FEP), thermodynamic integration (TI), or umbrella sampling/weighted histogram analysis, and metadynamics are computational techniques that are routinely used to estimate ligand binding affinities. However, the longstanding issue in MD or MC simulations is the difficulty of observing large number of protein-ligand conformations, in which case entropic contributions to the binding free energy are incorrectly calculated. The solution to this predicament is simple; the precision can usually be increased by running longer simulations. Although it still requires substantial computational resources, the use of special MD engines running on graphical processing units has greatly reduced these computational costs. Regardless, the MD approach is not relevant for high throughput screening for several reasons. First, the amount of calculations required for a single compound is prohibitively scalable to thousands of ligands and second, the degree of detail provided by the MD methodology is probably not needed during the screening phase. Instead, it is recommended during the lead optimization phase encompassing tens of potential lead compounds to completely understand their binding pathways before moving forward to the next stage of ligand developmentCitation2.
Molecular simulations employing FEP and TI methods, developed within the framework of statistical mechanics, constitute some of the best options from the standpoint of theoretical rigor required for evaluating receptor-ligand binding free energies. Moreover, the statistical noise associated with series of binding free energy simulations in FEP and TI methods, diminishes their discriminatory power in making a judicious choice among alternative candidate molecules. A recent study has reported that TI could be a useful tool in drug design only if all ligands involve a common scaffoldCitation3. Even though computationally intensive FEP and TI calculations are handled routinely by the computational power available today, they are not state–of–the–art approaches anymoreCitation4. Metadynamics is a free energy method that is based on biasing of the potential surface and is similar to umbrella sampling in this sense. However, in contrast to umbrella sampling, the bias or “hills” are dynamically placed on top of the underlying potential energy landscape and discourages the system from visiting the same points in the configurational space. The method was found to be robust and as accurate as umbrella sampling with weighted histogram techniques. However, it shares with the latter some drawbacks that make it scarcely practical for routine docking tasks. In particular, it requires a careful choice of collective variables (CVs). In order to reconstruct the free energy surface, the CVs must describe all the slow events that are relevant to the ligand binding/unbinding. When a relevant CV is neglected, the reconstruction of the free energy surface fails. This often implies the need to run multiple metadynamics simulations to understand which CVs are needed for each protein–ligand complex. Moreover, even when all the needed CVs are included, the reconstruction of the free energy surface in more than three dimensions requires huge computational resourcesCitation5.
Instead, several simplified methods have been developed that are still based on simulations but only sample the end-points of the reaction i.e. conformations of the free and bound species were generated and the binding free energy was computed by taking a difference e.g. linear interaction energy (LIE), molecular mechanics – Poisson–Boltzmann surface area (MM-PBSA) and molecular mechanics – generalized born surface area (MM-GBSA). Although early applications of MM-PBSA appeared promising, it is apparently difficult to converge the energy averages reliably for this method. One reason for this difficulty may be that they encompass energy fluctuations not only of the ligand and the binding site, but also of parts of the protein remote from the binding site, which are less relevant to bindingCitation4. Recent improvements in the description of intermolecular interactions using second-generation molecular mechanics force fields and the development of methodology like MM-GBSA has gained more interest because it is composed of physically well-defined terms and contains no adjustable parameters. The successful results of this new method have made it a more attractive alternative to the FEP, TI and MM-PBSA methodologies. In an earlier study, MM-GBSA approach was evaluated to rank the relative potencies of anilinoquinazoline, anilinopyrimidine chemical classes against kinases, with the potencies range from micromolar to nanomolar potency (IC50) valuesCitation6. MM-GBSA approach has also been applied to a range of pharmaceutically relevant targets, which include CDK-2, factor Xa, thrombin, and HIV-RT7. For the ERα receptor, scoring methods like MD simulations and LIE approaches were employed in the prediction of ligand binding affinity. At the same time, for the ERβ receptor, there are only a few physics-based scoring reports for the prediction of ligand binding affinity or selectivityCitation8,Citation9. Hence, MM-GBSA “state–of–the–art” free energy calculations could provide a convincing answer to medicinal chemists by rationalizing experimental observations, and in some instances, play a predictive role in ranking ligands for a specific target. Hence, in the present study, we intend to study the performance of MM-GBSA re-scoring and other docking parameters generated from two molecular docking programs in order to predict the binding affinity as well as selectivity of the ligands towards ERβ receptor.
Methodology
Data set
ERβ selective ligand series 7-substituted 2-phenyl-benzofuransCitation10, 4-hydroxy-biphenylsCitation11 and 4-hydroxy-biphenyl-carbaldehyde oxime derivativesCitation12 reported in the literature (a total of 116 molecules) were taken for the current study (). To have uniformity and minimize experimental error, the dataset of ERβ affinity and selectivity were taken from the same research groupCitation10–12 who adopted same methodology and experimental conditions. Competitive radio ligand binding assay methodology was adopted to measure the relative ERβ binding affinity (IC50) and the ERβ selectivity (ERα IC50/ ERβ IC50) of the compounds and the affinity/selectivity values were converted into logarithmic scale (pIC50) on a molar basis. The binding affinity values (as measured by IC50 values) vary from 0.35 to 78,000 nM whereas the ERβ selectivity values vary from 0 to 108 fold ().
Figure 1. Skeletal structure of estrogenic compounds. The moieties derivatized in each compound are indicated in supplementary table S1.
Table 1. Structural features, observed and predicted values for ERβ IC50 (nM) and ERα/β of 7-substituted 2-phenyl-benzofuran derivatives (1–21).
Table 2. Structural features, observed and predicted values for ERβ IC50 (nM) and ERα/β of biphenyl derivatives (22–61).
Table 3. Structural features, observed and predicted values for ERβ IC50 (nM) and ERα/β of 1-OH-3′-Chloro-biphenyl motif (62–76).
Table 4. Structural features, observed and predicted values for ERβ IC50 (nM) and ERα/β of 1-OH-3′,5′-Dichloro-biphenyl motif (77–83).
Table 5. Structural features, observed and predicted values for ERβ IC50 (nM) and ERα/β of 1-OH-3′-CF3-biphenyl motif (84–92).
Table 6. Structural features, observed and predicted values for ERβ IC50 and ERα/β of 4′-hydroxy-biphenyl-carbaldehyde oxime derivatives (93–110).
Table 7. Structural features, observed and predicted values for ERβ IC50 and ERα/β of Substituted biphenyl derivatives (111–116).
Table S1. Docking score and MM-GBSA output parameter values with experimental affinity and selectivity.
Ligand preparation
Compounds selected for the study were constructed using Maestro in sdf formatCitation13. Subsequently, these structures were optimized by means of OPLS-2005 (default settings) using LigPrepCitation14. LigPrep is a utility of Schrödinger software that combines tools for generating 3D structures from 2D (sdf) representation. LigPrep also assigns an appropriate bond order with correct chiralities for each successfully processed input structure and produces a number of structures from each input structure with various ionization states, tautomers, stereoisomers and ring conformations.
Protein preparation
The X-ray crystal protein structure of ERβ in complex with compound genistein (protein data bank (PDB) ID: 1QKMCitation15) obtained from the RCSB PDB was selected for the present study (resolution: 1.80 Å). Protein structure was prepared using the Maestro software package and aligned using the protein structure alignment module in PrimeCitation16. Bond orders and formal charges were added for heterogroups, and hydrogens were added to all atoms in the system. Protein was inspected visually for accuracy in the χ2 dihedral angle of Asn and His residues and the χ3 angle of Gln, and rotated by 180° when needed to maximize hydrogen bonding. The proper His tautomer was also manually selected to maximize hydrogen bonding. All Asp, Glu, Arg, and Lys residues were left in their charged state. Water molecules of crystallization were removed from the complex except in the active site. A brief relaxation was performed on structure using the protein preparation module in Maestro with the “refinement only” option. This is a two-part procedure that consists of optimizing hydroxyl and thiol torsions in the first stage followed by an all-atom constrained minimization carried out with the impact refinement module (ImprefCitation17) using the OPLS-2005 force field to alleviate steric clashes that may exist in the original PDB structures. The minimization was terminated when the energy converged or the root mean square deviation (RMSD) reached a maximum cut-off of 0.30 Å.
Grid generation and ligand docking
Grids were defined by centering them around the ligand in the crystal structure using the default box size setting in Glide: A scaling by 1.0 was done on the van der Waals radii of protein atoms having partial atomic charge of less then or equal to 0.25. Hydrogen bond constraints were not applied. Flips of 5- and 6-member rings were allowed and non-planar conformation of amide bonds was penalized. van der Waals radii of ligand atoms with partial atomic charge less than 0.15 were scaled by 0.8. The prepared ligands were docked against the ERβ receptor. All docking calculations were performed using the ‘‘extra precision’’ (XP) mode of Glide programCitation18. The obtained minimized poses of docking are rescored using Schrödinger’s proprietary GlideScore (GScore) scoring function. GScore is a modified version of ChemScore, but includes a steric-clash term and adds buried polar terms devised by Schrödinger to penalize electrostatic mismatches.
where vdW = van der Waal’s energy, Coul = Coulomb energy, Lipo = Lipophilic contact term, Hbond = Hydrogen-bonding term, Metal = Metal-binding term, BuryP = Penalty for buried polar groups, RotB = Penalty for freezing rotatable bonds, Site = Polar interactions at the active site; and the coefficients of vdW and Coul are: a = 0.065, b = 0.13018.
After docking, the ligand-receptor complex was subjected to MM-GBSA approach. All computations were carried out with the Linux OS (Red Hat Enterprise WS 5.0).
MM-GBSA
The pose with the lowest GScore of all the 116 ligands was rescored using Prime MM-GBSA approach. This approach is used to predict the free energy of binding for set of ligands to receptor. The docked poses were minimized using the local optimization feature in Prime and the ligand strain energies. Energies of the ligand-receptor complexes were calculated using Prime MM-GBSA technology with all receptor residues being held frozen.
The binding free energy ΔGbind is estimated as
ΔEMM = difference in energy between the complex structure and the sum of the energies of the ligand and unliganded protein using the OPLS force field.
ΔGsolv = difference in the GBSA solvation energy of the complex and the sum of the solvation energies for the ligand and unliganded protein.
ΔGSA = difference in the surface area energy for the complex and the sum of the surface area energies for the ligand and unliganded proteinCitation6.
ArgusLab
ArgusLab 4.0.1 (Mark Thompson and Planaria Software LLC) was used for performing molecular docking experimentsCitation19. The prepared protein and ligands from Maestro, Schrödinger was imported in ArgusLab for docking studies. The active site of the selected protein was defined to include residues within 3.5 Å radius to the complexed ligand. For docking, we have used the ArgusLab scoring function AScore, ArgusDock engine, grid resolution of 0.4 Å with a flexible mode of ligand docking. The docking score was calculated as best ligand-pose energy (kcal/mole) and the score was assessed for correlation between ERβ binding affinity and selectivity.
Correlation between docking, MM-GBSA parameters and affinity/selectivity of ERβ ligands
The docking output parameters from two different docking programs and MM-GBSA output parameters from Schrödinger were exported to StrikeCitation20, which is a collection of chemically aware statistical tools for examining correlations within data. Regression analysis was performed using Strike to provide the correlation models and was measured by the value of the correlation coefficient. In order to explore the reliability of the proposed models, cross validation method was used. The cross validation analysis was performed using the leave multiple out (LMO) method in which five compounds were removed from the data set and their activities predicted using the model derived from the rest of the data points. Significance of the obtained models was evaluated by their respective F values. F is the Fischer ratio which indicates the significance level for a particular multiple linear regression models.
Results and discussion
As evidence suggests that ERβ activation can promote cell differentiation in normal cells and inhibit cell proliferation in normal as well as in cancer cells; recent efforts have been focused on the discovery of ERβ selective ligands targeting cancer cells. Even though the overall tertiary struture of ERβ and ERα proteins remains relatively invariant, the size and shape of their respective ligand binding pockets show considerable variability and hence a range of ligands can able to bind with similar affinities. ERβ receptor is unique in which ligand binding is mediated not only by direct ligand-receptor interactions but also by indirect stabilization of water molecule. ERβ is easier to antagonize than ERα, which further complicates the design of an ERβ selective agonist because agonist induced orientation of helix12 in ERβ is more unstable than in ERαCitation15. In this study, ERβ receptor was taken as the molecular target, prediction model was built by considering the Glide output parameters, and ΔGbind as descriptors to assess ligand’s binding affinity and selectivity.
Glide extra precision docking mode
In the present study, XP mode of Glide custom scoring function was employed to identify ligand poses and ranking. Even though the standard-precision (SP) mode takes less computational time, we used XP mode since the XP scoring function (GScore) includes additional terms over the SP scoring function, and a more complete treatment of some of the SP terms. Further, a large database of co-crystallized structures has been used to optimize the parameters associated with the penalty terms. GScore XP specifically rewards ligand-receptor hydrophobic interactionsCitation18. Since all estrogen ligands bind in a hydrophobic cavity within the core of the receptor, the XP mode will be more suitable for our studies.
Docking of genistein with ERβ in Glide and ArgusLab
The original crystal structure of ERβ-genistein complex (PDB ID: 1QKM) was used to validate the Glide-XP docking protocol. This was carried out by removing the co-crystallized genistein ligand outside the active site and then re-docked in the active site. The RMSD was calculated for the best configuration in comparison to the co-crystallized genistein and it was found to be 1.60 Å. Our docking studies revealed that one of the phenolic hydroxyl groups of the genistein interacts with the side chains of Glu 305, Arg 346 and a buried water molecule in ERβ. The flavone portion of genistein occupies a position similar to that adopted by the C- and D-rings of E2 and the other hydroxyl makes a hydrogen bond with His 475 at the distal end, which is in accordance with the X-ray crystallographic data. The two proximal residues to the bound agonists that differ between ERα and ERβ are the ERα Leu384, which is replaced by Met336 in ERβ, and the ERα Met421 which is replaced by Ile373 in ERβ. These proximal residues are the obvious candidates for the selective effect of genistein and it is explained in detail elsewhereCitation15. The best-docked configuration of genistein in Glide-XP was superimposed with the X-ray crystal 1QKM (see supplementary Figure A). Similar interactions between genistein and ERβ were observed in the docking study using ArgusLab. The RMSD calculated for the best configuration in comparison to the co-crystallized genistein was found to be 1.87 Å. The best-docked configuration of genistein in ArgusLab was also superimposed with the X-ray crystal 1QKM (see supplementary Figure B).
Binding of ligand data set with ERβ
ERβ selective ligand series: 7-substituted 2-phenyl-benzofurans, 4-hydroxy-biphenyl derivatives and 4-hydroxy-biphenyl-carbaldehyde oxime derivatives have been docked in the genistein binding site of ERβ (PDB ID: 1QKM). Then the molecules were subjected to MM-GBSA approach. From the docking results, it is evident that all the molecules have a common binding mode in the binding pocket of ERβ i.e phenolic hydroxyl interacts with the side chains of Glu 305, Arg 346 and a buried water molecule in ERβ and the other O2 hydroxyl binds with the His 475. This may be because all the three nonsteroidal ligand series taken in the present study have used a common design strategy to synthesize molecules containing two hydroxyls with a hydrophobic spacerCitation10–12. In 4-hydroxy-biphenyl-carbaldehyde oximes alone, oxime derivatives extended much further towards and interacted with His 475 of ERβ. The binding modes of all ligands within genistein binding site are superimposed and represented in . From the docking figure, it can be stated that all the ligands are well fitted to the defined binding pocket.
Figure 2. Superimposition of the best-docked configurations and binding modes of 116 ERβ ligands taken in the present study in the genistein binding pocket of ERβ receptor 1QKM.
Building models for prediction of pIC50 using Gscore and Prime/MM-GBSA energy
The three series of ERβ molecules selected for the present study were evaluated in competitive radioligand binding assay for their binding affinity and selectivity. It has been noted that 7-substituted 2-phenyl benzofurans had significantly better activity pIC50 in the range of 7.495–9.456 and selectivity ranges between 0.602–2.033. Biphenyl derivatives showed pIC50 in the range of 4.638–7.619 and selectivity 0.0–1.857, whereas, 4-hydroxy-biphenyl-carbaldehyde oxime derivatives had pIC50 range from 5.070 to 8.096 and selectivity 0.477–1.690.
The docking study in Glide and ArgusLab indicates no significant correlation (r2 = 0.513 and 0.382, respectively) between the docking score and the experimental IC50 for the prediction of activity of the compounds. However, the correlation between Glide docking score (GScore) and experimental IC50 (r2 = 0.513) was better than earlier reports made with different targetsCitation21,Citation22. It has been widely, although not universally, believed among molecular modelers that docking and scoring software are not good enough to rank order compounds with respect to their activity. This poor correlation may be due to the sensitivity of the scoring functions, which can substantially amplify small docking inaccuracies. The rigidity of the protein can also negatively affect the accuracy of these results; hence, incorporation of receptor flexibility could yield more accurate predictions. Scoring functions like Gscore contain simplifications and approximations that limit their ability to fully account for entropic and solvent effects and hence might not correlate well with binding affinitiesCitation23. In addition, the docking score is more useful to measure how well a ligand’s pose complements the binding site rather than correlation with the experimental IC50. Regardless, reasonable correlation was achieved using Glide, Schrödinger in our present studies. However, an enrichment of actives is noted as the score decreases and this agrees with the usefulness of the Glide score for virtual screening.
When compared to docking scoring functions, the MM-GB/SA procedure provided improved enrichment of active ligands and better correlation between calculated binding affinities and experimental data in earlier studiesCitation6,Citation7. However, in the present study, MM-GBSA re-scoring of docked poses does not improve the correlation with the experimental IC50 (r2 = 0.491) when compared to GScore. The less correlation of MM-GBSA re-scoring and experimental IC50 may be related to the fact that in these calculations translational, rotational and vibrational entropy changes are ignored. Glide energy is a specially constructed Coulomb-van der Waals interaction energy score that is formulated to avoid overly rewarding charge-charge interactions and it is also a component of GScore scoring functionCitation18. Plotting of calculated Glide energy versus pIC50 for ERβ (as in ) produced high correlation (r2 = 0.667) between the calculated and experimental values. Inspection of vdw and coul terms, which are van der Waals interaction energy and electrostatic interaction energy predicted by Glide-XP docking, revealed that all ligands share higher values of van der Waals interaction energy, which is logical considering the lipophilic distal part substitutions. In addition, a striking observation was that for all the most active compounds, the van der Waals term is higher than the inactive compounds. In another study, it was stated that placing an aryl or heteroaryl ring at a similar position and orientation as the γ-pyrone of genistein are able to achieve better contact and stronger van der Waals interactions with ERβ Met336. Moreover, Met-aromatic interactions were thought to stabilize the protein. So it is possible that the interaction of the ligand with the ERβ Met336 could lead to increased affinity which was in accord with previous study Citation24. Specifically, ligands that engage in more favourable van der Waals interactions with the protein tend to place nonpolar groups more optimally in to the most hydrophobic pockets. Scoring function like Gscore emphasizes van der Waals and hydrophobic interactions and works better for the binding site with a large nonpolar surface area Citation18, which is in accordance with the present results.
Figure 3. Scatter plot showing the correlation between the values of experimental pIC50 and predicted pIC50 of 116 ERβ ligands using equation [1].
![Figure 3. Scatter plot showing the correlation between the values of experimental pIC50 and predicted pIC50 of 116 ERβ ligands using equation [1].](/cms/asset/e06016fd-f867-4aab-b272-1d3996f02ea5/ienz_a_618990_f0003_b.gif)
N = 116; r2 = 0.667; SD = 0.611; F (1,114) = 227.9; r2CV = 0.655
The root mean square error (RMSE) between the experimental pIC50 values and the predicted pIC50 values obtained by the regression model was 0.616, which is an indicator of the robustness of the fit and the LMO cross validated value (r2CV) suggested that the calculated pIC50 is reliable. The statistical significance of the prediction model was evaluated by the correlation coefficient r2, standard error, F-test value and LMO cross validation coefficient (r2cv).
MM-GBSA out-performed GScore and ArgusLab calculated energy, which did not correlate with experimental selectivity of ERβ ligands data ().
Figure 4. Scatter plot showing the correlation between the values of experimental selectivity and predicted selectivity of ERβ ligands using equation [2].
![Figure 4. Scatter plot showing the correlation between the values of experimental selectivity and predicted selectivity of ERβ ligands using equation [2].](/cms/asset/0e2e086b-ba78-49cc-bf38-6ce4a5ef9415/ienz_a_618990_f0004_b.gif)
N = 116; r2 = 0.736; SD = 0.272; F (1,114) = 319.1; r2CV = 0.725
The obtained RMSE between the experimental and the predicted selectivity using regression model was 0.177. In the present study, MM-GBSA studies were carried out with all receptor residues being held frozen and not allowing flexibility to spherical regions falling within 4, 8, and 12 Å from the ligand. In case of ERβ, the binding process of ligands is mainly dominated by electrostatic effects. For ligands with similar size and shape but differing in their charge distributions, it may be stated that MM-GBSA score provides the correct balance between the ligand desolvation penalty and the electrostatic interactions with the protein. MM-GB/SA describes different contributions to binding mainly by the loss of ligand-solvent interactions, the gain of protein-ligand interactions and the ligand intramolecular strain upon binding with the receptor Citation25. Hence, it is not surprising that MM-GB/SA provides better correlation between calculated binding affinities and experimental data.
The most selective compound in the series is compound 14, which displayed a log selectivity of 2.033. From the docking studies, it has been observed that the 7-acetonitrile derivative of the benzofuran series is in close proximity to ERβ Ile373, which in turn is substituted by Met421 in ERα (). The observed selectivity could be due to the development of a partial negative charge in acetonitrile moiety exerting a repulsive interaction with the electronegative sulphur atom of Met421 in ERα. Other studies also suggest that electronegative and nonpolarizable groups would be more likely to experience a repulsive interaction with ERα Met421 thereby contributing ERβ selectivity Citation9,Citation24. It is to be noted that some functional groups that are not electronegative/nonpolarizable showed modest selectivity. This may be due to differential steric repulsion or due to optimal contact made with ERβ Ile373 but not with ERα Met421. In addition, incorporation of electronegative groups ortho to the 2-phenyl hydroxyl moiety of benzofuran could have increased the selectivity.
Figure 5. The acetonitrile moiety of benzofuran derivative and its best-docked pose with ERβ (PDB ID: 1QKM).
Conclusion
Even though in the present study the experimental affinities and selectivity of ERβ ligands differ in several orders of magnitude, acceptable correlations with the estimated Glide energies and MM-GBSA re-scoring were achieved. Glide docking program from Schrödinger, outperformed ArgusLab and showed better correlation with binding affinities and selectivity of ERβ ligands. The current study found a good correlation between the calculated Glide energy (a component of Gscore function) and the experimental IC50 of ERβ ligands. Compared to affinity, the present study very well explained the selectivity using physics-based scoring by MM-GBSA method and the results were superior to docking. However, Glide energy and MM-GBSA re-scoring results cannot be generalized to predict binding affinity and selectivity of other scaffolds. Inaccuracies in the force field and the ligand/protein induced fit effects are still likely to play important roles when dealing with different scaffolds. Since recent force fields are able to accurately predict the conformational changes in the protein, the ligand/protein induced fit effect is likely to play an important role when dealing with different scaffolds. Further investigation using induced fit docking methodology is required to yield more information. The obtained results will aid the discovery of novel ERβ ligands with isoform selectivity.
Supplementary Material
Download PDF (176.4 KB)Acknowledgment
The authors would also like to thank Dr. Madhavi Sastry and Dr. Shashidhar N. Rao from Schrödinger, India for their valuable suggestions and technical assistance.
Declaration of interest
The authors are grateful for the financial support from the Department of Science and Technology, SERC division India.
References
- Heldring N, Pike A, Andersson S, Matthews J, Cheng G, Hartman J et al. Estrogen receptors: How do they signal and what are their targets. Physiol Rev 2007;87:905–931.
- Buch I, Giorgino T, De Fabritiis G. Complete reconstruction of an enzyme-inhibitor binding process by molecular dynamics simulations. Proc Natl Acad Sci USA 2011;108:10184–10189.
- Genheden S, Nilsson I, Ryde U. Binding affinities of factor Xa inhibitors estimated by thermodynamic integration and MM/GBSA. J Chem Inf Model 2011;51:947–958.
- Gilson MK, Zhou HX. Calculation of protein-ligand binding affinities. Annu Rev Biophys Biomol Struct 2007;36:21–42.
- Fidelak J, Juraszek J, Branduardi D, Bianciotto M, Gervasio FL. Free-energy-based methods for binding profile determination in a congeneric series of CDK2 inhibitors. J Phys Chem B 2010;114:9516–9524.
- Lyne PD, Lamb ML, Saeh JC. Accurate prediction of the relative potencies of members of a series of kinase inhibitors using molecular docking and MM-GBSA scoring. J Med Chem 2006;49:4805–4808.
- Guimarães CR, Cardozo M. MM-GB/SA rescoring of docking poses in structure-based lead optimization. J Chem Inf Model 2008;48:958–970.
- Shen J, Tan C, Zhang Y, Li X, Li W, Huang J et al. Discovery of potent ligands for estrogen receptor β by structure-based virtual screening. J Med Chem 2010;53:5361–5365.
- Tuccinardi T, Bertini S, Martinelli A, Minutolo F, Ortore G, Placanica G et al. Synthesis of anthranylaldoxime derivatives as estrogen receptor ligands and computational prediction of binding modes. J Med Chem 2006;49:5001–5012.
- Collini MD, Kaufman DH, Manas ES, Harris HA, Henderson RA, Xu ZB et al. 7-Substituted 2-phenyl-benzofurans as ERβ selective ligands. Bioorg Med Chem Lett 2004;14:4925–4929.
- Edsall RJ Jr, Harris HA, Manas ES, Mewshaw RE. ERβ ligands. Part 1: The discovery of ERβ selective ligands which embrace the 4-hydroxy-biphenyl template. Bioorg Med Chem 2003;11:3457–3474.
- Yang C, Edsall R Jr, Harris HA, Zhang X, Manas ES, Mewshaw RE. ERβ ligands. Part 2: Synthesis and structure-activity relationships of a series of 4-hydroxy-biphenyl-carbaldehyde oxime derivatives. Bioorg Med Chem 2004;12:2553–2570.
- Maestro, Version 9.0, Schrödinger, LLC, New York, NY, 2009.
- Ligprep, Version 2.3, Schrödinger, LLC, New York, NY, 2009.
- Pike AC, Brzozowski AM, Hubbard RE, Bonn T, Thorsell AG, Engström O et al. Structure of the ligand-binding domain of oestrogen receptor β in the presence of a partial agonist and a full antagonist. EMBO J 1999;18:4608–4618.
- Prime, Version 2.1, Schrödinger, LLC, New York, NY, 2009.
- Impact, Version 5.5, Schrödinger, LLC, New York, NY, 2009.
- Glide, Version 5.5, Schrödinger, LLC, New York, NY, 2009.
- Thompson MA (2004) ArgusLab 4.0.1. Planaria Software LLC, Seattle, WA.
- Strike, 1.8, Schrödinger, LLC, New York, NY, 2009.
- Leach AR, Shoichet BK, Peishoff CE. Prediction of protein-ligand interactions. Docking and scoring: Successes and gaps. J Med Chem 2006;49:5851–5855.
- Tirado-Rives, Jorgensen WL. Contribution of conformer focusing to the uncertainty in predicting free energies for protein-ligand binding. J Med Chem 2006;49:5880–5884.
- Coupez B, Lewis RA. Docking and scoring–theoretically easy, practically impossible? Curr Med Chem 2006;13:2995–3003.
- Hsieh RW, Rajan SS, Sharma SK, Guo Y, DeSombre ER, Mrksich M et al. Identification of ligands with bicyclic scaffolds provides insights into mechanisms of estrogen receptor subtype selectivity. J Biol Chem 2006;281:17909–17919.
- Guimarães CR, Mathiowetz AM. Addressing limitations with the MM-GB/SA scoring procedure using the WaterMap method and free energy perturbation calculations. J Chem Inf Model 2010;50:547–559.