Structure-based virtual screening and molecular docking for the identification of potential multi-targeted inhibitors against breast cancer

Abstract

Breast cancer is characterized by an uncontrolled growth of cells in breast tissue. Genes that foster cell growth in breast cells are overexpressed, giving rise to breast tumors. The identification of effective inhibitors represents a rational chemopreventive strategy. The current in silico study provides a pharmacoinformatic approach for the identification of active compounds against a co-chaperone HSP90 and the human epidermal growth factor receptors EGFR and HER2/neu receptor. The elevated levels of expression of these target proteins have been documented in breast cancer. The utilization of drug-likeness filters helped to evaluate the pharmacological activity of potential lead compounds. Those fulfilling this criterion were subjected to energy minimization for 1000 steepest descent steps at a root means square gradient of 0.02 with an Amber ff12SB force field. Based on molecular docking results and binding interaction analysis, this study represents five chemical compounds (S-258282355, S-258012947, S-259417539, S-258002927, and S-259411474) that indicate high binding energies that range between −8.7 to −10.3 kcal/mol. With high cytochrome P inhibitory promiscuity activity, these multi-targeted potential hits portray not only good physiochemical interactions but also an excellent profile of absorption, distribution, metabolism, excretion, and toxicity, which hypothesizes that these compounds can be developed as anticancer drugs in the near future.

Keywords:

• EGFR
• metastasis
• TNBC
• HER2
• heat shock proteins

Introduction

Breast cancer is characterized by the uncontrolled growth of cells forming a hard painless lump in the breast tissue, usually in the milk ducts or lobules that provide them with milk.^1,2 The most common practice to classify breast tumors is according to the status of three specific cell surface receptors: the estrogen receptor (ER), the progesterone receptor, and the human epidermal growth factor (EGF) receptor HER2/neu receptor. Approximately 75% of all breast cancers are hormone receptor-positive.³ HER2-positive breast cancer accounts for 20–30% of hormone receptor-positive breast cancer that relates to the overexpression of HER2/neu protein. A rare form of breast cancer, triple-negative breast cancer (TNBC), comprises tumors cells that lack receptors for estrogen and progesterone and do not overexpress the HER2 protein.

One of the factors that has been shown to increase a woman’s risk of developing breast cancer is age.⁴ Breast cancer in women has been observed to occur at or over the age of 50. A positive personal or family history is also known to pose an increased risk of developing breast cancer. In addition, a long menstrual life or the use of hormone replacement therapy after menopause increases the risk of developing breast cancer.⁵

Growth factors and their receptors play pivotal roles in the regulation of epithelial cell growth and differentiation.⁶ The HER/erbB family of receptors includes the EGF receptor (EGFR, HER1, and erbB-1), the orphan HER2/neu (erbB-2), and the neuregulin/neuregulin receptors HER3 (erbB-3) and HER4 (erbB-4). Increased knowledge about the pathophysiological mechanisms of breast cancer has led to exponential growth in the identification of biomolecular markers. In this study, a rational approach has been utilized to research the more widely recognized molecular targets in breast cancer. Although there are FDA-approved drugs for breast cancer treatment, their use is accompanied by various side effects. Therefore, it is important to devise an alternative treatment that not only treats the tumor but also has no side effects.

Members of EGF family specifically interact with the EGF receptor (EGFR), which is a cell surface receptor.⁷ Binding of specific ligands such as EGF and transforming growth factor α activates EGFR. PI3 kinase, Ras-Raf-MAPK, JNK, and PLCγ mediate the major EGFR signaling pathways.^8,9 Loss of cell polarization is one of the main features observed when most of the members of the EGFR family are activated, leading to scattering and invasion of breast epithelial cells.¹⁰ Metastasis and angiogenesis processes have been found to be linked with the dysregulation of EGFR pathways. Furthermore, in many human malignancies, this has been observed to be the reason behind the poor prognosis.^11–13 One of the mechanisms of EGFR overexpression is an amplification of the EGFR gene, which has been implicated in oligodendroglioma,¹⁴ glioblastoma, lung cancer,¹⁵ gastric cancer, and breast cancer.¹⁶ In breast cancers, EGFR gene amplification has been observed in 0.8–14% of tumors.^17,18 However, gene amplification has been shown in ~25% of cases of metaplastic breast cancer – a specific phenotype of TNBC.^19,20

The neu gene encodes a protein HER2, similar in structure as human EGFR.^21,22 The HER2 receptor plays a significantly important role in the cell growth and differentiation process, an overexpression of which is associated with the development of human cancers including breast, ovarian, and gastrointestinal (GI) tract cancers.^23,24 HER2 is overexpressed in 15–30% of invasive breast cancers, which has both prognostic and predictive implications. A 40- to 100-fold increase in the expression of HER2 protein results in the tumor cell surface containing ~2 million receptors. Approximately 20–50 HER2-gene copies have been found in breast cancers.²⁵HER2 gene amplification is associated with shorter disease-free and overall survival in breast cancer. The prognostic significance of HER2 amplification has been established in 189 human breast cancers.²⁶ HER2 amplification is one of the early events in human breast tumorigenesis. As per evidence, during progression to invasive disease and nodal or distant metastasis, HER2 status is maintained.²⁷ An increased resistance to certain hormonal agents, enhanced propensity to metastasize to the brain, and high sensitivity to certain cytotoxic chemotherapeutic agents have been found to be associated with breast cancers resulting from an overexpression of HER2.²⁸

Heat shock proteins (HSPs) are members of the molecular chaperones that play an essential role in the folding of cellular proteins.^29,30 Furthermore, through inhibition of programmed cell death and cell senescence during hyperthermia, HSPs directly participate in cell survival.^31–33 In carcinogenesis, HSPs have a role in helping cells (a) to escape tumor suppression pathways, (b) in becoming treatment resistant, (c) in progressing to an advanced stage of cancer, and (d) for facilitating metastasis.³⁴ Tumor suppressor protein (p53), ER, HIF-1alpha, Raf-1 MAP kinase, and a variety of receptor tyrosine kinases, such as HER2, bind with HSP90. These interactions play a key role in breast neoplasia.³⁵ Breast ductal carcinomas indicate an increased HSP90 expression,^36–38 whereas lobular carcinomas and lobular neoplasia manifest significantly decreased HSP90 expression.³⁹ Breast cancer cells become resistant to stress stimuli through HSP90 overexpression. Therefore, therapeutic opportunities in treating cancer can be devised by the pharmacological inhibition of these targets.

In silico approaches have paved the way to solve many biological problems,^40–42 which have led to the identification of novel inhibitors against numerous diseases.^43–47 In this investigation, active compounds against these three targets were screened by structure-based virtual screening (VS) to identify potential virtual hits. The molecular docking tool, AutoDockVina (AD Vina),⁴⁸ was used to dock 50 filtered compounds against EGFR, HER2, and HSP90. The ligands were also analyzed for their profile of absorption, distribution, metabolism, excretion, and toxicity (ADMET), which determined the ADMET efficiency of the drug. Potential hits that indicate good pharmacokinetic (PK) and pharmacodynamic (PD) properties have a better chance to be future drugs. The results of the current study concluded that five multi-targeted compounds with high binding energies as well as a good ADMET profile against all three targets be taken into consideration, suggesting them as potential hits for drug development against breast cancer after testing through in vitro experiments.

Materials and methods

Binding-site analysis

High-resolution crystals with experimentally identified drug sites for EGFR, HSP90, and HER2 were used. Binding pockets of these proteins were examined from crystal structures and were further evaluated using the CASTp server (Computed Atlas of Surface Topography of Proteins).⁴⁹ It locates all likely binding pockets, and the algorithm critically determines the binding pocket and possible cavities in a solvent-accessible surface area.

Protein dataset

X-ray crystallographic structures of EGFR (Protein Data Bank [PDB] ID: 1M17) in complex with the 4-anilinoquinazoline inhibitor erlotinib (resolution: 2.6 Å; R-value free: 0.295), HSP90 complexed with ganetespib (PDB ID: 3TUH; resolution: 1.8 Å; R-value free: 0.213), and HER2 in complex with TAK-285 (PDB ID: 3RCD; resolution: 3.21 Å; R-value free: 0.294) were retrieved from PDB. The criteria for choosing PDBs were (a) minimum resolution and (b) conformation of docked ligand being the same as in the crystallized structure after the redocking procedure.

Protein preparation

To prepare the selected proteins for docking, co-crystallized water molecules, small molecules, nonpolar hydrogens, lone pairs, and nonstandard residues were deleted, and hydrogens and Gasteiger charges were added and merged. Energy minimization and geometry optimization of all structures were performed using Dock Prep, a built-in tool for preparing structures before docking, in UCSF Chimera 10.1 (Resource for Biocomputing, Visualization, and Informatics, University of California San Francisco, San Francisco, CA, USA). The energy was minimized for the 1000 steepest descent steps at a root means square gradient of 0.02, an update interval of 10, and with an Amber ff12SB force field using UCSF Chimera 10.1.⁵⁰ Furthermore, incomplete side chains were repaired using the Dunbrack rotamer library. Hydrogens were added to reasonably generate protonation states at physiological pH. Potentially ambiguous or rare protonation states, especially in binding sites and nonstandard residues, were verified and corrected before charges were added. Charges were added to assign partial charges to atoms. Partial charges were assigned as an atom attribute named charge and included in the Mol2 file (molecule model 3D atomic coordinate files and metadata).

Preparation of chemical library

E-molecule databases and chemical companies including Pub-Chem and Active ZINC were ultilized for literature analysis. All information on chemical compounds was merged, thus obtaining an overall database of ~3 million chemical compounds. Duplicated structures were removed to acquire new scaffolds through InChlKey generated by open babel. Drugs presenting poor PK properties failed to continue progression in the drug development process. The compounds were checked for their oral bioavailability alongside PK parameters such as the blood–brain barrier (BBB), better human intestinal absorption (HIA), good solublity, and low toxicity through the application of Lipinski’s Ro5 (LRo5).⁴⁰ Compounds violating this rule were eliminated from the study. Using Open Bable, the two-dimensional (2D) files of structures that followed this criteria were transformed into three-dimensional (3D) format and saved as molecule files. Discovery Studio^® (DS) visualizer (Accelrys Inc., San Diego, CA, USA)⁵¹ was used for energy minimization during the preparation of ligand molecules.

Virtual screening

VS is a computational technique used to screen potential compounds against specific target protein from chemical compound libraries. The aim was to reduce the number of druggable compounds, inhibiting target proteins more efficiently, to a manageable number. Structure-based virtual screening (SBVS) was utilized to screen substantial databases in the quest of finding novel and potential lead compounds. SBVS estimates the likelihood of a ligand to bind efficiently to the target protein. It involves the application of a scoring function preceeding docking of candidate ligands into a protein target with high binding energies.

In the current study, the drug discovery platform Mcule was utilized to employ SBVS.⁵² For this, an Mcule database of filtered compounds was created manually. Mcule⁴⁸ has a built-in AD Vina tool, which was used to perform molecular docking. In SBVS, molecular docking predicts the binding conformation of small-molecule ligands to the binding site of the selected target. Characterization of the binding behavior plays an essential role in the rational design of drugs as well as to elucidate fundamental biochemical processes. Before VS, an in-depth analysis of the protein structures was performed to understand the binding pocket of target proteins. Using the AutoGrid program of Chimera, a grid map with dimensions of 30 × 30 × 30 Å was generated that covered the binding pocket of each protein. For SBVS, each ligand was docked in the respective binding site of EGFR, HER2, and HSP90, and was scored by maximum binding affinity.

In silico ADMET assessment

Prediction of the ADMET profile for drug candidates and environmental chemicals plays a major role in drug discovery and environmental hazard assessment. Top virtual hits from VS against each protein were further filtered through PK and PD parameters and evaluation of favorable ADMET parameters. To identify the possible adverse effects of these compounds in humans, the ADMET properties of the filtered compounds were predicted using various tools that include the OSIRIS property explorer (Actelion Pharmaceuticals Ltd., Allschwil, Switzerland) which highlights undesired effects if posed by the compound being analyzed, admetSAR (http://lmmd.ecust.edu.cn:8000/) for prediction of quantitative structure–activity relationship-based ADMET properties, and Molsoft (Molsoft LLC, La Jolla, CA, USA) for calculation of drug-likeness properties. ADMET parameters include the BBB, HIA, P-glycoprotein substrates and inhibitors, renal organic cation transporter (ROCT), cytochrome P (CYP) inhibitory promiscuity, and toxicity risks.

Results and discussion

The aim of the current study was to elucidate alternative inhibitory compounds against EGFR, HER2, and HSP90 proteins. In silico VS of ~3 million compounds was performed by establishing a pipeline of LRo5 and PK properties to assess drug likeness. Compounds that showed a strong binding affinity for EGFR, HER2, and HSP90 proteins were selected for further investigations. Further, the ADMET profile was analyzed to filter top hits to identify suitable virtual hit candidate compounds. The current in silico study was undertaken to identify efficient anti-breast cancer compounds. The five multi-targeted top compounds were also evaluated for their PK and PD properties. The virtual hits identified in this study can be used as an alternative targeting agent for breast cancer after being tested through in vitro experiments.

Binding-site evaluation

The crystal structure of the kinase domain of the EGFR in complex with 4-anilinoquinazoline inhibitor erlotinib (PDB ID: 1M17) was obtained from PDB. It was found that residues Leu694, Leu768, Met769, Gly772, and Leu820 formed significant interactions with the inhibitor. The crystal structure of HER2 complexed with TAK-285 (PDB ID: 3RCD) indicated Ala751, Leu800, Met801, Leu852, and Asp863 to be among the significant interacting residues. The X-ray co-crystallized structure of TAK-285 with HER2 demonstrated that it interacts with the expected residues in the respective ATP pocket. HSP90 complexed with ganetespib (PDB ID: 3TUH) indicated that the residues Asn51, Lys58, Asp93, Gly97, and Thr184 form significant interactions at the binding site. The binding-site residues of all three target proteins were also confirmed by using CASTp. Binding-site interacting residues, as reported in the PDB, are shown in Figure 1.

Figure 1 (A) Epidermal growth factor receptor (EGFR) tyrosine kinase domain with 4-anilinoquinazoline inhibitor erlotinib. (B) Novel HER2/EGFR dual. Inhibitors bearing a pyrrolo[3,2-d] pyrimidine scaffold. (C) The crystal structure of the N-terminal domain of an HSP90 in the presence of an inhibitor ganetespib.

Notes: Green dashed lines indicate hydrogen bonds with distance in Angstrom (Å), spoked red arcs indicate hydrophobic contacts, atoms are shown in black for carbon, blue for nitrogen, red represents oxygen, yellow represents sulfur, and green represents fluorine.

Structure-based virtual screening

To identify anti-breast cancer compounds, a database of ~3 million compounds was screened by the application of a series of filters (LRo5, drug likeness, PK filters), which led us to select the top 66, 832 best compounds. These compounds were subjected to SBVS against EGFR, HER2, and HSP90 separately to identify multi-targeted compounds common among all these targets. The flowchart of multistep VS is shown in Figure 2. The docking compounds were ranked based on their highest binding energy with the corresponding protein. To identify potential hits having binding affinities greater than the reported co-crystallized inhibitor, binding energy cutoffs of −8.9, −8.5, and −8.3 kcal/mol were set against EGFR, HER2, and HSP90, respectively, after redocking with the reported co-crystallized inhibitor as deposited in the PDB to select only those hits which were above that specific threshold. Only top-ranked hits with higher binding energy for all three proteins were considered and investigated further.

Figure 2 Schematic workflow summarizing the screening of active compounds against breast cancer.

Abbreviations: AD Vina, AutoDockVina; PD, pharmacodynamic; PDB, Protein Data Bank; PK, pharmacokinetic.

Binding energy cutoff value was applied for VS to select only the high-ranked compounds. A total of 96 hits with binding energies lower than −8.9 kcal/mol for EGFR, 114 hits with <−8.5 kcal/mol for HER2, and 79 hits with <−8.3 kcal/mol for HSP90 were chosen based upon the AD Vina docking score (∆G). Visual inspection was carried out to screen these 289 compounds further. As a result, 71 compounds with promising high binding energies against all three targets were selected. An extensive ADMET analysis was performed to screen these virtual hits, and their molecular interactions with binding-site residues were carefully analyzed.

Prediction of ADMET properties

Molecular descriptors are the deciding factor for PK properties and toxicity of a compound. ADMET properties predicted in silico identify the likelihood of compounds to be used as human therapeutic agent.⁵³ The admetSAR online server was utilized for calculating the ADMET properties of 71 common compounds of all three proteins. It is important for compounds to have a promising ADMET profile. The BBB,⁵⁴ HIA,⁵⁵ aqueous solubility,⁵⁶ Caco-2 cell permeability, CYP450 inhibition,⁵⁷ and Ames toxicity were calculated for 71 compounds. It was observed that only 15 virtual hits were successful at passing through these ADMET filters. These 15 virtual hits were further filtered through a series of drug safety profiling parameters. Undesirable moieties and substructures causing potential toxicity were eliminated through a series of PAINS (pan assay interference compounds) filters. Among the 15, only five compounds were unidentified by PAINS-1, -2, and -3.

The admetSAR predictions indicated that the ability to penetrate the BBB and HIA was portrayed by all compounds; in addition, all the compounds were revealed as non-inhibitors of the P-gp inhibitor except S-258002927. None of the compounds indicated any inhibitory effects on the ROCT. CYP enzymes, including various CYP450 substrates and inhibitors, play a fundamental role in drug metabolism. The results showed that most compounds were non-inhibitors of CYP450 enzymes. A significant finding was observed for compounds S-258012947 and S-258002927 that indicated low CYP-inhibitory promiscuity, being non-inhibitors of all CYP450 enzymes including 1A2, 2C9, 2D6, 2C19, and 3A4. Further drug metabolism analyses showed these compounds were non-substrates of two CYP450 substrates (2C9 and 2D6). The ADMET properties of the top five best virtual hits as common against all three targets are tabulated in Table 1, and their 2D structures are presented in Figure 3.

Figure 3 Chemical structures of lead compounds: (A) S-258002927; (B) S-258012947; (C) S-258282355; (D) S-259411474; (E) S-259417539.

Table 1 ADMET profile of potential five multi-targeted virtual hits

ADMET profile	S-258282355	S-258012947	S-259417539	S-258002927	S-259411474
BBB penetration	+	+	+	+	+
HIA	+	+	+	+	+
Caco-2 permeabilty	+	−	−	+	−
Aqueous solubility	−3.9498	−2.1264	−3.7602	−4.5343	−3.7813
P-gp
Substrate	−	−	+	−	−
Inhibitor	−	−	−	+	+
CYP450 substrate
CYP450 2C9	−	−	−	−	−
CYP450 2D6	−	−	−	−	−
CYP450 3A4	+	+	+	+	+
CYP450 inhibitor
CYP450 1A2 inhibitor	−	−	+	−	−
CYP450 2C9 inhibitor	−	−	+	−	−
CYP450 2D6 inhibitor	+	−	−	−	−
CYP450 2C19 inhibitor	−	−	+	−	−
CYP450 3A4 inhibitor	−	−	−	−	−
CYP IP	High	Low	High	Low	High
ROCT	−	−	−	−	−

Abbreviations: ADMET, absorption, distribution, metabolism, excretion, and toxicity; BBB, blood–brain barrier; HIA, human intestinal absorption; CYP, cytochrome P; IP, inhibitory promiscuity; ROCT, renal organic cation transporter; +, present; −, not present.

Molecular docking studies of EGFR

All five compounds showed high binding energies with EGFR ranging from −8.9 to −9.7 kcal/mol with AD Vina, as presented in Table 2. All compounds established a network of molecular interactions (H-bonds, Van der Waals [VdW], alkyl, π-alkyl, and π-sigma bonds) with the active-site residues of EGFR when analyzed in 2D plots as shown in Figure 4. S-259417539 formed various binding interactions including H-bonds (Thr766), π-alkyl (Ala719), alkyl (Val702), and VdW (Leu694, Cys751, Asp831, Lys721, Leu764, Met742, Glu738, Asn818, Arg817). S-259411474 established H-bond and alkyl bond with Met769 and Leu820, π-alkyl with Ala719 and Val702, and VdW interactions with Ile720, Ile765, Leu768, Leu764, Leu753, Thr830, Asp831, Gly695, Glu738, Cys773, Phe699, and Gly772, respectively. Likewise, S-258012947 was involved in various binary interactions including H-bond (Asp831), π-alkyl (Leu694), and VdW (Ile720, Ile765, Cys773, Leu834, Gly833, and Thr765), whereas S-258282355 and S-258002927 were found to interact with Asn818 and Asp831 through H-bond, with residues Leu694, Leu820, and Val702 by forming alkyl bonds, and with residues Arg817, Thr766, Gln767, Leu768, Met769, Gly772, and Cys773, and Gly695, Gly772, Arg817, Cys773, Leu820, Ala719, Met769, Leu768, Phe699, Thr830, Leu764, Glu738, Gln767, and Met742 through VdW interactions, respectively.

Figure 4 Molecular surface representation of EGFR with top-ranked multi-targeted virtual hits: (A) S-258002927; (B) S-258012947; (C) S-258282355; (D) S-259411474; (E) S-259417539.

Notes: Molecular surface representation of the EGFR binding pocket (in blue) with respective ligands in stick format (in purple). Alongside each 3D complex, 2D interaction plots indicate important binding-site interactions between respective ligands and binding-site residues.

Table 2 Binding energy and molecular interaction of top five multi-targeted compounds with EGFR

Compound molecule ID	Binding energy (kcal/mol)	Van der Waals	H-bond	π-alkyl	Alkyl	π-sigma
S-258282355	−9.7	Arg817, Thr766, Gln767, Leu768, Met769, Gly772, Cys773	Asn818	Leu694, Leu820	Lys721, Val702, Ala719	–
S-258012947	−9.3	Ile720, Leu764, Ile765, Thr766, Gly772, Cys773, Lys721, Leu834, Gly833, Gln767, Leu768, Gly695	Asp831	Ala719, Leu820, Met769	Leu694	Val702
S-259417539	−9.1	Leu694, Cys751, Asp831, Lys721, Leu764, Met742, Glu738, Asn818, Arg817	Thr766	Ala719	Cys773	Leu820
S-259411474	−8.9	Ile720, Ile765, Leu768, Leu764, Leu753, Thr830, Asp831, Gly695, Glu738, Cys773, Phe699, Gly772	Met769	Ala719, Val702	Leu820	Leu694
S-258002927	−9	Gly695, Gly772, Arg817, Cys773, Leu820, Ala719, Met769, Leu768, Phe699, Thr830, Leu764, Glu738, Gln767, Met742	Asp831	Val702	–	Leu694

Molecular docking studies of HER2

The top hits represented high binding energies with HER2 ranging from −9.4 to −10.3 kcal/mol, as shown in Table 3. Post-docking analysis of HER2 revealed that all five compounds (S-259417539, S-259411474, S-258012947, S-258002927, and S-258282355) were also found to bind inside the kinase domain, surrounded by interacting residues with noticeably high binding energy value that ranged from −9.4 to −10.3 kcal/mol (Figure 5). The ligand S-258282355 was found not to form any H-bond, but alkyl bonds with Leu694 and Leu820 were observed as were π-alkyl bonds with Lys721, Val702, and Ala719 residues, whereas VdW with Ile752, Leu726, Leu800, Thr798, Ser783, Arg784, Thr862, Leu852, Asp863, Phe864, Leu755, Gly865, and Glu770 were found. The S-258012947 established H-bonds (Asp863), π-alkyl (Leu694), alkyl (Ala719, Leu820, and Met769), and VdW (Asn850, Phe731, Arg849, Gly729, Phe864, Gly727, Ser783, Cys805, Gly804, Leu800, Met801, Thr798, Leu785, Val797, and Ile752). Similarly, H-bond for S-258002927 was found to form H-bonds with Asp863 and Cys805, π-alkyl with Ala719, π-sulfur with Met774, alkyl bonds with Leu764 and Lys721, and VdW with Ile752, Thr798, Thr862, Ser783, Arg784, Gln799, Val734, Met801, Leu800, Leu726, Phe1004, and Gly804. S-259411474 and S-259417539 were found to form H-bonds with Thr862, Thr798, Leu785, Ser783, and Asp863, and Lys753, respectively, π-alkyl with Val702, Ala719, and Cys773, alkyl with Leu820, Val702, Ala719, and VdW with Met801, Gly804, Leu800, Ile752, Val734, Asp863, Lys753, Val797, and Asn850, and Phe731, and Arg849, Gly729, Phe864, Gly727, Ser783, Cys805, Gly804, Leu800, Met801, Thr798, Leu785, Val797, and Ile752 residues, respectively.

Figure 5 Molecular surface representation of HER2 with respective ligands: (A) S-258002927; (B) S-258012947; (C) S-258282355; (D) S-259411474; (E) S-259417539.

Notes: Molecular surface representation of HER2 binding pocket (in blue) with respective ligands in stick format (in purple). Alongside each 3D complex, 2D interaction plots indicate important binding-site interactions between respective ligands and binding-site residues.

Table 3 Binding energy and molecular interaction of top five multi-targeted compounds with HER2

Compound Mcule ID	Binding energy (kcal/mol)	Van der Waals	H-bond	π-sulfur	π-alkyl	Alkyl	π-sigma
S-258282355	−10.3	Ile752, Leu726, Leu800, Thr798, Ser783, Arg784, Thr862, Leu852, Asp863, Phe864, Leu755, Gly865, Glu770	–	–	Lys721, Val702, Ala719	Leu694, Leu820	Leu785, Lys753, Leu796
S-259411474	−9.8	Met801, Gly804, Leu800, Ile752, Val734, Asp863, Lys753, Val797	Thr862, Thr798, Leu785, Ser783	–	Ala719, Val702	Leu820	Leu852
S-258002927	−9.5	Ile752, Thr798, Thr862, Ser783, Arg784, Gln799, Val734, Met801, Leu800, Leu726, Phe1004, Gly804	Asp863, Cys805	Met774	Ala719	Leu764, Lys721	Leu785, Leu852
S-259417539	−9.4	Gly804, Phe1004, Cys805, Asn850, Leu755, Leu726, Phe731, Gly727, Thr862, Leu785, Thr798	Asp863, Lys753	–	Cys773	Val702, Ala719	Leu820
S-258012947	−9.4	Asn850, Phe731, Arg849, Gly729, Phe864, Gly727, Ser783, Cys805, Gly804, Leu800, Met801, Thr798, Leu785, Val797, Ile752	Asp863	–	Leu694	Ala719, Leu820, Met769	Val702

Molecular docking studies of HSP90

The binding energies of top hits with HSP90 ranged from −8.7 to −9.4 kcal/mol, as represented in Table 4, and binary interactions are displayed in Figure 6. The binding interactions of S-258002927 with HSP90 included H-bond with Thr184, π-alkyl with Ala55, π-sulfur with Met98, and VdW with residues Val186, Leu48, Val150, Ser52, Leu107, Phe138, Asn51, Val136, Tyr139, Gly135, Gly97, Ile96, Lys58, Asp54, Asn106, Lys112, Ala111, and Ser113. H-bond for S-258012947 was found with Ala111, alkyl bond with Ala111 residue, and VdW with Ile91, Ile26, Val186, Ala55, Thr184, Met98, Ser113, Tyr139, Gly135, Val150, Asn106, and Leu107 residues. S-258282355 did not form any H-bonds; π-alkyl bond was formed with residues Met98 and Ala55, alkyl with Val136 and Ala111, π-sulfur with Met98, and VdW with Thr109, Ile110, Lys112, Leu48, Ser113, Ile26, Leu107, Tyr139, Gly137, Val186, Phe138, Thr184, Asn106, Lys58, Ile96, and Gly97. S-259411474 and S-259417539 were observed to form H-bonds with Asn51 and Lys58, respectively, π-alkyl with Ala111, alkyl bond with Leu107, Phe138, Val186, Leu48, and Ala55, π-sulfur with Met98, and VdW with Ile26, Thr109, Tyr139, Val150, Val136, Asn106, Gly135, Ser113, Ile110, Lys112, Lys58, Gly97, Ser52, Asp93, Ile96 and Val186, and Ser52, Ala55, Met98, Val150, Asn51, Asp54, Thr184, Asp93, Phe138, Leu107, Gly137, Val136, Tyr139, Ser113, Lys112, and Ile110, respectively.

Figure 6 Molecular surface representation of HSP90 with respective ligands: (A) S-258002927; (B) S-258012947; (C) S-258282355; (D) S-259411474; (E) S-259417539.

Notes: Molecular surface representation of HSP90 binding pocket (in blue) with respective ligands in stick format (in purple). Alongside each 3D complex, 2D interaction plots indicate important binding-site interactions between respective ligands and binding-site residues.

Table 4 Binding energy and molecular interaction of top five multi-targeted compounds with HSP90

Compound Mcule ID	Binding energy (kcal/mol)	Van der Waals	H-bond	π-sulfur	π-alkyl	Alkyl	π-sigma
S-259411474	−9.4	Ile26, Thr109, Tyr139, Val150, Val136, Asn106, Gly135, Ser113, Ile110, Lys112, Lys58, Gly97, Ser52, Asp93, Ile96	Asn51	Met98	Ala111	Leu107, Phe138, Val186, Leu48, Ala55	Thr184
S-259417539	−9.2	Val186, Ser52, Ala55, Met98, Val150, Asn51, Asp54, Thr184, Asp93, Phe138, Leu107, Gly137, Val136, Tyr139, Ser113, Lys112, Ile110	Lys58	–	Ala111	–	–
S-258002927	−9.2	Val186, Leu48, Val150, Ser52, Leu107, Phe138, Asn51, Val136, Tyr139, Gly135, Gly97, Ile96, Lys58, Asp54, Asn106, Lys112, Ala111, Ser113	Thr184	Met98	Ala55	–	Thr184
S-258282355	−9.1	Thr109, Ile110, Lys112, Leu48, Ser113, Ile26, Leu107, Tyr139, Gly137, Val186, Phe138, Thr184, Asn106, Lys58, Ile96, Gly97	–	Met98	Met98, Ala55	Val136, Ala111	Ala55
S-258012947	−8.7	Ile91, Ile26, Val186, Ala55, Thr184, Met98, Ser113, Tyr139, Gly135, Val150, Asn106, Leu107	Ser52, Asp93, Thr109, Ile110	Ala111	–	Ala111	–

The current study is an effort to identify anti-breast cancer compounds that may be considered for drug development to treat breast cancer. The SBVS of ~3 million compounds against three EGFR, HER2, and HSP90 proteins has assisted in robust screenings of most potential virtual hits. Based on the experimentally determined EGFR, HSP90, and HER2 inhibitors from previous studies, binding energy cutoff to select for possible 96, 79, and 114 virtual hits that could act on EGFR, HSP0, and HER2 proteins, respectively, was set. These virtual hits followed: (a) LRo5 and drug likeness; the compounds violating more than one drug-like parameter indicated that they might lead to issues with bioavailability, and hence, such compounds were not included in the study; (b) good HIA and BBB for moderately lipophilic compounds to cross the BBB, whereas polar molecules are poor central nervous system drugs; and (c) comprehensive ADMET profile. The best possible orientations forming stable ligand–target protein complexes, through the process of molecular docking, were achieved.

The binding energy of the ligand and its half-life quantifies the efficiency of the ligand–protein complex. The rate of dissociation is rapid for weaker interactions; however, the rate of dissociation slows down with an increase in the value of binding energy.⁵⁸ Compounds with a longer half-life have strong binding energies, and hence, they take a longer time to dissociate. Strong interactions between residues imply that binding with inhibitory compounds may be stable, leading to an inhibitory reaction. The top five best virtual hits in complex with EGFR, HER2, and HSP0 revealed strong binding affinities and highlighted several H-bonds and hydrophobic interactions between functional groups, and side chains of essential residues. Post-docking analysis of EGFR, HER2, and HSP90 with these hits disclosed more in-depth details. The EGFR that complexed with S-259429764 showed the highest binding energy (−10.7 kcal/mol). Other compounds showed binding energy in the range from −10.7 to −8.7 kcal/mol. Moreover, common interacting residues that displayed VdW and H-bond with different compounds were Gln767, Leu768, Gly772, Cys773, Leu764, Gly695, and Thr830, respectively (Table 2). However, the HER2 complexed with S-258282355 showed the highest binding energy (−10.3 kcal/mol), and other compounds showed binding energy in the range from −9.8 to −8.9 kcal/mol. Together with this, the common interacting residues that displayed VdW and H-bond with different compounds were Ile752, Leu726, Leu800, Thr798, Ser783, Thr862, Met801, and Asp863, respectively (Table 3). Likewise, HSP90 complexed with S-259411474 showed the highest binding energy (−9.4 kcal/mol) and displayed H-bonds with Asn106 and Ser50. Other compounds showed binding energy in the range from −9.4 to −8.3 kcal/mol. Interestingly, all compounds interacted equally in terms of VdW interactions with Tyr139 and Ser113 of HSP90. Other common interacting residues were Val150, Asn106, Lys112, Val136, Gly135, Ile110, Lys58, Ser52, and Gly97 that displayed VdW interactions with different compounds (Table 4).

Furthermore, all screened virtual hits predicted promising ADMET profile, and their strong binding affinity stipulated the multi-targeted potential of these top virtual hits.

Conclusion

Breast cancer is one of the leading cancer types in women across the globe. Computational tools have been widely used for drug development and discovery of multi-targeted inhibitors of many overexpressed proteins induced in breast cancer. This study reveals five multi-targeted compounds that possess high binding energies against most common target proteins that are involved in breast cancer. With excellent PK and PD properties as predicted, these virtual hits may be considered for early drug development against breast cancer after being tested through in vitro and in vivo studies.

Disclosure

The authors report no conflicts of interest in this work.