Identification of Novel Natural Dual HDAC and Hsp90 Inhibitors for Metastatic TNBC Using e-Pharmacophore Modeling, Molecular Docking, and Molecular Dynamics Studies

Breast cancer (BC) is one of the main types of cancer that endangers women’s lives. The characteristics of triple-negative breast cancer (TNBC) include a high rate of recurrence and the capacity for metastasis; therefore, new therapies are urgently needed to combat TNBC. Dual targeting HDAC6 and Hsp90 has shown good synergistic effects in treating metastatic TNBC. The goal of this study was to find potential HDAC6 and Hsp90 dual inhibitors. Therefore, several in silico approaches have been used. An e-pharmacophore model generation based on the HDAC6-ligand complex and subsequently a pharmacophore-based virtual screening on 270,450 natural compounds from the ZINC were performed, which resulted in 12,663 compounds that corresponded to the obtained pharmacophoric hypothesis. These compounds were docked into HDAC6 and Hsp90. This resulted in the identification of three compounds with good docking scores and favorable free binding energy against the two targets. The top three compounds, namely ZINC000096116556, ZINC000020761262, and ZINC000217668954, were further subjected to ADME prediction and molecular dynamic simulations, which showed promising results in terms of pharmacokinetic properties and stability. As a result, these three compounds can be considered potential HDAC6 and Hsp90 dual inhibitors and are recommended for experimental evaluation.


Introduction
Breast cancer (BC) is the most common malignancy in women and the second leading cause of cancer mortality in women [1]. Despite the significant developments in cancer treatment, the American Cancer Society (ACS) predicted 284,200 new BC cases and 44,130 deaths in 2021 [2]. In terms of clinical and biological behavior, BC is a very complex, heterogeneous disease with multiple categories [3,4]. It can be classified according to the expression of proteins and genes into Three categories: hormone receptor-positive (HR+), human epidermal growth factor receptor 2 (HER2)-enriched, and triple-negative breast cancer (TNBC) [5]. TNBC is an aggressive and dangerous of all BC types. TNBC is characterized by the absence of the progesterone receptor (PR), estrogen receptor (ER), and HER2 expression and amplification, and it accounts for approximately 15% of all BC cases [6,7]. Compared to all other subsets of BC, TNBC malignancies have a much higher Figure 1 summarizes the study's workflow.

Results and Discussion
a good way to prevent BC metastasis [36].
Given the importance of HDAC6 and Hsp90 in TNBC metastasis and the ongoing development of breast cancer-targeting drugs, this study was conducted to find novel dual HDAC6 and Hsp90 inhibitors using e-pharmacophore modeling, screening, molecular docking, and MM-GBSA calculations to calculate the binding affinity towards the two targets. In addition, the results were augmented using molecular dynamics (MD) simulations to investigate the stability of the compound-protein complexes. Figure 1 summarizes the study's workflow.

E-Pharmacophore Modeling and Screening
Numerous studies have demonstrated that inhibiting HDAC6 and Hsp90 simultaneously has a synergistic anticancer effect [37][38][39]. This study aimed to identify natural compounds targeting the HDAC6 and Hsp90 proteins using an in silico approach. A pharmacophoric hypothesis containing seven features (one donor (D), one acceptor (A), negative ionic, and four aromatic rings) was obtained using the Phase module of Maestro ( Figure  2). These features were used to test a library comprising 270,540 natural chemicals against the pharmacophoric hypothesis. Of those natural chemicals, 13,629 compounds matched five out of seven pharmacophoric hypotheses and were thus selected for the upcoming molecular docking.

E-Pharmacophore Modeling and Screening
Numerous studies have demonstrated that inhibiting HDAC6 and Hsp90 simultaneously has a synergistic anticancer effect [37][38][39]. This study aimed to identify natural compounds targeting the HDAC6 and Hsp90 proteins using an in silico approach. A pharmacophoric hypothesis containing seven features (one donor (D), one acceptor (A), negative ionic, and four aromatic rings) was obtained using the Phase module of Maestro ( Figure 2). These features were used to test a library comprising 270,540 natural chemicals against the pharmacophoric hypothesis. Of those natural chemicals, 13,629 compounds matched five out of seven pharmacophoric hypotheses and were thus selected for the upcoming molecular docking.

Molecular Docking and MM-GBSA
Throughout the docking process, the 13,629 compounds were subjected to the HTVS mode of Glide against the HDAC6 protein. A total of 64 compounds showed docking scores better than the reference (<−8.782 Kcal/mol) and were then docked into Hsp90 through the XP docking mode. Thirteen compounds resulting from this XP had docking scores of ≤−6 Kcal/mol. The Prime module's MM-GBSA calculation was applied to the three compounds ( Figure 3) with the highest XP docking scores to determine the free binding energies to HDAC6 and Hsp90, as shown in Table 1. The three compounds showed favorable binding energies ranging from −42.45 to −23.31 Kcal/mol and from −52.69 to −24.66 Kcal/mol for HDAC6 and Hsp90, respectively. Co-crystallized ligands bound with HDAC6 protein showed lower docking scores of −6.14 kcal/mol. The top three compounds with satisfactory free binding energies and good docking scores for the two targets were further investigated, including visual examination of intermolecular interactions with HDAC6 (Table 2 and Figures 4 and 5) and Hsp90 (Table 2 and Figures 6 and 7).

Molecular Docking and MM-GBSA
Throughout the docking process, the 13,629 compounds were subjected to the HTVS mode of Glide against the HDAC6 protein. A total of 64 compounds showed docking scores better than the reference (<−8.782 Kcal/mol) and were then docked into Hsp90 through the XP docking mode. Thirteen compounds resulting from this XP had docking scores of ≤−6 Kcal/mol. The Prime module's MM-GBSA calculation was applied to the three compounds ( Figure 3) with the highest XP docking scores to determine the free binding energies to HDAC6 and Hsp90, as shown in Table 1           HDAC6 consists of a cap, a zinc-binding group (ZBG), and a linker that connects the first two components [40]. It belongs to the metal-protein family in which the stability of the catalytic center is crucially maintained by the zinc ion, and furthermore, the inhibitors of HDAC6 interact by chelating with this ion [41,42]. In this regard, compound A showed a metal coordination bond between its carbonyl oxygen, near the amide group, and the positively charged zinc ion. Moreover, the HDAC6-compound A complex showed three pi-cation interactions with HIE614, HIE463, and PHE583 residues and H-bonds with HIE614 residue (Figures 4A and 5A). Furthermore, compound A displayed a hydrophilic interaction through a water bridge with HIE614. Compound A also exhibited hydrophobic interaction through CYS584, TYR745, PHE583, PRO571, PRO464, LEU712, and PHE643 residues. Similar binding interactions were observed with other inhibitors [43][44][45][46], among which was belinostat, an FDA-approved HDAC6 inhibitor [43,44]. Further investigations on belinostat has been performed by calculating its docking score which found to be −9.511 which is lower than that of compound A (−9.968). Moreover, binding free energy (BFE) has been computed using MM-GBSA which revealed that compound A has BFE of −42.45 Kcal/mol which is higher than the reference belinostat (with a BFE value of −32.77 Kcal/mol).
Benzamide is a ZBG that is important in determining the potency of HDAC6 inhibitors [40]. Compound B contains the benzamide group as a ZBG; it showed interaction through the carbonyl group. Compound B showed three types of interactions with the binding pocket of HDAC6 which are two H-bonds with GLY582 and TYR745 residues; pipi interactions through PHE643, PHE583, and HIE614; and hydrophobic interactions with the PHE643, TYR745, PRO571, PHE583, CYS584, PRO464, and LEU712 amino residues (Figure 4B and 5B). HDAC6 binding affinity and van der Waals interactions are enhanced by pointing to the Phe643 and Leu 712 residues in the so-called L1-loop pocket, as studied by Bai [47].
The carbonyl group of the ZBG benzamide of compound C interacts with the zinc ion, which contributes to the stability of the complex. Moreover, compound C and HDAC6 interactions include three pi-pi interactions through PHE583, while hydrogen bonds were displayed with HIE614, LEU712, and GLY582 residues ( Figures 4C and 5C). It also showed hydrophilic interactions through water bridges with HIE614 and LEU712. ZINC000096116556, ZINC000020761262, and ZINC000217668954, labeled A-B-C, respectively, were the top-ranked compounds that achieved the highest docking scores in the library of natural products of the ZINC database with Hsp90 and satisfactory docking scores with HDAC6 ( Figure 3). The interaction patterns between compounds A-C, HDAC6, and Hsp90 are presented in Figures 4-7.
HDAC6 consists of a cap, a zinc-binding group (ZBG), and a linker that connects the first two components [40]. It belongs to the metal-protein family in which the stability of the catalytic center is crucially maintained by the zinc ion, and furthermore, the inhibitors of HDAC6 interact by chelating with this ion [41,42]. In this regard, compound A showed a metal coordination bond between its carbonyl oxygen, near the amide group, and the positively charged zinc ion. Moreover, the HDAC6-compound A complex showed three pication interactions with HIE614, HIE463, and PHE583 residues and H-bonds with HIE614 residue ( Figures 4A and 5A). Furthermore, compound A displayed a hydrophilic interaction through a water bridge with HIE614. Compound A also exhibited hydrophobic interaction through CYS584, TYR745, PHE583, PRO571, PRO464, LEU712, and PHE643 residues. Similar binding interactions were observed with other inhibitors [43][44][45][46], among which was belinostat, an FDA-approved HDAC6 inhibitor [43,44]. Further investigations on belinostat has been performed by calculating its docking score which found to be −9.511 which is lower than that of compound A (−9.968). Moreover, binding free energy (BFE) has been computed using MM-GBSA which revealed that compound A has BFE of −42.45 Kcal/mol which is higher than the reference belinostat (with a BFE value of −32.77 Kcal/mol).
Benzamide is a ZBG that is important in determining the potency of HDAC6 inhibitors [40]. Compound B contains the benzamide group as a ZBG; it showed interaction through the carbonyl group. Compound B showed three types of interactions with the binding pocket of HDAC6 which are two H-bonds with GLY582 and TYR745 residues; pi-pi interactions through PHE643, PHE583, and HIE614; and hydrophobic interactions with the PHE643, TYR745, PRO571, PHE583, CYS584, PRO464, and LEU712 amino residues ( Figures 4B and 5B). HDAC6 binding affinity and van der Waals interactions are enhanced by pointing to the Phe643 and Leu 712 residues in the so-called L1-loop pocket, as studied by Bai [47].
The carbonyl group of the ZBG benzamide of compound C interacts with the zinc ion, which contributes to the stability of the complex. Moreover, compound C and HDAC6 interactions include three pi-pi interactions through PHE583, while hydrogen bonds were displayed with HIE614, LEU712, and GLY582 residues ( Figures 4C and 5C). It also showed hydrophilic interactions through water bridges with HIE614 and LEU712. Leu712 is thought to be crucial for inhibitory actions according to Nguyen's research [44]. Furthermore, hydrophobic interactions through ALA641, PRO464, CYS584, TYR745, PHE643, PHE583, PRO571, PHE642, and LEU712 residues were noticed. Two residues (PHE583 and LEU712) were reported by H. Losson et al. and were considered important residues for the inhibition of HDAC6, which were also with compound C [43].
It is essential to emphasize that the HDAC6 catalytic domain, where the acetylated lysine is deacetylated, contains a typical narrow hydrophobic channel formed by residues PRO464, GLY582, PHE583, PHE643, and LEU712 [48]. Interactions with one or more residues in the catalytic site were observed with compounds A, B, and C. An interaction with PHE643 was reported with R-trichostatin (TSA), a histone deacetylase inhibitor [49]. In this regard, compound B formed a hydrogen bond with TYR745 which emphasized bond stabilizes the reaction intermediate to eventually release the product with deacetylated lysine residues. This interaction was observed with the bioactive conformer of TSA, but it was absent in the inactive one [49].
On the other hand, the three compounds interacted with Hsp90's essential amino acid residues. For example, the Hsp90 compound A complex conveyed a pi-cation interaction through PHE138 and TYR139. PHE138 was found to be involved in the inhibition of Hsp90 by the Sanchez group [50]. Furthermore, it formed a hydrophilic interaction through a water bridge with ASN51. Asn51 was discovered to be an essential residue to interact stably with Hsp90, as reported by Magwenyane [51]. The hydrophobic interactions were noticed with ALA111, LEU107, LEU103, VAL150, TYR139, PHE138, TRP162, VAL136, VAL186, MET98, ILE96, ALA55, and LEU48 residues ( Figure 6A and 7A). M. Abbasi et al., who reported these findings, highlighted the significance of interactions with PHE138, ALA55, LEU107, VAL186, TYR139, and VAL136, in addition to other amino acid residues [52,53].
Many interactions were formed between compound C and the Hsp90 binding site ( Figures 6C and 7C). Pi-pi interactions with PHE138 and TYR139 residues which have been revealed to play key role in their complex formation were also observed by Rezvani et al. [55]. Moreover, LYS58 and TYR139 made H-bond interactions, whereas TRP162, LEU103, TYR139, PHE138, LEU107, VAL150, ALA55, MET98, ILE96, VAL186, VAL136, and ALA111 residues were observed in hydrophobic interactions. Additionally, compound C showed hydrophilic interactions through water bridges with LYS58 and TYR139, and these two residues were also reported by Jia [56]. PHE138 had been documented by El-Shafey et al. [57] as one of the residues crucial to Hsp90's active binding site's recognition of their reference ligand, geldanamycin.
The top three compounds, ZINC000096116556, ZINC000020761262, and ZINC000217668954, were further subjected to ADME prediction and MD simulations.

ADMET Prediction
The pharmacokinetic profile of the effective clinical candidates comprises its ADMET properties (absorption, distribution, metabolism, elimination, and toxicity), which are crucial for assessing their pharmacodynamics effect; therefore, these molecules must pass the ADMET test to be considered drug-like [54]. In this study, ADMET analysis was conducted on the three best-found hits with higher molecular docking and free binding energy using the QikProp module of Maestro (Table 3). According to their ADMET profile, these compounds passed "the rule of 5" and other pharmacokinetic tests with a (0) value, indicating that they were drug-like in agreement with their ADMET properties. Pharmacokinetic parameters such as QPPCaco and QPlogBB show the likelihood of transporting a molecule through the gut and the brain-blood barrier, respectively, thus allowing the prediction of the compound's permeation through several barrier models. ZINC000096116556 has a greater gastrointestinal absorption, with a QPPCaco value of 1480.139, compared to ZINC000217668954, ZINC000020761262, and the reference, which have values of 105.569, 21.522, and 100.046, respectively. Limited CNS penetration is predicted by low QPlogBB values of (−0.618, −1.797, −2.276) for ZINC000096116556, ZINC000217668954, and ZINC000020761262, respectively, which are almost in the same range as the reference value of (−2.123). These top three hits have conformational independent aqueous values of (CIQplogs;  Regarding Percent Human-Oral Absorption, none of the top three compounds deviates from the recommended values (0 to 100%). ZINC000096116556 has a greater HPOA of (100%) in comparison with ZINC000217668954, ZINC000020761262, and the reference HPOA, which have values of 71.807, 56.718, and 69.793, respectively. This property fully agrees with human oral absorption values since both measure the same parameter. None of them had a low HOA value (HOA values coded 1, 2, or 3 for low, medium, and high, respectively). Furthermore, our findings revealed that the QPlog HERG value of ZINC000020761262 was much lower than the other compounds and the reference compound. This demonstrated that the reference compound could be more cardiotoxic than our hits (Table 3).

Molecular Dynamic (MD)
In order to assess the flexibility and stability of the protein HDAC-Ligands (ZINC000096116556, ZINC000020761262, ZINC000217668954) complexes, molecular dynamics simulations for 200 ns were carried out using the Desmond software [56]. Such MD results included the root mean square deviation (RMSD), which was evaluated to confirm the stability of the protein system during the simulation period by measuring the variable distance between atoms where a distance range of 1-3 Å is generally acceptable for small globular proteins [57,58]. The HDAC-ZINC000096116556 complex fluctuated throughout the simulation with a maximum peak of less than 3 Å; however, it remained lower than the other ligands and even the protein backbone until it showed sudden stability in~25 ns and remained bound until~100; after that, the ligand showed significant fluctuations until the end of the simulation time. ZINC000020761262 and ZINC000217668954 exhibited similar fluctuation patterns, since they showed deviations until~100 ns and maintained steadiness until 200 ns, with little deviations as observed in Figure 8. The RMSF (root mean square fluctuation) of each residue can be used to assess w residues are involved in protein or complex changes in structure. Low RMSF residu more stable since the motions of the residue atoms during modeling are limited, and vice [59,60], ZINC000096116556 bound to HDAC6 has a notably low RMSF average value o Å , indicating its interaction stability as illustrated in (Figure 9). By con The RMSF (root mean square fluctuation) of each residue can be used to assess which residues are involved in protein or complex changes in structure. Low RMSF residues are more stable since the motions of the residue atoms during modeling are limited, and vice versa [59,60], ZINC000096116556 bound to HDAC6 has a notably low RMSF average value of 2.04 Å, indicating its interaction stability as illustrated in (Figure 9). By contrast, ZINC000020761262 and ZINC000217668954 have relatively high values of 4.89 and 3.38 Å respectively.  Generally, protein-ligand interactions are facilitated by hydrogen bonds, and their strength in a water environment (Figure 10 and Figure S1 Supplementary Materials) provides further provisions for the docking results for HDAC6 and ZINC000096116556, ZINC000020761262, ZINC000217668954 interactions. ZINC000096116556 formed strong hydrophobic interactions with PHE643 (130%), HIS574 (~48%), and PHE642 (~38%). It also had hydrogen bonds with TYR745 (~100%) and water bridge interactions with HIS615 (78%) and ASP612 (88%). ZINC000020761262 interacted through water bridge bonds with HIS462 (25%); in addition, it formed hydrophobic interactions with PHE533 (48%), HIS462 (30%), PHE583 (19%), and HIS463 (11%) and formed hydrogen bonds with HIS463 (53%), HIS462 (23%), and HIS574 (10%). ZINC000217668954 primarily interacted with PHE642 (>50%) through hydrophobic bonds and formed hydrogen bonds with GLU638 (13%), ALA641 (38%), and HIS462 (21%). Through water bridge bonds, meanwhile, it formed water bridge interactions with HIS614 (35%), GLU638 (22%), HIS462 (13%), and ALA641 (<10%). These interactions confirm the docking results, since they show the same type of Generally, protein-ligand interactions are facilitated by hydrogen bonds, and their strength in a water environment (Figure 10 and Figure S1 Supplementary Materials) provides further provisions for the docking results for HDAC6 and ZINC000096116556, ZINC000020761262, ZINC000217668954 interactions. ZINC000096116556 formed strong hydrophobic interactions with PHE643 (130%), HIS574 (~48%), and PHE642 (~38%). It also had hydrogen bonds with TYR745 (~100%) and water bridge interactions with HIS615 (78%) and ASP612 (88%). ZINC000020761262 interacted through water bridge bonds with HIS462 (25%); in addition, it formed hydrophobic interactions with PHE533 (48%), HIS462 (30%), PHE583 (19%), and HIS463 (11%) and formed hydrogen bonds with HIS463 (53%), HIS462 (23%), and HIS574 (10%). ZINC000217668954 primarily interacted with PHE642 (>50%) through hydrophobic bonds and formed hydrogen bonds with GLU638 (13%), ALA641 (38%), and HIS462 (21%). Through water bridge bonds, meanwhile, it formed water bridge interactions with HIS614 (35%), GLU638 (22%), HIS462 (13%), and ALA641 (<10%). These interactions confirm the docking results, since they show the same type of interactions with (PHE 643) and (PHE583) through hydrophobic bonds and with ALA641 through the water bridge bond.  As shown in Figure S2 (Supplementary Materials), the calculated ligand parameters for the top three compounds included RMSD, which for the three ligands stabilized at 0.9659, 2.406, and 1.701 A˚, respectively, throughout the simulation period, demonstrating the better stability of ZINC000096116556 compared to the other ligands. Moreover, the stabilities of the complexes were examined using the radius of gyrus (Rg) measurement to determine the compactness. The most stable ligands were those with the lowest Rg value. It can be seen from Table 4 that the average Rg values obtained from the complexes HDAC-ZINC000096116556 complex, HDAC-ZINC000020761262 complex, and HDAC-ZINC000217668954 complex were about 5.630, 5.990, and 5.294 Å , respectively, indicating that the HDAC-ZINC000217668954 complex has best compactness compared to the other complexes. Meanwhile, the SASA analysis was carried out to complete the stability analysis of each hit ligand complex, and the identification of SASA was conducted to anticipate the protein conformational changes that water molecules may access throughout the simulation. Commonly, the stability of the ligand-receptor complex increases with decreasing SASA values. As shown in Table 4, the HDAC-ZINC000096116556 complex has the lowest SASA average value, i.e., 259 Å , compared to the other two hits. The molecular surface (MolSA) is equivalent to the Van der Waals surface zone, and the three compounds had values of 358.5, 429.5, and 553.1 Å 2 , respectively. Intramolecular H-bonds display the number of intramolecular hydrogen bonds in ligand atoms. The three compounds showed no significant intramolecular H-bonds, according to molecular dynamics data. Furthermore, each protein-ligand's polar surface area (PSA), which is defined as the total number of polar atoms on all molecule surfaces, was measured. Less occupied space makes the complex more stable. ZINC000096116556 had a PSA average of 130 Å 2 . The As shown in Figure S2 (Supplementary Materials), the calculated ligand parameters for the top three compounds included RMSD, which for the three ligands stabilized at 0.9659, 2.406, and 1.701 Å, respectively, throughout the simulation period, demonstrating the better stability of ZINC000096116556 compared to the other ligands. Moreover, the stabilities of the complexes were examined using the radius of gyrus (Rg) measurement to determine the compactness. The most stable ligands were those with the lowest Rg value. It can be seen from Table 4 that the average Rg values obtained from the complexes HDAC-ZINC000096116556 complex, HDAC-ZINC000020761262 complex, and HDAC-ZINC000217668954 complex were about 5.630, 5.990, and 5.294 Å, respectively, indicating that the HDAC-ZINC000217668954 complex has best compactness compared to the other complexes. Meanwhile, the SASA analysis was carried out to complete the stability analysis of each hit ligand complex, and the identification of SASA was conducted to anticipate the protein conformational changes that water molecules may access throughout the simulation. Commonly, the stability of the ligand-receptor complex increases with decreasing SASA values. As shown in Table 4, the HDAC-ZINC000096116556 complex has the lowest SASA average value, i.e., 259 Å, compared to the other two hits. The molecular surface (MolSA) is equivalent to the Van der Waals surface zone, and the three compounds had values of 358.5, 429.5, and 553.1 Å 2 , respectively. Intramolecular H-bonds display the number of intramolecular hydrogen bonds in ligand atoms. The three compounds showed no significant intramolecular H-bonds, according to molecular dynamics data. Furthermore, each protein-ligand's polar surface area (PSA), which is defined as the total number of polar atoms on all molecule surfaces, was measured. Less occupied space makes the complex more stable. ZINC000096116556 had a PSA average of 130 Å 2 . The simula-tion process for the three molecules also revealed the torsional trend of rotatable bonds ( Figure S3, Supplementary File). As a result, each rotatable bond in these three compounds maintains its orientation during a simulation that lasts 200 nanoseconds, proving again the stability of their conformation.

Materials and Methods
In silico studies were conducted in Maestro v12.8 of the Schrödinger suite. Molecular dynamics simulations were performed using academic Desmond v6.5 by D.E. Shaw Research.

Preparation of Proteins and Ligands
The crystallographic structures of HDAC6 and Hsp90 with PDB IDs 6PYE and 3D0B, respectively, were retrieved from the protein data bank (https://www.rcsb.org/ accessed on 25 May 2022). After that, with the Protein Preparation Wizard tool of Maestro [61], the two proteins' structures were prepared for the upcoming procedures. Firstly, they underwent structural pre-processing where bond orders were assigned, absent hydrogen atoms were added, zero-order bonds to metals and disulfide bonds were created, incomplete side chains and loops were filled, water molecules beyond 5 Å were deleted, and het states were generated at 7 ± 2 pH via Epik tool. Secondly, the pre-processed structures underwent refinement. In this step, the crystallized water molecules' orientations were assigned and the residues' protonation states were checked and corrected via PROPKA. Finally, the refined structures were passed through a restrained minimization process under the OPLS4 force field, in which heavy atoms converged to an RMSD of 0.30 Å.
A library containing 270,540 natural product structures was obtained from the ZINC database (https://zinc.docking.org/ accessed on 25 May 2022), then, using the MacroModel tool with the default parameters [62], the library was subjected to energy minimization under the OPLS4 force field to ensure that the compounds' structures were acceptable for further computational calculations.

E-Pharmacophore Modeling and Virtual Screening
In the Phase module, the E-Pharmacophore generation panel was opened and the option "create a pharmacophore model using receptor-ligand complex" was selected to create a pharmacophore model from the HDAC6-ligand complex [63]. Subsequently, the prepared natural products' library was screened against the generated pharmacophoric features via the Phase Ligand Screening panel, resulting in a group of compounds that matched the features [64].

Grid Generation and Molecular Docking
One of the main approaches in computational chemistry is molecular docking. Schrödinger's Glide module includes three docking methodologies: HTVS (high-throughput virtual screening), SP (standard precision), and XP (extra precision). These three methodologies differ in terms of speed and accuracy [65,66].
To carry out molecular docking, the Receptor Grid Generation panel was used to create the grid boxes of the two proteins where the docking occurs [67]. Since the structures from which grids will be created are receptors with co-crystallized ligands, the co-crystallized ligands molecules were identified to be distinguished from the receptor. The van der Waals radius scaling settings were left as default: scaling factor = 1 and partial charge cutoff = 0.25. Then, the receptor grid generation process was allowed to run. Thereafter, the molecular docking was conducted using the Ligand Docking panel of the Glide module with the default settings (only the docking methodology was changed each time). The group of compounds obtained from the screening against the predicted pharmacophoric features was docked to the HDAC6-generated grid through the HTVS mode of Glide. The compounds that achieved high docking scores were further docked via the XP mode. In the end, the resulting compounds with high docking scores on HDAC6 were docked to the Hsp90 grid through the XP mode.

Free Binding Energy Calculations
The MM-GBSA panel of the Prime module was used to estimate the free binding energies of the top compounds according to their docking scores with HDAC6 and Hsp90. [68]. The free binding energy of the receptor-ligand complexes is determined by the following equation: where ∆E is the free binding energy, E c is the target/ligand complex energy, E R is the receptor energy, and E L is the ligand energy [61]. The solvation model was set to be VSGB and the force field was OPLS4.

ADMET Analysis
The top compounds according to the docking score were also analyzed to predict their pharmacokinetic properties by the QikProp tool of Schrödinger [62]. This initial prediction may be useful in decreasing the failure probability in the future stages of drug development.

Molecular Dynamics (MD) Simulations
After ADME analysis, the top compounds bound to HDAC6 were subjected to the MD simulations utilizing the Desmond software to further investigate the binding interactions and stability of the receptor-ligand complexes [69]. Before starting the simulation process, it is necessary to install a suitable biological system by using the system builder panel, in which the compounds were immersed in 11663 TIP3P molecules in an orthorhombic figured box with diameters of 10 × 10 × 10 Å. Salt was added in a specific concentration of Na + and Cl − charge: 59.239 mM (Total charge + 38) for Na + and 51.445 mM (Total charge-33) for Cl − . Then, the system underwent energy minimization using OPLS4. The system was equilibrated in a two-phase sequential equilibration simulation: isothermal-isochoric (NVT) and isothermalisobaric (NPT) ensembles. After that, the simulation process was allowed to run for 200 ns using an NPT ensemble at a temperature of 1 K and pressure of 1 bar which was kept constant throughout the simulation duration. The Nose-Hoover chain and the Martyna-Tobias-Klein methods were used as a thermostat and a barostat, respectively. The cutoff radius = 9.0 Å. The trajectory of the recording interval was every 100 p. A total of 2000 frames were obtained during the simulation. The Simulation Interaction Diagram tool of Desmond was used to analyze the resulting frames, such as RMSD, RMSF, and receptor-ligand contacts. The RMSD for a frame X is calculated by the following equation (Equation (1)): where N is the number of atoms in the atom selection; t ref is the reference time (typically the first frame is used as the reference and it is regarded as time t = 0); and r is the position of the selected atoms in a frame x after superimposing on the reference frame, where frame x is recorded at time t x . The procedure is repeated for every frame in the simulation trajectory. The RMSF for a residue i is calculated by the following formula (Equation (2)): where T is the trajectory time over which the RMSF is calculated; r i is the position of residue I; r is the position of atoms in residue i after superposition on the reference; and the angle brackets indicate that the average of the squared distance is taken over the selection of atoms in the residue.

Conclusions
This study aimed to identify new dual inhibitors to combat one of the most aggressive types of breast cancer, TNBC, by considering two therapeutic targets: HDAC6 and Hsp90. Several computational chemistry methods, including e-pharmacophore modeling, molecular docking, MM-GBSA calculations, ADMET analysis, and molecular dynamics, were applied to screen 270,450 natural compounds from the ZINC database for their dual inhibition capability on the two targets. In this study, e-pharmacophore-based virtual screening followed by molecular docking have identified three compounds (ZINC000096116556, ZINC000020761262, and ZINC000217668954) that exhibited the highest docking scores with Hsp90 and satisfactory docking scores with HDAC6. In addition, the investigation of the interaction patterns of the three compounds showed many favorable interactions that agree with previous studies in the literature. Among these interactions, the chelation of the zinc ion in the HDAC6 active site plays a major role in HDAC6 inhibition. Further analysis of the three compounds revealed that they have acceptable free binding energies, favorable ADMET properties, and good interaction stabilities. These three hits may be potential inhibitors against HDAC6 and Hsp90, and if experimentally examined, they may serve as promising hits against TNBC in the near future.

Data Availability Statement:
The datasets generated during and/or analyzed during the current study are available from the corresponding author upon reasonable request. Schrödinger suite is a commercial software (www.schrodinger.com). Academic Desmond by D.E. Shaw Research is available at https://www.deshawresearch.com/, accessed on 8 April 2022.