Next Article in Journal
Reliability Analysis of Military Vehicles Based on Censored Failures Data
Previous Article in Journal
Impact of Radial Lands on the Reduction of Torque/Force Generation of a Heat-Treated Nickel-Titanium Rotary Instrument
 
 
Search for Articles:
Title / Keyword
Author / Affiliation / Email
Journal
Article Type
 
 
Advanced Search
 
Section
Special Issue
Volume
Issue
Number
Page
 
Journals
Applied Sciences
Volume 12
Issue 5
10.3390/app12052621
applsci-logo
Submit to this Journal Review for this Journal Propose a Special Issue
► Article Menu

Article Menu

share Share announcement Help format_quote Cite question_answer Discuss in SciProfiles thumb_up
...
Endorse textsms
...
Comment
first_page
settings
Order Article Reprints
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

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

by
Siripat Chaichit
1,2,*,
Pathomwat Wongrattanakamon
2,
Busaban Sirithunyalug
2,
Piyarat Nimmanpipug
3 and
Supat Jiranusornkul
2,*
1
PhD Degree Program in Pharmacy, Faculty of Pharmacy, Chiang Mai University, Chiang Mai 50200, Thailand
2
Laboratory of Molecular Design and Modeling (LMDS), Department of Pharmaceutical Sciences, Faculty of Pharmacy, Chiang Mai University, Chiang Mai 50200, Thailand
3
Faculty of Science, Chiang Mai University, Chiang Mai 50200, Thailand
*
Authors to whom correspondence should be addressed.
Appl. Sci. 2022, 12(5), 2621; https://doi.org/10.3390/app12052621
Submission received: 21 January 2022 / Revised: 24 February 2022 / Accepted: 25 February 2022 / Published: 3 March 2022
Browse Figures
Versions Notes

Abstract

:
Osteoporosis is a complex bone disease indicating porous bone with low bone mass density and fragility. Cathepsin K, V-ATPase, and αVβ3 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β3 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β3 type of integrin receptor is highly expressed in osteoclasts. The inhibition of the αVβ3 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β3 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β3 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β3 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 Å3. 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 Å3 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β3 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 Å3 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−16, 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−11. 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β3 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β3 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β3 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β3 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β3 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β3 Integrin Receptor

Ligands of integrin αvβ3 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β3 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β3 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β3 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β3 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.

Author Contributions

Conceptualization, S.C. and S.J.; methodology, S.C., S.J. and P.W.; data curation, P.W., B.S., P.N. and S.J.; writing, S.C.; visualization, S.C.; supervision, S.J. All authors have read and agreed to the published version of the manuscript.

Funding

This work is supported by the National Research Council of Thailand (fiscal year 2018) and partially supported by TA/RA Scholarship from the Graduate School, Chiang Mai University.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Not applicable.

Acknowledgments

The author would like to thank Yang Peng for providing a guidance regarding analysis using R.

Conflicts of Interest

The authors declare no conflict of interest.

Abbreviation

V-ATPase: Vacuole-H+ ATPase.

References

  1. Tanaka, Y.; Nakayamada, S.; Okada, Y. Osteoblasts and osteoclasts in bone remodeling and inflammation. Curr. Drug Targets Inflamm. Allergy 2005, 4, 325–328. [Google Scholar] [CrossRef] [PubMed]
  2. Boyce, B.F. Advances in the regulation of osteoclasts and osteoclast functions. J. Dent. Res. 2013, 92, 860–867. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  3. Novack, D.V.; Mbalaviele, G. Osteoclasts-Key Players in Skeletal Health and Disease. Microbiol. Spectr. 2016, 4, 3. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  4. Hartman, G.D.; Duggan, M.E. αvβ3 Integrin antagonists as inhibitors of bone resorption. Expert Opin. Investig. Drugs 2000, 9, 1281–1291. [Google Scholar] [CrossRef] [PubMed]
  5. Horton, M.A.; Helfrich, M.H. Integrins and development: Integrins in skeletal cell function and development. In Madame Curie Biosci. Database [Internet]; Landes Bioscience: Austin, TX, USA, 2013. [Google Scholar]
  6. Forgac, M. Vacuolar ATPases: Rotary proton pumps in physiology and pathophysiology. Nat. Rev. Mol. Cell Biol. 2007, 8, 917–929. [Google Scholar] [CrossRef]
  7. Huss, M.; Wieczorek, H. Inhibitors of V-ATPases: Old and new players. J. Exp. Biol. 2009, 212, 341–346. [Google Scholar] [CrossRef] [Green Version]
  8. Zhao, H. Membrane Trafficking in Osteoblasts and Osteoclasts: New Avenues for Understanding and Treating Skeletal Diseases. Traffic 2012, 13, 1307–1314. [Google Scholar] [CrossRef] [Green Version]
  9. Lacombe, J.; Karsenty, G.; Ferron, M. Regulation of lysosome biogenesis and functions in osteoclasts. Cell Cycle 2013, 12, 2744–2752. [Google Scholar] [CrossRef] [Green Version]
  10. Hou, W.S.; Li, Z.; Gordon, R.E.; Chan, K.; Klein, M.J.; Levy, R.; Keysser, M.; Keyszer, G.; Brömme, D. Cathepsin k is a critical protease in synovial fibroblast-mediated collagen degradation. Am. J. Pathol. 2001, 159, 2167–2177. [Google Scholar] [CrossRef] [Green Version]
  11. Brömme, D.; Lecaille, F. Cathepsin K inhibitors for osteoporosis and potential off-target effects. Expert Opin. Investig. Drugs 2009, 18, 585–600. [Google Scholar] [CrossRef] [Green Version]
  12. Stoch, S.A.; Wagner, J.A. Cathepsin K Inhibitors: A Novel Target for Osteoporosis Therapy. Clin. Pharmacol. Ther. 2008, 83, 172–176. [Google Scholar] [CrossRef] [PubMed]
  13. Yuan, H.; Ma, Q.; Cui, H.; Liu, G.; Zhao, X.; Li, W.; Piao, G. How Can Synergism of Traditional Medicines Benefit from Network Pharmacology? Molecules 2017, 22, 1135. [Google Scholar] [CrossRef] [PubMed]
  14. Caesar, L.K.; Cech, N.B. Synergy and antagonism in natural product extracts: When 1 + 1 does not equal 2. Nat. Prod. Rep. 2019, 36, 869–888. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  15. Zhou, X.; Seto, S.W.; Chang, D.; Kiat, H.; Razmovski-Naumovski, V.; Chan, K.; Bensoussan, A. Synergistic Effects of Chinese Herbal Medicine: A Comprehensive Review of Methodology and Current Research. Front. Pharmacol. 2016, 7, 201. [Google Scholar] [CrossRef] [Green Version]
  16. D’Archivio, M.; Filesi, C.; Di Benedetto, R.; Gargiulo, R.; Giovannini, C.; Masella, R. Polyphenols, dietary sources and bioavailability. Ann. Ist. Super. Sanita 2007, 43, 348–361. [Google Scholar]
  17. Scalbert, A.; Johnson, I.T.; Saltmarsh, M. Polyphenols: Antioxidants and beyond. Am. J. Clin. Nutr. 2005, 81, 215S–217S. [Google Scholar] [CrossRef]
  18. Scalbert, A.; Manach, C.; Morand, C.; Rémésy, C.; Jiménez, L. Dietary Polyphenols and the Prevention of Diseases. Crit. Rev. Food Sci. Nutr. 2005, 45, 287–306. [Google Scholar] [CrossRef]
  19. Weaver, C.M.; Alekel, D.L.; Ward, W.E.; Ronis, M.J. Flavonoid Intake and Bone Health. J. Nutr. Gerontol. Geriatr. 2012, 31, 239–253. [Google Scholar] [CrossRef]
  20. Medina-Franco, J.L.; Giulianotti, M.A.; Welmaker, G.S.; Houghten, R.A. Shifting from the single to the multitarget paradigm in drug discovery. Drug Discov. Today 2013, 18, 495–501. [Google Scholar] [CrossRef] [Green Version]
  21. Ramsay, R.R.; Popovic-Nikolic, M.R.; Nikolic, K.; Uliassi, E.; Bolognesi, M.L. A perspective on multi-target drug discovery and design for complex diseases. Clin. Transl. Med. 2018, 7, 3. [Google Scholar] [CrossRef] [Green Version]
  22. Ma, X.H.; Shi, Z.; Tan, C.; Jiang, Y.; Go, M.L.; Low, B.C.; Chen, Y.Z. In-silico approaches to multi-target drug discovery: Computer aided multi-target drug design, multi-target virtual screening. Pharm. Res. 2010, 27, 739–749. [Google Scholar] [CrossRef] [PubMed]
  23. Sliwoski, G.; Kothiwale, S.; Meiler, J.; Lowe, E.W., Jr. Computational methods in drug discovery. Pharmacol. Rev. 2013, 66, 334–395. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  24. Hartigan, J.A.; Wong, M.A. A K-means clustering algorithm. J. R. Stat. Soc. Ser. C (Appl. Stat.) 1979, 28, 100–108. [Google Scholar]
  25. MacQueen, J. Some methods for classification and analysis of multivariate observations. In Proceedings of the 5th Berkeley Symposium on Mathematical Statistics and Probability; Cambridge University Press: Oakland, CA, USA, 1967; Volume 1, pp. 281–297. [Google Scholar]
  26. Ward Jr, J.H. Hierarchical grouping to optimize an objective function. J. Am. Stat. Assoc. 1963, 58, 236–244. [Google Scholar] [CrossRef]
  27. Štulíková, K.; Karabín, M.; Nešpor, J.; Dostálek, P. Therapeutic Perspectives of 8-Prenylnaringenin, a Potent Phytoestrogen from Hops. Molecules 2018, 23, 660. [Google Scholar] [CrossRef] [Green Version]
  28. Che, C.-T.; Wong, M.S.; Lam, C.W.K. Natural Products from Chinese Medicines with Potential Benefits to Bone Health. Molecules 2016, 21, 239. [Google Scholar] [CrossRef] [Green Version]
  29. Szklarczyk, D.; Franceschini, A.; Wyder, S.; Forslund, K.; Heller, D.; Huerta-Cepas, J.; Simonovic, M.; Roth, A.; Santos, A.; Tsafou, K.P.; et al. STRING v10: Protein-protein interaction networks, integrated over the tree of life. Nucleic Acids Res. 2015, 43, D447–D452. [Google Scholar] [CrossRef]
  30. Shannon, P.; Markiel, A.; Ozier, O.; Baliga, N.S.; Wang, J.T.; Ramage, D.; Amin, N.; Schwikowski, B.; Ideker, T. Cytoscape: A software environment for integrated models of biomolecular interaction networks. Genome Res. 2003, 13, 2498–2504. [Google Scholar] [CrossRef]
  31. Frisch, A. Gaussian 09W Reference; Gaussian: Wallingford, CT, USA, 2009; 25p. [Google Scholar]
  32. Hehre, W.; Ditchfield, R.; Stewart, R.F.; Pople, J.A. self-consistent molecular orbital methods. iv. use of Gaussian expansions of Slater-type orbitals. Extension to second-row molecules. J. Chem. Phys. 1970, 52, 2769–2773. [Google Scholar] [CrossRef]
  33. Morris, G.M.; Huey, R.; Lindstrom, W.; Sanner, M.F.; Belew, R.K.; Goodsell, D.S.; Olson, A.J. AutoDock4 and AutoDockTools4: Automated Docking with Selective Receptor Flexibility. J. Comput. Chem. 2009, 30, 2785–2791. [Google Scholar] [CrossRef] [Green Version]
  34. Borišek, J.; Vizovišek, M.; Sosnowski, P.; Turk, B.; Turk, D.; Mohar, B.; Novič, M. Development of N-(Functionalized benzoyl)-homocycloleucyl-glycinonitriles as Potent Cathepsin K Inhibitors. J. Med. Chem. 2015, 58, 6928–6937. [Google Scholar] [CrossRef] [PubMed]
  35. BIOVIA. Discovery Studio Visualizer; Dassault Systèmes: San Diego, CA, USA, 2017; 936p. [Google Scholar]
  36. Suzuki, K.; Mizutani, K.; Maruyama, S.; Shimono, K.; Imai, F.L.; Muneyuki, E.; Kakinuma, Y.; Ishizuka-Katsura, Y.; Shirouzu, M.; Yokoyama, S.; et al. Crystal structures of the ATP-binding and ADP-release dwells of the V1 rotary motor. Nat. Commun. 2016, 7, 13235. [Google Scholar] [CrossRef] [PubMed]
  37. Dong, X.; Mi, L.-Z.; Zhu, J.; Wang, W.; Hu, P.; Luo, B.-H.; Springer, T.A. αVβ3 Integrin Crystal Structures and Their Functional Implications. Biochemistry 2012, 51, 8814–8828. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  38. R: A Language and Environment for Statistical Computing; R Foundation for Statistical Computing: Vienna, Austria, 2020.
  39. Charrad, M.; Ghazzali, N.; Boiteau, V. NbClust: An r package for determining the relevant number of clusters in a data set. J. Stat. Softw. 2014, 61, 1. [Google Scholar] [CrossRef] [Green Version]
  40. Rousseeuw, P.J. Silhouettes: A graphical aid to the interpretation and validation of cluster analysis. J. Comput. Appl. Math. 1987, 20, 53–65. [Google Scholar] [CrossRef] [Green Version]
  41. Schrodinger, L.; DeLano, W. The PyMOL Molecular Graphics System Version 1.8. Am. J. Infect. Dis. Microbiol. 2015, 4, 61–71. [Google Scholar]
  42. Kolde, R.; Kolde, M.R. Package ‘pheatmap’. R Packag. 2015, 1, 790. [Google Scholar]
  43. Gu, Z.; Eils, R.; Schlesner, M. Complex heatmaps reveal patterns and correlations in multidimensional genomic data. Bioinformatics 2016, 32, 2847–2849. [Google Scholar] [CrossRef] [Green Version]
  44. Kassambara, A.; Mundt, F. Factoextra: Extract and visualize the results of multivariate data analyses. R Packag. 2017, 1, 337–354. [Google Scholar]
  45. Wolber, G.; Langer, T. LigandScout:  3-D Pharmacophores Derived from Protein-Bound Ligands and Their Use as Virtual Screening Filters. J. Chem. Inf. Model. 2005, 45, 160–169. [Google Scholar] [CrossRef]
  46. Sørensen, M.G.; Henriksen, K.; Neutzsky-Wulff, A.V.; Dziegiel, M.H.; Karsdal, M.A. Diphyllin, a novel and naturally potent V-ATPase inhibitor, abrogates acidification of the osteoclastic resorption lacunae and bone resorption. J. Bone Miner. Res. 2007, 22, 1640–1648. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  47. Buommino, E.; Baroni, A.; Canozo, N.; Petrazzuolo, M.; Nicoletti, R.; Vozza, A.; Tufano, M.A. Artemisinin reduces human melanoma cell migration by down-regulating αVβ3 integrin and reducing metalloproteinase 2 production. Invest. New Drugs 2009, 27, 412–418. [Google Scholar] [CrossRef] [PubMed]
  48. Imamura, H.; Funamoto, S.; Yoshida, M.; Yokoyama, K. Reconstitution in vitro of V1 complex of Thermus thermophilus V-ATPase revealed that ATP binding to the A subunit is crucial for V1 formation. J. Biol. Chem. 2006, 281, 38582–38591. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  49. Arai, S.; Saijo, S.; Suzuki, K.; Mizutani, K.; Kakinuma, Y.; Ishizuka-Katsura, Y.; Ohsawa, N.; Terada, T.; Shirouzu, M.; Yokoyama, S.; et al. Rotation mechanism of Enterococcus hirae V1-ATPase based on asymmetric crystal structures. Nature 2013, 493, 703–707. [Google Scholar] [CrossRef]
  50. Craig, D.; Gao, M.; Schulten, K.; Vogel, V. Structural Insights into How the MIDAS Ion Stabilizes Integrin Binding to an RGD Peptide under Force. Structure 2004, 12, 2049–2058. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  51. Talevi, A. Multi-target pharmacology: Possibilities and limitations of the “skeleton key approach” from a medicinal chemist perspective. Front. Pharmacol. 2015, 6, 205. [Google Scholar] [CrossRef] [Green Version]
  52. Roux, S. New treatment targets in osteoporosis. Jt. Bone Spine 2010, 77, 222–228. [Google Scholar] [CrossRef]
  53. Gimeno, A.; Ojeda-Montes, M.J.; Tomás-Hernández, S.; Cereto-Massagué, A.; Beltrán-Debón, R.; Mulero, M.; Pujadas, G.; Garcia-Vallvé, S. The Light and Dark Sides of Virtual Screening: What Is There to Know? Int. J. Mol. Sci. 2019, 20, 1375. [Google Scholar] [CrossRef] [Green Version]
  54. Belkadi, A.; Kenouche, S.; Melkemi, N.; Daoud, I.; Djebaili, R. K-means clustering analysis, ADME/pharmacokinetic prediction, MEP, and molecular docking studies of potential cytotoxic agents. Struct. Chem. 2021, 32, 2235–2249. [Google Scholar] [CrossRef]
  55. Bouvier, G.; Evrard-Todeschi, N.; Girault, J.-P.; Bertho, G. Automatic clustering of docking poses in virtual screening process using self-organizing map. Bioinformatics 2010, 26, 53–60. [Google Scholar] [CrossRef] [Green Version]
  56. Wang, Q.-L.; Huo, X.-C.; Wang, J.-H.; Wang, D.-P.; Zhu, Q.-L.; Liu, B.; Xu, L.-L. Rutin prevents the ovariectomy-induced osteoporosis in rats. Eur. Rev. Med. Pharmacol. Sci. 2017, 21, 1911–1917. [Google Scholar]
  57. Gera, S.; Pooladanda, V.; Godugu, C.; Swamy Challa, V.; Wankar, J.; Dodoala, S.; Sampathi, S. Rutin nanosuspension for potential management of osteoporosis: Effect of particle size reduction on oral bioavailability, in vitro and in vivo activity. Pharm. Dev. Technol. 2020, 25, 971–988. [Google Scholar] [CrossRef] [PubMed]
  58. Zhang, H.-F.; Zhang, X.; Yang, X.-H.; Qiu, N.-X.; Wang, Y.; Wang, Z.-Z. Microwave assisted extraction of flavonoids from cultivated Epimedium sagittatum: Extraction yield and mechanism, antioxidant activity and chemical composition. Ind. Crops Prod. 2013, 50, 857–865. [Google Scholar] [CrossRef]
  59. Sze, S.C.W.; Tong, Y.; Ng, T.B.; Cheng, C.L.Y.; Cheung, H.P. Herba Epimedii: Anti-oxidative properties and its medical implications. Molecules 2010, 15, 7861–7870. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  60. Wang, L.; Li, Y.; Guo, Y.; Ma, R.; Fu, M.; Niu, J.; Gao, S.; Zhang, D. Herba Epimedii: An Ancient Chinese Herbal Medicine in the Prevention and Treatment of Osteoporosis. Curr. Pharm. Des. 2016, 22, 328–349. [Google Scholar] [CrossRef]
  61. Zhao, B.; Wang, J.; Song, J.; Wang, C.; Gu, J.; Yuan, J.; Zhang, L.; Jiang, J.; Feng, L.; Jia, X. Beneficial Effects of a Flavonoid Fraction of Herba Epimedii on Bone Metabolism in Ovariectomized Rats. Planta Med. 2016, 82, 322–329. [Google Scholar] [CrossRef] [Green Version]
  62. Wang, Z.; Wang, D.; Yang, D.; Zhen, W.; Zhang, J.; Peng, S. The effect of icariin on bone metabolism and its potential clinical application. Osteoporos. Int. 2018, 29, 535–544. [Google Scholar] [CrossRef]
  63. Wang, J.; Tao, Y.; Ping, Z.; Zhang, W.; Hu, X.; Wang, Y.; Wang, L.; Shi, J.; Wu, X.; Yang, H.; et al. Icariin attenuates titanium-particle inhibition of bone formation by activating the Wnt/β-catenin signaling pathway in vivo and in vitro. Sci. Rep. 2016, 6, 23827. [Google Scholar] [CrossRef]
  64. Kim, B.; Lee, K.Y.; Park, B. Icariin abrogates osteoclast formation through the regulation of the RANKL-mediated TRAF6/NF-κB/ERK signaling pathway in Raw264.7 cells. Phytomedicine 2018, 51, 181–190. [Google Scholar] [CrossRef]
  65. Sun, P.; Liu, Y.; Deng, X.; Yu, C.; Dai, N.; Yuan, X.; Chen, L.; Yu, S.; Si, W.; Wang, X.; et al. An inhibitor of Cathepsin K, icariin suppresses cartilage and bone degradation in mice of collagen-induced arthritis. Phytomedicine 2013, 20, 975–979. [Google Scholar] [CrossRef]
  66. Li, H.; Zhang, X.; Zhu, X.; Qi, X.; Lin, K.; Cheng, L. The Effects of Icariin on Enhancing Motor Recovery Through Attenuating Pro-Inflammatory Factors and Oxidative Stress via Mitochondrial Apoptotic Pathway in the Mice Model of Spinal Cord Injury. Front. Physiol. 2018, 9, 1617. [Google Scholar] [CrossRef] [PubMed]
  67. Ming, L.-G.; Chen, K.-M.; Xian, C.J. Functions and action mechanisms of flavonoids genistein and icariin in regulating bone remodeling. J. Cell. Physiol. 2013, 228, 513–521. [Google Scholar] [CrossRef] [PubMed]
  68. Lv, X.; Chen, K.-M.; Ge, B.-F.; Ma, H.-P.; Song, P.; Cheng, K. Comparative study on effect of 8-prenylnaringenin and narigenin on activity of osteoclasts cultured in vitro. Zhongguo Zhong Yao Za Zhi 2013, 38, 1992–1996. [Google Scholar] [PubMed]
  69. Ma, X.-Q.; Han, T.; Zhang, X.; Wu, J.-Z.; Rahman, K.; Qin, L.-P.; Zheng, C.-J. Kaempferitrin prevents bone lost in ovariectomized rats. Phytomedicine 2015, 22, 1159–1162. [Google Scholar] [CrossRef]
  70. Ye, Q.; Ma, X.-Q.; Hu, C.-L.; Lin, B.; Xu, L.-S.; Zheng, C.-J.; Qin, L.-P. Antiosteoporotic activity and constituents of Podocarpium podocarpum. Phytomedicine 2015, 22, 94–102. [Google Scholar] [CrossRef]
Applsci 12 02621 g001 550
Figure 1. Schematic representation of this study. The potential targets in different osteoclast stages were identified, including Cathepsin K, V-ATPase, and αVβ3 integrin receptor. These targets were investigated to find the interaction with phytochemicals from the previous literature through molecular docking. The docking results were clustered to find the potential cluster based on the AutoDock binding interactions. The compounds in the selected cluster were subclustered using the clustering method incorporated with the heatmap analysis to find the promising compounds for further studies.
Figure 1. Schematic representation of this study. The potential targets in different osteoclast stages were identified, including Cathepsin K, V-ATPase, and αVβ3 integrin receptor. These targets were investigated to find the interaction with phytochemicals from the previous literature through molecular docking. The docking results were clustered to find the potential cluster based on the AutoDock binding interactions. The compounds in the selected cluster were subclustered using the clustering method incorporated with the heatmap analysis to find the promising compounds for further studies.
Applsci 12 02621 g001
Applsci 12 02621 g002 550
Figure 2. The input genes present in the different box styles such as ellipse for Cathepsin K, triangle for integrin receptor, and round rectangle for V-ATPase. Nodes represent protein, and edges represent interactions. (a) The gene expression in bone tissues are mapped to the node color; nodes with high expression are colored dark blue, while those of low expression are colored gray. The blue gradient is used to define the location of genes in bone tissue. A total of 25 enriched modes with PPI enrichment p-value of <1.0 × 10−16. (b) The most significant genes from the PPI network (node color: yellow indicates the Cathepsin K gene, green indicates the V-ATPase genes, and red indicates the αVβ3 integrin genes).
Figure 2. The input genes present in the different box styles such as ellipse for Cathepsin K, triangle for integrin receptor, and round rectangle for V-ATPase. Nodes represent protein, and edges represent interactions. (a) The gene expression in bone tissues are mapped to the node color; nodes with high expression are colored dark blue, while those of low expression are colored gray. The blue gradient is used to define the location of genes in bone tissue. A total of 25 enriched modes with PPI enrichment p-value of <1.0 × 10−16. (b) The most significant genes from the PPI network (node color: yellow indicates the Cathepsin K gene, green indicates the V-ATPase genes, and red indicates the αVβ3 integrin genes).
Applsci 12 02621 g002
Applsci 12 02621 g003 550
Figure 3. The binding energies of phytochemicals. A heatmap represents top 20 molecular docking score of herbal compounds against Cathepsin K, V-ATPae, and αvβ3 integrin receptor. Color values correspond to binding energy assessed with AutoDock scoring function. The normalized value (z-score) of each docked compound represents an individual dot plot in the scatter plot.
Figure 3. The binding energies of phytochemicals. A heatmap represents top 20 molecular docking score of herbal compounds against Cathepsin K, V-ATPae, and αvβ3 integrin receptor. Color values correspond to binding energy assessed with AutoDock scoring function. The normalized value (z-score) of each docked compound represents an individual dot plot in the scatter plot.
Applsci 12 02621 g003
Applsci 12 02621 g004 550
Figure 4. Clustering analysis of docking results. The numbers of clusters are identified using (a) average silhouette width and (b) frequency among all indices. The clustering shows two groups, indicating Cluster 1 (blue) and Cluster 2 (yellow). The groups of clusters are visualized in the (c) cluster plot generated by the K-means clustering method, and (d) a dendrogram generated by the hierarchical clustering method.
Figure 4. Clustering analysis of docking results. The numbers of clusters are identified using (a) average silhouette width and (b) frequency among all indices. The clustering shows two groups, indicating Cluster 1 (blue) and Cluster 2 (yellow). The groups of clusters are visualized in the (c) cluster plot generated by the K-means clustering method, and (d) a dendrogram generated by the hierarchical clustering method.
Applsci 12 02621 g004
Applsci 12 02621 g005 550
Figure 5. Subclustering analysis. The compounds in the selected cluster were subclustered using k-means and the hierarchical approach. (a) Proposed number of clusters, (b) cluster silhouette calculation for internal clustering validation, (c) cluster dendrogram, and (d) a heatmap clustering plot defined by color indicates the top 4 compounds.
Figure 5. Subclustering analysis. The compounds in the selected cluster were subclustered using k-means and the hierarchical approach. (a) Proposed number of clusters, (b) cluster silhouette calculation for internal clustering validation, (c) cluster dendrogram, and (d) a heatmap clustering plot defined by color indicates the top 4 compounds.
Applsci 12 02621 g005
Applsci 12 02621 g006 550
Figure 6. The 2D and 3D binding poses of candidates in a binding site of Cathepsin K. Ligands (green) are bound with amino acids (white) at the binding site. The pharmacophore features are represented as green arrows (HBD), red arrows (HBA), and yellow spheres (hydrophobic property).
Figure 6. The 2D and 3D binding poses of candidates in a binding site of Cathepsin K. Ligands (green) are bound with amino acids (white) at the binding site. The pharmacophore features are represented as green arrows (HBD), red arrows (HBA), and yellow spheres (hydrophobic property).
Applsci 12 02621 g006
Applsci 12 02621 g007 550
Figure 7. The 2D and 3D binding poses of candidates in a binding site of V-ATPase. Ligands (green) are bound with amino acids in chain A (white) and chain B (pink). The pharmacophore features are represented as green arrows (HBD), red arrows (HBA), and yellow spheres (hydrophobic property).
Figure 7. The 2D and 3D binding poses of candidates in a binding site of V-ATPase. Ligands (green) are bound with amino acids in chain A (white) and chain B (pink). The pharmacophore features are represented as green arrows (HBD), red arrows (HBA), and yellow spheres (hydrophobic property).
Applsci 12 02621 g007
Applsci 12 02621 g008 550
Figure 8. The 2D and 3D binding poses of candidates in a binding site of integrin. Ligands (green) are bound with amino acids in chain A (white) and chain B (pink). The pharmacophore features are represented as green arrows (HBD), red arrows (HBA), and yellow spheres (hydrophobic property).
Figure 8. The 2D and 3D binding poses of candidates in a binding site of integrin. Ligands (green) are bound with amino acids in chain A (white) and chain B (pink). The pharmacophore features are represented as green arrows (HBD), red arrows (HBA), and yellow spheres (hydrophobic property).
Applsci 12 02621 g008
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Chaichit, S.; Wongrattanakamon, P.; Sirithunyalug, B.; Nimmanpipug, P.; Jiranusornkul, S. Multitarget-Based Virtual Screening for Identification of Herbal Substances toward Potential Osteoclastic Targets. Appl. Sci. 2022, 12, 2621. https://doi.org/10.3390/app12052621

AMA Style

Chaichit S, Wongrattanakamon P, Sirithunyalug B, Nimmanpipug P, Jiranusornkul S. Multitarget-Based Virtual Screening for Identification of Herbal Substances toward Potential Osteoclastic Targets. Applied Sciences. 2022; 12(5):2621. https://doi.org/10.3390/app12052621

Chicago/Turabian Style

Chaichit, Siripat, Pathomwat Wongrattanakamon, Busaban Sirithunyalug, Piyarat Nimmanpipug, and Supat Jiranusornkul. 2022. "Multitarget-Based Virtual Screening for Identification of Herbal Substances toward Potential Osteoclastic Targets" Applied Sciences 12, no. 5: 2621. https://doi.org/10.3390/app12052621

APA Style

Chaichit, S., Wongrattanakamon, P., Sirithunyalug, B., Nimmanpipug, P., & Jiranusornkul, S. (2022). Multitarget-Based Virtual Screening for Identification of Herbal Substances toward Potential Osteoclastic Targets. Applied Sciences, 12(5), 2621. https://doi.org/10.3390/app12052621

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

Article Metrics

Back to TopTop