Multitarget-Based Virtual Screening for Identification of Herbal Substances toward Potential Osteoclastic Targets

Abstract

Osteoporosis is a complex bone disease indicating porous bone with low bone mass density and fragility. Cathepsin K, V-ATPase, and α_Vβ₃ integrin are exhibited as novel targets for osteoporosis treatment. Our preliminary study uses a state-of-the-art method, including target-based virtual screening and clustering methods to determine promising candidates with multitarget properties. Phytochemicals with osteoprotective properties from the literature are used to elucidate the molecular interactions toward three targets. The binding scores of compounds are normalized and rescored. The K-means and hierarchical clustering methods are applied to filter and define the promising compounds, and the silhouette analysis is supposed to validate the clustering method. We explore 108 herbal compounds by virtual screening and the cluster approach, and find that rutin, sagittatoside A, icariin, and kaempferitrin showed strong binding affinities against Cathepsin K, V-ATPase, and α_Vβ₃ integrin. Dockings of candidates toward three targets also provide the protein-ligand interactions and crucial amino acids for binding. Our study provides a straightforward and less time-consuming approach to exploring the new multitarget candidates for further investigations, using a combination of in silico methods.

1. Introduction

Osteoporosis is the most common bone disorder arising from an imbalance between bone formation by osteoblasts and bone resorption by osteoclasts [1]. Reducing the rate of bone degradation becomes the therapeutic target for osteoporosis treatment due to its short duration, within 2–4 weeks. Osteoclasts are multinucleated cells derived from hematopoietic stem cells in the bone marrow to form preosteoclast precursors, and are activated by various factors, such as inflammatory cytokines, chemokines [2,3].

Focusing on osteoclast function, the pre-osteoclasts from hematopoietic stem cells (MSCs) adhere to the bone mineral matrix using the cell adhesion integrin receptor [4,5]. The α_Vβ₃ type of integrin receptor is highly expressed in osteoclasts. The inhibition of the α_Vβ₃ integrin receptor resulted in decreasing bone resorption [5]. After osteoclast activation, the hydrogen ions were transported to the bone-resorbing compartment via the vacuolar-type H+-ATPase (V-ATPase) [6]. V-ATPase is a proton pump responsible for the regulation of intracellular and extracellular pH environments. V-ATPase shows the capability to maintain the acidic environment of various organelles, such as lysosomes, secretory granules. The acidifying compartment and extracellular environment are the predominant roles of V-ATPase, leading to bone resorption and tumor metastasis [6]. Recently, many evidence-based studies proved the correlation of V-ATPase function and various diseases, especially cancer and osteoporosis [7]. In bone resorption, secretory lysosomes with V-ATPase in their membrane fuse with the specialized membrane of osteoclasts and then release the degrading enzymes and proton ion into the bone-resorbing compartment [8,9]. Cathepsin K is a lysosomal protease enzyme that mainly degrades type I collagen and elastin in the bone matrix at a low pH environment [10,11]. The inhibition of Cathepsin K suppresses the osteoclast function without impairing its viability, indicating the benefits over current medicines [12]. The current medicines used in the treatment of osteoporosis are bisphosphonates, such as alendronate, which are targeted in the reduction of the bone resorption rate. However, bisphosphonates exhibited several side effects and interrupted the bone remodeling process by induced cell death in osteoclasts. The alternative medicines for osteoporosis treatment are herbal medicines, especially traditional Chinese medicines (TCMs).

The use of medicinal plants and natural products in traditional and modern therapy are gaining more interest nowadays. Many studies exhibited that medicinal plants exert their therapeutic potential due to a synergistic effect from a combination of phytochemicals [13,14,15]. It was suggested that natural products tend to possess multitarget activity, and less toxicity than synthetic compounds and therefore are potentially a useful starting point for multitarget drug discovery. Additionally, the use of current medicines in osteoporosis treatment may be limited with the treatment period. Polyphenols were extensively investigated, according to their wide variety of efficacies on human health, including antioxidative, anti-inflammatory, anticancer properties, and prevention of cardiovascular diseases, neurodegenerative, and osteoporosis [16,17,18]. Polyphenols, specifically flavonoids, possess a positive impact on bone health. The possible mechanisms of polyphenols to balance bone metabolism are due to their anti-oxidative and anti-inflammatory properties, relevant to the number of phenol groups [19].

Multitarget drug discovery is a new paradigm in drug design and development [20]. The multitarget approach has been increasingly interesting for complex disease drug development recently, particularly cancer and Alzheimer’s disease [21]. Conceptually, multitarget agents possessing modest activity at more clinically relevant targets may produce similar or better in vivo bioactivities than high target-selective compounds [22]. Several in silico methods were applied to find multitarget candidates. Virtual screening is an in silico method used to randomly screen a vast amount of conformations from a large chemical library to identify the high probability candidates bound to the target of interest [23]. Target-based virtual screening is carried out by docking small molecules to the 3D structure of the targeted site, followed by optimizing binding configurations and evaluating binding feasibility based on molecular interactions. However, identifying promising candidates from many docked poses generated from virtual screening is still critical for successful drug discovery. Clustering is an unsupervised machine learning method for partitioning the dataset into a set of groups or clusters. Ideally, the similarity within a same cluster should be higher than distinct clusters. Several clustering approaches are extensively used in the field of molecular modeling, especially k-means [24,25] and hierarchical algorithm [26]. However, each clustering method has its limitations, such as sensitivity to outliers, and identifying the number of clusters in advance. To avoid these limitations, we used a combination of these methods for analyzing the molecular results and visualizing the relationships between top predictive features.

Our study aimed to find the multitarget candidates from the herbal substances that reported the osteoprotective effects. To find an insight into the osteoclastic activity at a molecular level, we focused on studying proteins involving bone resorption in the different stages of osteoclasts, including Cathepsin K, V-ATPase, and α_Vβ₃ integrin receptor, by performing the virtual screening. The top 26 hits from each protein–ligand docking were filtered by several clustering methods to find the promising candidates for these three target proteins. The final candidates were defined from our proposed method. Our method would provide a straightforward and reasonable strategy to identify the promising multitarget ligands for further use in other investigations by being less time-consuming. An overview of the study for screening the multitarget compounds with anti-osteoclastic activity was shown in Figure 1. The phytochemicals with the osteoprotective activity from the previous literature [27,28] were selected to be used in this work. The in silico methods, such as a molecular docking, and a cluster approach were utilized to explore the promising multitarget candidates for further studies.

2. Materials and Methods

2.1. Target Selection and Protein-Protein Interaction Analysis

The targets used in this study were selected based on the specificity and necessity for the osteoclast function. With regard to the stages of osteoclast function, α_Vβ₃ integrin receptor was a potential target in the polarizing stage, whereas V-ATPase and Cathepsin K were presented as the crucial targets in the resorbing stage. After that, we analyzed the connectivity of these proteins through the protein–protein interaction (PPI) network analysis.

We performed the PPI network analysis using the STRING (Search Tool for the Retrieval of Interacting Genes) database [29] by assigning the gene targets of particular proteins, including Cathepsin K (CTSK), V-ATPase (V0 domain: ATP6V0A1, ATP6V0A2, TCIRG1, ATP6V0A4, ATP6V0B, ATP6V0C, ATP6V0D, ATP6V0D2, ATP6V0E1, ATP6V0E2; V1 domain: ATP6V1A, ATP6V1B1, ATP6V1B2, ATP6V1C1, ATP6V1C2, ATP6V1D, ATP6V1E1, ATP6V1E2, ATP6V1F, ATP6V1G1, ATP6V1G2, ATP6V1G3, and ATP6V1H), and α_Vβ₃ integrin receptor (ITGAV, ITGB3) with the species defined to “Homo sapiens”, and a confidence cut-off of 0.4. STRING was utilized to search for the neighboring interactions and then generate the PPI network based on the interactions occurring from the input proteins. The PPI analysis was visualized by Cytoscape software [30].

2.2. Data Collection and Ligand Preparation

A total of 108 herbal substances with the osteoprotective activity from previous literatures [27,28] were selected for use in this work, including flavones, flavonoids, phenylpropanol, terpenoids, flavanones, isoflavones, and miscellaneous compounds. The three-dimensional structures of selected herbal substances were downloaded from the PubChem database (https://pubchem.ncbi.nlm.nih.gov/, (accessed on 14 January 2021)). All ligands were geometrically optimized using the Gaussian09w [31] program with the B3LYP method at 6–31 G(d,p) levels [32]. The hydrogen atoms were assigned. Aromatic carbons and rotatable bonds of each ligand were identified through AutoDockTools 1.5.6 program [33].

2.3. Target Preparation and Binding Site Identification

The three-dimensional structures of all targets were retrieved from the Protein Data Bank (RCSB database). The structures with a high resolution and domain completeness were selected in this study. All native ligands, water molecules, and metal ions were removed. Hydrogen atoms and Gasteiger partial charge were assigned by AutoDockTools 1.5.6. The binding site of the native ligand was defined as the ligand-binding site for the molecular docking study.

In this study, a PDB structure of human Cathepsin K was 4 × 6 H with a resolution of 1.0 Å [34]. Prior to performing the molecular docking study, the protonation state of protein was adjusted by PROPKA to provide a structure of amino acids at pH 5.5 and then checked by Biovia Discovery Studio Visualizer [35]. Subsequently, AutoGrid 4.0 was used to define the binding site by creating a grid box with a size of 60 × 60 × 60 ų. A grid center was set at 17.552, 8.018, −19.098 at XYZ-dimension with a grid spacing of 0.375.

The protein structure of V-ATPase used in this study was 5KNB, which cocrystallized with ADP ligands [36]. This structure is a V1 domain which contains subunit A and B. The active site of ADP is located between the two subunits, A and B. A grid box size of 60 × 60 × 60 ų centered at −11.539, 23.975, −18.976 with a grid spacing of 0.375 Å was set as the binding site of this docking set.

The α_vβ₃ integrin receptor used in this experiment was a PDB code of 4G1M [37]. A center of the native ligand was set as a grid box center at dimensions −19.145, 44.449, and 43.944. The grid box size of 80 × 80 × 80 ų with a grid spacing of 0.375 Å was set as the parameters for this study.

2.4. Target-Based Virtual Screening

In all experiments, the genetic algorithm (GA) search method was applied to determine the best pose of each ligand within the active site of the receptor. All ligands were docked into the defined binding site. Top 100 docked conformations generated from 300 population sizes were performed by the Lamarckian genetic algorithm (LGA). Simulations were performed with a maximum of 2,500,000 evaluations and a maximum of 250,000 generations. Overall conformations were ranked by order of binding energy. The redocking of native ligand with the root–mean–square deviation (RMSD) tolerance within 2.0 was used to validate the docking method. The docked conformations with the lowest binding score were selected as a model for cluster analysis.

2.5. Cluster and Subcluster Analyses

The clustering analyses of a data set were performed using RStudio version 3.6.3 [38]. This experiment used the k-means and hierarchical methods to classify the clusters. The first clustering method is a k-means algorithm based on the clustering through nearest distances within a clustered group and the largest distance between different cluster groups. The second algorithm is a hierarchical clustering algorithm, using an algorithm of average distances, which generates a dendrogram that offers the possibility to explore sequentially the groups that demonstrate the best candidates with the high binding affinities against three receptors. Lastly, a combination of hierarchical and k-means called hybrid hierarchical k-means (hkmeans) was used in this study. This approach showed advantages over the two methods above due to not necessarily defining the optimal number of clusters in advance. Additionally, a hierarchical method enables the interactive visualization groups through a heatmap incorporated with the dendrogram. The heatmap plot with the implemented gradient color bar efficiently identified the clusters with greater affinities toward targets and determined which clusters hold a greater number of observations. Interestingly, the dendrogram showed detailed analyses of the hierarchical clustering, which also helped to differentiate each cluster and conformation according to the docking results. The optimal numbers of clusters were defined from the Nbclust function [39], which was implemented in RStudio. This function was used to identify the best number of clusters and suggest the best clustering scheme.

The clusters were validated using the silhouette method [40] and demonstrated in terms of an average silhouette width. Further, the existing data from the cluster analysis were performed the principal component analyses (PCA) in order to show the isolated cluster. The color coding and dendrograms in the clustered heatmap were calculated from the hierarchical clustering algorithms.

2.6. Data Interpretation and Visualization

All structural proteins and the binding modes were generated from Pymol [41] and Discovery Studio Visualizer 2017R [35]. The binding affinities of compounds were presented in the box plot. The clustering heatmap was done in RStudio [38], using the pheatmap [42], ComplexHeatmap [43] and Factoextra [44] packages. The LigandScout version 4.4 [45] was used to generate the 2D and 3D structure-based interactions of ligands toward an isolated target.

3. Results

3.1. Protein–Protein Interaction (PPI) Network

The PPI network was constructed with 25 proteins (Figure 2a). As defined by color, the CTSK gene is mostly found in bone tissues represented in dark blue color, related to the highly expressed in osteoclast from bones. The PPI network using STRING v11.0 showed molecular cross talk from all genes in this study. Evidence suggesting a functional link is made from various data, particularly in co-expression, co-mentioned in PubMed data, and the curated database. We first constructed the PPI network from 25 nodes, 237 edges with an average node degree of 19, and an average local clustering coefficient of 0.936. The PPI enrichment p-value is less than 1 × 10⁻¹⁶, suggesting that these proteins have more interactions among themselves than the expectation for a random set of proteins. Such an enrichment indicates that the proteins are at least partially biologically connected as a group. Further, the PPI network was reconstructed by collecting the nodes showing the intramolecular edges with the STRING score cut-off of 0.5 (Figure 2b). This network was generated from particular nodes, such as CTSK, ITGAV, ITGB3, TCIRG1, ATP6V0D2, and ATP6V1D, with a p-value of 4.78 × 10⁻¹¹. Based on the interaction score calculated from the data mining and co-expression data, CTSK showed strong connectivity with ATP6V0D2, followed by TCIRG1. These two genes belong to the V-ATPase family. CTSK also showed the connectivity with ITGAV and ITGB3 genes, which belong to the α_Vβ₃ integrin receptor. Additionally, ITGAV also exhibited strong connectivity to ATP6V1D. According to these results, it was suggested that this PPI network could be considered small relative to the random one, and the proteins might be biologically relevant.

3.2. Target-Based Virtual Screening Using AutoDock

Traditional Chinese medicines (TCMs) areused in the treatment of several diseases. The polyphenols possess activities that inhibit the adhesion of the osteoclast to the bone matrix, affect the proton transportation relevant to a less acidic environment, and inhibit the degrading enzymes, which should be considered multitarget candidates.

All 108 selected compounds were docked into the active site of three targets: Cathepsin K, V-ATPase, and α_Vβ₃ integrin receptor. The binding energies of docked compounds were presented in terms of AutoDock score. A set of docking scores toward the individual enzyme showed no difference among groups in a range of −3.0 to −14.5 kcal/mol. The validations of the docking set were made by redocking the native ligands, and the rmsd tolerance of native ligands was less than 2.0 Å.

Considering the binding affinities of compounds against the individual target, the top 20 docked scores of compounds are depicted in Figure 3. Kushenol F was used as a positive control for Cathepsin K with a dock score of −8.3 kcal/mol. According to the result, rutin exhibited the highest binding affinity at −13.75 kcal/mol, followed by icariin, kaempferitrin, acteoside, and salvialic acid with the binding affinities of −11.93, −11.35, −11.23, and −11.22 kcal/mol, respectively.

As a docking result of V-ATPase, the compounds were docked into the binding site of ATP (or ADP), as mentioned in the Material and Method section. Diphyllin was used as a positive control for V-ATPase [46] with a dock score of −8.56 kcal/mol. Among the top five candidates, kaempferitrin showed the highest binding affinity of −14.00 kcal/mol toward V-ATPase, followed by rutin, sagittatoside A, and icariin with the binding affinities of −13.91, −13.67, and −13.26 kcal/mol, respectively.

The top 20 binding energies of compounds bound to an active site of α_Vβ₃ integrin receptor are presented in Figure 3. A binding site is defined according to the location of the native ligand. Artemisinin was used as a positive control for α_Vβ₃ integrin receptor [47] with a dock score of −6.41 kcal/mol. The docking result indicated that rutin exhibited the strongest binding affinity of −13.76 kcal/mol toward an integrin receptor, followed by sagittatoside A, icariin, kaempferitrin, and oleuropein with the binding affinities of −12.93, −11.76, −11.57, and −11.21 kcal/mol, respectively.

We rearranged and combined these three docking results from Cathepsin K, V-ATPase, and α_Vβ₃ integrin receptor. The top 20 energy-based compounds were visualized in a heatmap plot (Figure 3). The cells with high binding score were given a steel blue color, whereas those of a low binding score were shown in white color. Further, these data sets were normalized and illustrated in the scatter plot. Each point represented the mean normalized score of an individual compound toward three targets, called the z-score. We found that some compounds having the highest binding affinities were placed in different ranks and orders in each receptor. The selection of candidates for further analysis is a critical issue. To avoid selection biases, we used the unsupervised machine learning cluster methods to help in understanding and identifying the promising compounds as multitarget candidates.

3.3. Cluster Methods Help Identifying the Multitarget Candidates

All sets of the docking result were used to sequentially perform the cluster analysis using the ‘factoextra’ package in RStudio. First, the docked scores of all compounds were normalized to equal the weight of all groups. To define the optimal number of clusters (k), we used the average silhouette width and the most frequent indices as the parameters. The average silhouette method is used for evaluating the quality of clustering by computing the average observations (n) for different values of k. Generally, the optimal number of clusters is the one that maximizes the average silhouette width over a range of possible values of clusters. Our results showed that two clusters maximize the average silhouette coefficients, with three clusters coming in as the second optimal number of clusters (Figure 4a). This study was in agreement with the result of the frequency among indices (Figure 4b). The bar plot indicated that the most frequent optimal number of clusters for all indices is two. No difference was observed in the result of both validation methods. We reasoned that the optimal number of clusters was defined into two clusters that possessed 82 and 26 in Cluster 1 (blue) and 2 (yellow), respectively. Further, the clustering data were visualized in the dimension plot (Figure 4c). A clustering data set performed the principal component analysis (PCA) to demonstrate the dominant features of this model and were plotted; the data points followed the first two principal components. These two PCs explain 98.5% of the point variability, which obtained the high variations of data. In addition, the hierarchical clustering algorithm was also used to confirm the accuracy of our clustering result. A hierarchical dendrogram helps in illustrating the clustering data, for which it was suggested that our data could be defined into two clusters (Figure 4d). Interestingly, the compounds with high binding affinities toward three receptors were placed in Cluster 2. Thus, we selected Cluster 2 for further analysis.

Subsequently, all the observations in Cluster 2 were subjected to the subcluster analysis using the hierarchical k-means algorithm. This method uses a combination of the hierarchical and k-means algorithm to determine the cluster groups. The optimal numbers of clusters were defined as previously described. In the subclustering analysis, the proposed number of clusters was 3, according to the majority rules of all indices calculated from the nbclust package in R (Figure 5a). Afterward, the silhouette method was used to evaluate the quality of the clustering result (Figure 5b). This method calculates the distance between the clusters, called the silhouette width (Si). In general, a high value of Si suggests the well-clustered observation, whereas the negative value of the observation indicates that the observation may be placed in the wrong cluster. Then, the cluster groups of the observations were identified (Figure 5c).

Further, we performed a clustering heatmap plot. This method provides data visualization by color coding and correlation by using a dendrogram from the clustering method. The data from the subcluster analysis were prepared for a distance matrix prior to clustering into several groups using a hierarchical method. The subcluster data were clustered hierarchically by binding scores and by targets with color coding of dark green indicating high binding affinity and dark brown indicating low binding affinity (Figure 5d). It can be noted from the subcluster results that several compounds tended to cluster together. Particularly in the green color group, the compound including rutin, sagittatoside A, icariin, and kaempferitrin showed the potent binding affinities toward the individual target. According to the heatmap plot and subcluster analysis, only the top four compounds represented the promising candidates distinguished from our studies.

3.4. Binding Mode Analysis of Candidates to Cathepsin K

The binding interactions of docked complexes were identified from LigandScout according to the manufacture’s guideline. Briefly, the distance of hydrogen bond forming is in the range of 2.2 to 3.8 angstrom and that of the hydrophobic interaction is 1.0 to 5.9 angstrom.

The binding modes of the promising compounds toward the active site of Cathepsin K are depicted in Figure 6. The active site showed various amino acids involved in ligand binding. Cys25 and His162 represented conserved crucial residues for the active site of Cathepsin K.

Rutin was located along the V-shape active cleft formed by the R- and L-domains. It was observed that rutin can interact with Gly19, Gly23, Cys25, Gly66, Tyr67, and Leu160 through the hydrogen bonding formation. The hydroxyl groups of benzene A and C interact with the sidechain of Gly19, Gly23, and Cys25. Additionally, the glycoside moiety of rutin can bind to Leu160 using the hydrophobic interaction. Sagittatoside A forms the hydrogen bonding through the binding of Cys25, Gly66, Tyr67, Asp158, Asn159, and Leu160. A benzene ring of sagittatoside A core forms the hydrophobic interaction with Tyr67 and Ala163. Icariin also forms the hydrogen bond through Cys25, Gly66, Met68, and Asn161. The methyl group of glucoside connected with C3 can interact with Trp184 via hydrophobic interaction. Kaempferitrin can interact with Cys25, Gly66, and His162 via the hydrogen bonding formation. This compound also binds to residues Tyr67, Ala134, Leu160, Trp184, and Leu209 via hydrophobic interactions.

Accordingly, the important interactions were made of the strong hydrogen bond between the candidates and the residue Cys25 and Gly66, which represents the amino acid forming catalyzing site of Cathepsin K. Interestingly, the volume of the binding site of Cathepsin K is limited due to the narrow active cleft formed by the R- and L-domains, especially the bulky configuration of the phenyl ring of Tyr67. Regarding these binding modes, we observed that all candidates fitted well in the active pocket and preferred the binding at the R-domain more than the L-domain.

3.5. Binding Mode Analysis of Candidates to V-ATPase

Inhibitors of osteoclast V-ATPase are expected to become therapeutic agents for the treatment of several bone diseases that are caused by an increase in bone resorption. The assembly of AB complexes requires ATP binding to help the protein fold correctly [48]. The hydrolysis of ATP within AB complex provides the energy for the rotation of the V0 domain, related to the physiological function of V-ATPase. The nucleotide binding site of V-ATPase is located between subunits A and B. The ATP hydrolysis occurs depending on the configuration of Arg-finger (Arg350) [49]. The amino acids that involved ligand binding include Gly235, Ala236, Gly237, Lys238, Tyr239, Val240, Phe425, Gln503, and Ala505 from subunit A, and Arg350 from subunit B [36]. As depicted in Figure 7, kaempferitrin exhibited the highest binding affinity toward V-ATPase, followed by rutin, sagittatoside A, and icariin, respectively. Kaempferitrin can form six hydrogen bonds to amino acids within subunit A, including Gly235, Lys238, Thr239, Glu261, Arg262, and Arg333, and three hydrogen bonds to amino acids within subunit B, such as Tyr123, Tyr321, and Arg350. Rutin can form hydrogen interactions with Gly235 and Gln503 in subunit A, and Arg350 and Asp353 in subunit B. Sagittatoside A attaches with Asn504 and Ala505 in subunit A, and Tyr123 in subunit B, through hydrogen bond formation. Icariin can form the hydrogen bonds with Phe425 in subunit A, and Tyr123, Arg350, and Lys352 in subunit B. Regarding these hydrogen bond formations, Tyr123 and Arg350 represent the crucial amino acids in subunit B. Additionally, these candidates also interact with several residues involving the binding cavity, especially Phe506 from subunit A and Leu348 from subunit B.

3.6. Binding Mode Analysis of Candidates to α_Vβ₃ Integrin Receptor

Ligands of integrin α_vβ₃ are expected to have a major impact on the treatment of several human diseases, such as osteoporosis, rheumatoid arthritis, and cancer. The RGD binding site is located at the interface between the β-propeller domain and the β I-like domain. Amino acids involving ligand interactions include Asp150, Asp218, and Tyr178 from chain A, and Ser121, Ser123, Asn215, and Arg216 from chain B [50]. Our docked ligands were located in the same location of RGD-native ligand. As in the docking result, rutin also presented the high binding affinity against the integrin receptor, followed by sagittatoside A, icariin, and kaempferitrin, respectively. The binding poses of candidates are depicted in Figure 8. Rutin can interact with integrin through hydrogen bond formation with Ala149 from chain A, and Ser121, Ser123, Tyr166, Arg214, Asn215, Arg216, and Asp251 from chain B. Sagittatoside A forms two hydrogen bonds with Ser123, and Asn215 from chain B. Icariin forms a hydrogen bond via interaction with Tyr178 from chain A, and Tyr166, Arg214, Asn215, Arg216 from chain B. Meanwhile, kaempferitrin forms hydrogen bonds with residues from chain A, including Asp148, Asp150, Tyr178, and Asp219, and those from chain B including Tryr166, Arg214, Asn215, Arg216, and Asp217. Interestingly, the important hydrogen interactions can be found between all candidates and Asn215 from chain A, suggesting that Asn215 is a crucial amino acid for this study.

4. Discussion

Osteoporosis is a complex disorder triggered by several pathways. Even though the individual target studies are well documented, the investigation on phytochemicals toward multiple osteoclastic targets has not been yet elucidated. In this present study, we first studied the herbal substances that reported osteoprotective activities from previous studies [27,28] with multitarget activity on osteoclasts by performing molecular docking and cluster analysis.

Defining the potential target is the critical issue for drug discovery nowadays [51]. Therefore, we selected the promising multiple osteoclastic targets [52] in advance, including Cathepsin K, V-ATPase, and α_Vβ₃ integrin receptor to explore the multitarget compounds using the target-based virtual screening. Network analysis provides a useful tool to study the complexity of diseases. Here, we reported an approach to study the PPI network based on the gene expressed in targeted proteins of osteoclasts. The result we presented here indicates that proteins, such as Cathepsin K, α_Vβ₃ integrin receptor, and V-ATPase, share common characteristics and direct interactions. We then hypothesized that the inhibition of these proteins might promote the synergistic effect in the reduction in osteoclast function.

We further performed the virtual screening of herbal substances toward three osteoclastic targets. However, the limitation of molecular docking is due to the low accuracy of the scoring function [53]. Thus, we performed a predictive model for the selection of the candidates, utilizing a combination of molecular docking and clustering approaches. The use of the clustering algorithm in filtering the virtual screening in the preliminary study helps to reduce the number of interested compounds with an unbiased algorithm. Several studies demonstrated the advantages of utilizing the clustering method as an additional method that efficiently helps in identifying the active compounds [54,55]. Our study used both k-means and hierarchical clustering methods. These methods allow to find the reliable cluster and define the phytochemical candidates.

As highlighted in our results, we defined four compounds, namely, rutin, sagittatoside A, icariin, and kaempferitrin, obtained the multitarget activity against osteoclastic receptors. Our finding is in agreement with the previous experiments that showed that rutin exhibited an osteoprotective effect by improving the bone mineral density (BMD) and bone strength in the ovariectomized (OVX) rat [56,57]. Additionally, rutin reduced the levels of acid phosphatases secreted during bone resorption in osteoclasts and increased the osteocalcin activity, a biomarker enzyme for bone formation found in osteoblasts [57].

Sagittatoside A and icariin have been reported as the bioactive compounds from Herba epimedii extract. The herb extract has gained more interest in various studies due to its beneficial activities, such as anti-oxidative, anti-inflammatory, and anti-tumor properties [58,59]. This herb has been used for a decade to treat the cardiovascular disease, impotence, neurological disorder, and osteoporosis [60]. Given that Epimedii extract to the OVX rats exhibited significantly increased bone strength via potent inhibiting the alkaline and acid phosphatase suggesting that the extract of Epimedii helps to improve postmenopausal osteoporosis [61].

Icariin is a prenylated flavonol glycoside that presents as a main bioactive compound retrieved form Epimedium extract. Icariin possesses the osteoprotective property through the stimulating osteogenic differentiation of mesenchymal stem cells (MSCs), which related to the bone formation, and inhibiting the osteoclastogenic differentiation along with the bone resorption by osteoclasts [62,63,64]. Interestingly, icariin showed 64% inhibition toward Cathepsin K at the concentration of 1 mM [65]. Icariin dramatically increased the ATP level which was related to the inhibition of mitochondria apoptosis [66]. Insight into the structure, the presence of a prenyl moiety at the position 8 suggested the relevance to osteogenic activity [67], and enhanced the anti-osteoclastic activity [68].

Sagittatoside A (or icariin A) is a phytochemical compound from Epimedii extract. Even a small amount recovered from extract, but sagittatoside A possess the potential effects on several tissues [61]. Treatment with sagittatoside A reported a strong osteoprotective effect, as well as icariin, in terms of increasing the proliferation and differentiation of osteoblast and suppressing osteoclastogenesis.

Kaempferitrin or lespedin is a glycosyloxyflavone, constituted by two rhamnose molecules bound to kaempferol at positions 3 and 7 through glycosidic linkage. This compound can be isolated from several plants and fungus. In the OVX rat, the given kaempferitrin considerably improved the bone mass and microarchitecture and inhibited the osteoclastic function from in vitro assay [69]. Ye et al. investigated the administration of Podocarpium extract consisted of kaempferitrin, the most abundant phytochemical, demonstrated the decreasing of bone resorption markers, acid phosphatase and Cathepsin K, in osteoclastic cells [69,70].

Our predictive model indicated that rutin, sagittatoside A, icariin, and kaempferitrin were chosen on the basis of their high energy-based docking scores against osteoclasts on specific Cathepsin K, V-ATPase, and α_Vβ₃ integrin receptors. We assume that these compounds might show a positive effect on the bone mechanism, particularly anti-osteoclastic activity. In particular, rutin demonstrated the highest energetically favorable interactions against three osteoclastic targets, according to molecular docking and clustering studies. Further investigations are needed to clarify the effects of each promising candidate on individual targets, especially molecular dynamics (MD), and the evaluation of anti-osteoclastic activity from in vitro and in vivo assays should be performed.

5. Conclusions

In this study, we carried out a virtual screening of 108 herbal compounds and performed the cluster analyses to identify the multitarget herbal compounds against Cathepsin K, V-ATPase, and α_Vβ₃ integrin receptors. A combination of several in silico methods helped to reduce the number of interesting compounds with the unbiased algorithm. Thus, we found that rutin, sagittatoside A, icariin, and kaempferitrin showed to be preferably binding to selected proteins involving osteoclast function, suggesting the idealistic candidates with anti-osteoclastic activity. These candidates, particularly rutin, should be considered for further evaluation using in vitro and in vivo studies. To sum up, our present work provides a reliable and straightforward method to identify multitarget candidates for further experimental analyses.