Screening of Phytochemicals against Osteoporosis: Molecular Docking and Simulation-Based Computational Approaches

Bandar Hamad Aloufi¹*

¹Department of Biology, College of Science, University of Ha’il, Ha’il 2440, Saudi Arabia.

*Email: [email protected]; [email protected]

ABSTRACT

Osteoporosis is a skeletal disorder characterized by low bone mineral density and microarchitecture deterioration, which may result in fragility and fracture risk. Osteoporosis mostly affects postmenopausal women but elderly men may also be affected. It is a silent disease and reveals at the time of fracture. It is the most prevalent bone disease in the world. Currently, there is no cure available for osteoporosis. Through a literature survey, an association is found between the overexpressed genes and osteoporosis. E2 ubiquitin ligase TRAF3IP2 is reported to be overexpressed in osteoporotic vertebral fractures. This study was designed to develop the drug targets against osteoporosis with more therapeutic efficacy and least side effects. Proteins of the overexpressed genes were searched from the PDB database. After inhibitory protein was searched from literature and used as reference protein. A library of 13000 phytochemicals was docked against novel drug targets through Molecular Operating Environment (MOE). Absorption, distribution, metabolism, excretion, and toxicological analysis were done through ADMETsar. After careful analysis top four compounds Adenosine, 2'-deoxy-, 2-Amino-1-hydroxyoctadecan-3-one, (3,5-Dihydroxyoxolan-2-yl) methyl dihydrogen phosphate and 2-(3,7,11,15,19,23,27,31-Octamethyldotriaconta-2,6,10,14,18,22,26,30 -octaenyl) naphthalene-1,4-dione were reported. The top leading compound showed the least docking score of -14.32 kcal/mol. Then simulation was done for the top two compounds for better analysis and understanding of the process. This study provides novel insights into osteoporosis. The outcomes of this study are helpful to identify better drug candidates for the treatment of osteoporosis that will inhibit the function of the target protein thus suppressing the disease at the sub-clinical stage.

Key words:Osteoporosis, Vertebral fractures, Phytochemicals, Drug targets, Molecular docking

INTRODUCTION

Osteoporosis means porous bones. It is a systemic pathology of the skeleton due to which bones become fragile and the risk of fractures enhance. It is a major health problem in aging societies. It affects more than 200 million people worldwide and almost 8.9 million fractures are caused by osteoporotic fracture [1]. According to World Health Organization (WHO), osteoporosis is a depletion of bone mineral density (BMD) [2]. Basically, it is a silent disease and reveals only when the fractures appear. In fact, even after bones breakdown, about 80% of affected individuals still remain undiagnosed and not treated for osteoporosis, the responsible disease which has resulted in the fracture. Most people ignore the back pain and treated them as just a symbol of getting older and do not undergo the diagnosis and proper treatment. This disease often remains undiagnosed until the fractures appear so the only focus given to its medication therapies was to reduce the incidences of further fractures [3]. This disease is most common in menopause women but elderly men may also be affected by it. With the increasing age and menopause, the balance between the rate of bone restoration and formation is disturbed and restoration becomes higher than absorption so the risk of fractures becomes high [4].

It is characterized by two main types: primary and secondary. Primary osteoporosis is the common type of osteoporosis that is divided into postmenopausal osteoporosis (type 1) and age-associated osteoporosis (type 2) [5]. The causes of primary osteoporosis involve the loss of estrogen and androgen resulting in increased bone turnover and systemic senescence. Secondary osteoporosis may be due to a large and diverse group of medical disorders like diabetes, hyperparathyroidism, scurvy, liver disease, and malabsorption. The signs and symptoms of osteoporosis include bone pain or tenderness, fractures with little or no trauma, loss of height, depression, neck or lower back pain due to fractures, and stooped posture. These symptoms increase with the passage of time and age. Older women with bending back are commonly seen in our society [6].

Bones are essential in providing strength, structure, and protection to the body and organs. Normally in the period from childhood to adulthood approximately 30 years of age, the peak bone mass is reached. Then bone mass gradually starts to decline with the passage of age, and in women, this decline process is fast after menopause. Although bone mass is severely affected by the food qualities, diseases, and several adverse medications [7].

The human skeleton reinstates itself completely after every ten years. Osteoblasts, osteoclasts, and osteocytes are three major cells responsible for maintaining bone strength [8]. Osteoblasts and osteoclasts, coupled together, repair and renew the bones while osteocytes produce the regulatory signals by sensing. For a rigid skeleton, the balance between these units is very important. In osteoporosis, signaling in these three units is imbalanced or disturbed. Over time, net bone loss takes place because the osteoplastic activity increases with increasing age at the rate of about 1% every year after the age of 30 years [9]. The bones of the spine and hips are mostly affected by osteoporosis. The vertebrae of osteoporotic patients often become weak and shrink with time. As a result, height decreases and leading to a rounded back. Osteoporosis is diagnosed by measuring the BND of the spine and hip through dual-energy X-ray absorptiometry [10]. Large-scale genome-wide association studies (GWAS) provide a better understanding of the genetic makeup and mechanisms of this disease. The betterment of GWAS and post-GWAS techniques will make fruitful participation of genetics in clinical practice that will give better disease diagnosis and prevention [11].

TRAF3IP2 is the necessary adapter molecule in the IL-17 mediated signaling. IL-17 mediated signaling is important for the development of autoimmunity and provides a defense against bacteria and fungi [12]. TRAF3IP2 and TRAF5 were found to be involved in some autoimmune diseases like Vogt Koyanagi Harada (VKH) syndrome and Behcet's disease (BD) [13]. Through Genome-wide association (GWA) studies it was found that TRAF3IP2 is involved in psoriasis (a type of skin disease that causes redness and itchy patches). The perianal CD is a chronic disease that affects the gastrointestinal tract. Association between TRAF3IP2 and Perianal Crohn’s disease was also seen [14].

The highly expressed gene TRAF3IP2 was selected as a target to develop the treatment against osteoporosis [15]. Recently no efforts are made to make the drug by suppressing that overexpress gene. No data is available for the inhibition of these expressed genes reported in osteoporotic vertebral fracture. New drug candidates can be identified by using docking software MOE in which a library of 13000 phytochemicals will be docked against the receptor protein. In silico inhibition of over-expressed genes will pave the way for the development of a new drug candidate against this disease.

MATERIALS AND METHODS

A graphical representation of the methodology is shown in Figure 1.

Retrieval & preparation of protein

The Protein Data Bank was used to obtain the three-dimensional structure of E2 ubiquitin ligase TRAF3IP2 [PDB: 1YLA]. The MMFF94x force-field was chosen for energy minimization using the Molecular Operating Environment tool, version 2018.01. All bound water molecules and ligands were removed from the protein before docking. If explicit hydrogen, bond orders, hybridization, and charges were missing, they were assigned to the protein structure.

Ligand database preparation

Using in silico approaches, a thousand known phytochemicals were gathered from different databases, including PubChem, MPD3, and Zinc, to assess the possible inhibitory effect on E2 ubiquitin ligase TRAF3IP2 protein [16, 17]. According to the literature review, the plant-based compounds were chosen based on their medicinal capabilities [18]. Alkaloids and sterols were the most common phytochemicals chosen. The MOE tool was used to construct a fully prepared library of the compounds that were chosen [19]. Chem Draw was used to sketch the two-dimensional (2D) chemical composition of the chosen ligands. Before using the MOE ligand database, the ligands were optimized with Protonate3D and the energy was decreased.

Docking protocol

For predicting the binding affinities of a variety of ligands, molecular docking methods are commonly utilized. The current study sought to investigate the likelihood of an existing association between the experimental bioactive components of the inhibitors under investigation and the docking scores. MOE was used to locate the active pocket on the receptor protein molecule. The MOE tool was used to screen a library of 13000 phytochemicals against the interacting residues of the E2 ubiquitin ligase TRAF3IP2 protein using a molecular docking procedure. All docking experiments were run with the default parameters to obtain reliable results. The London dG scoring algorithm in MOE was used to rescore simulated poses. The nominated constraints used to calculate the interaction and score of ligand molecules with catalytic triad or compounds were; rescoring function: London dG, Rescoring 2: London dG, Retain 10, Placement: Triangle matcher, Refinement: Force field.  Based on RMSD values and S-score binding affinity, phytochemicals having top and best poses were chosen after docking. To visualize the best-docked complexes and analyze the 2D plots of ligand-receptor interactions, the MOE LigX tool was utilized. Three-dimensional pictures of protein–inhibitor complexes were also generated using MOE [19].

Evaluation of inhibitors’ druglikeness

Evaluating the drug-likeness properties of a potential drug is a crucial step in the drug development process. Molecular weight, hydrogen bond acceptors, the coefficient logP (miLogP), and hydrogen bond donors were among the physical and chemical factors evaluated. The druggability of the top-docked ligands was determined using the Molinspiration online tool (https://www.molinspiration.com (accessed on May 8, 2021)) [20].

ADMET profiling

The SwissADME and admetSAR algorithms were used to further assess phytochemical pharmacokinetic features [21, 22]. The pharmacokinetic properties of a drug include its cytotoxicity, absorption, absorption, and distribution in the human body (pharmacokinetic properties).

Energy calculations MMGBSA analysis

The binding free energy (Prime/MM-GBSA) of docking complexes was calculated using the Schrodinger Suite Release 2020 [23]. To get binding free energies, the best poses of inhibitory drugs associated with osteoporosis were chosen. In Prime, the docked complexes were minimized using the local optimization technique. The Binding free energies are calculated using Prime MM-GBSA, which combines the OPL-SAA force field, EMM (Molecular Mechanics Energies), GSGB (solvation model for polar solvation), and nonpolar solvation term (GNP), which is made up of solvent accessible surface area (SASA) and van der walls interactions.

MD simulation

MD simulations of the top two hits for the E2 ubiquitin ligase TRAF3IP2 target were done to understand the dynamics better. The simulation studies were carried out using the AMBER20 program [24]. In order to preprocess E2 ubiquitin ligase TRAF3IP2, the antechamber program was employed. The GAFF force field was employed for both targets, while the ff14SB force field was used for the enzyme [25]. The topology of enzymes and inhibitors was recorded using LEaP, with counter ions added for electrostatic neutrality.

The systems were contained within a TIP3P box [26] filled with water molecules. For 1500 steps, the sharpest descent approach was used [27]. The lengths of hydrogen bonding were limited using the SHAKE algorithm. The temperature was kept constant using the Langevin coupling integration procedure. A time step of 2 fs was used to solve Newton's equations, and the trajectory data were collected every 1 ps for further study. All MD trajectory experiments were conducted out with AmberTools20's CPPTRAJ module, and visual assessment was done with Visual Molecular Dynamics software [28].

RESULTS AND DISCUSSION

Structural and Ligand library retrieval

The three-dimensional structure of E2 ubiquitin ligase TRAF3IP2 was retrieved from the PDB database having PDB ID: 1YLA respectively. The structure was chosen as a target because of its high resolution. The resolution of 1YLA has 2.40Å. The structure was optimized and then used as a receptor. Docking was performed using the ligands library of PubChem, ZINC, and MPD3 containing 13000 phytochemicals.

Molecular docking

This section involves the results obtained through docking the receptor protein structures with the phytochemicals library using MOE software. Ten different conformations were obtained for each compound. All these compounds’ conformations were sorted based on S scoring, RSMD values, and bonding interaction with the active sites of the proteins. The top 4 compounds were selected from each receptor protein for further analysis based on the lowest S value. These selected compounds have shown strong interactions with the binding pockets of the proteins and have minimum binding energies with the scoring function of each docked ligand as shown in Table 1. Zoledronic acid was chosen as a positive control to cross-check the efficacy of our drug candidates. Zoledronic Acid Anhydrous is an antiresorptive anhydrous version of a synthetic imidazole third generation. This result in the loss of downstream metabolites required for osteoclast function, as well as the stimulation of apoptosis and, eventually, the death of osteoclast cells. Zoledronate reduces bone turnover and stabilizes the bone matrix by inhibiting osteoclast-mediated bone resorption. Although our results showed that the efficacy of our reported drug candidates was more than the standard drug. The binding affinity score for the standard drug was -7.05 kcal/mol^-1 which was far less than our reported drug candidates.

Table 1. Top drug candidates Binding affinity along with interacting residues

When E2 Ubiquitin ligase TRAF3IP2 protein was docked with the ready to dock library of 13000 phytochemicals, Adenosine,2'-deoxy- and 2-Amino-1-hydroxyoctadecan-3-one was ranked at the top due to the lowest docking score of -14.32 kcal mol-1 and -12.90 kcal mol-1 have shown potential interaction with Lys A72, Phe A75, Arg B40, Leu B31, Arg A74, Tyr A152, IIe A91, Thr A77, Lys D28, and Phe A75 as shown in Figures 2a and 2b.

 (3,5-Dihydroxyoxolan-2-yl) methyl dihydrogen phosphate had shown interaction with lleA91, Leu B31, Arg A34 and Arg B40 with the binding score of -12.80 kcal mol^-1 as shown in the Figure 3c. 2-(3,7,11,15,19,23,27,31-Octamethyldotriaconta-2,6,10,14,18,22,26,30-octaenyl) naphthalene-1,4-dione demonstrated the binding affinity of -12.66 kcal mol-1 and interaction with the sites LysA72, Lys D28 and lleA91as shown in the Figure 3d.

Mol-inspiration results

Molecular descriptors and drug likeliness properties of top compounds obtained through interactions with TRAF3IP2 were examined through the Mol-inspiration server based on Lipinski Rules of five. The results are shown in Table 2. All the compounds' values depict that they have zero volition. The low bioavailability of the compounds is a mystery. Understanding the microbiota, which works as an effective bioreactor in the human GI tract, allows for bioavailability to be redefined. The bioavailability score of all the reported inhibitors was calculated through the swissadme tool. Overall bioavailability score for all the inhibitors showed good results as mentioned in Table 2.

Table 2. Top compounds following Lipinski Rules of five

-2.05

0.55

299.5

4.57

0.55

214.11

-2.6

0.56

403.1

-3.54

0.45

Standard Drug

272.09

-4.3

0.55

Several pharmacokinetic factors were evaluated using ADME and AdmetSAR. Pharmacokinetic parameters can be used to estimate the ADME and toxicity of the top therapeutic candidate drugs. The ADMET characteristics of derived phytochemicals for both targets are shown in Table 3. Many drugs do not use this mechanism in their development due to poor pharmacokinetic properties and toxicity. To identify active lead compounds, early drug discovery relies on high-performance and rapid ADMET profiling studies.

Table 3. ADMET profiling of top compounds indicating several pharmacokinetic

MMGBSA energy calculations

The binding free energy (Delta Gbind) of inhibitory drugs with a strong potential to inhibit the activity of the osteoporosis was calculated using MMGBSA. 20ns frame was used to perform MMGBSA analysis by using the AMBER20 program. To perform the MMGBSA energy estimates, docking complexes with high energy function scores were obtained. By minimizing protein-ligand complexes, salvation energy, and surface area energy, the total free binding energy is determined. Adenosine, 2'-deoxy-/ E2 ubiquitin ligase TRAF3IP2 showed good stability with a total binding free energy of −70.04 kcal mol−1, while 2-Amino-1-hydroxyoctadecan-3-one/ E2 ubiquitin ligase TRAF3IP2 demonstrated strong binding with a binding free energy of −57.76 kcal mol−1 as mentioned in Table 4. Although the same analysis for the compounds (3,5-Dihydroxyoxolan-2-yl)methyl dihydrogen phosphate and 2-(3,7,11,15,19,23,27,31-Octamethyldotriaconta-2,6,10,14,18,22,26,30-octaenyl)naphthalene-1,4-dione in complex with the receptor protein E2 ubiquitin ligase TRAF3IP2 showed the total binding energies of −66.09 kcal mol−¹ and −88.71 kcal mol−¹.

MD simulation

A 20-ns molecular dynamic simulation was done to better understand the molecular mechanisms involved in the binding of Adenosine, 2'-deoxy, 2-Amino-1-hydroxyoctadecan-3-one, (3,5-Dihydroxyoxolan-2-yl)methyl dihydrogen phosphate and the ligand 2-(3,7,11,15,19,23,27,31-Octamethyldotriaconta-2,6,10,14,18,22,26,30-octaenyl)naphthalene-1,4-dione in the active pocket of the 1YLA protein.

RMSD analysis

We were able to determine the stabilization of the docking complex (Adenosine, 2'-deoxy- and the ligand 2-Amino-1-hydroxyoctadecan-3-one) by measuring the root mean squared deviations (RMSD) during the production run. The complex Adenosine, 2'-deoxy- root mean square deviation (RMSD) indicated a moderate fluctuation of roughly 0.2 Å up to 20ns and a constant trajectory during the production run, within 2.5 Å to 3Å as shown in Figure 4a. While the second ligand 2-Amino-1-hydroxyoctadecan-3-one showed stability throughout the time period of 20 ns as shown in Figure 4b.

RMSF analysis

During the simulation interval, the atoms in the Adenosine, 2'-deoxy-, and ligand 2-Amino-1-hydroxyoctadecan-3-one complexes displayed acceptable and steady RMSF variations, with the ligand contact maximum at 4.8 and 4.5, respectively (Figures 4c and 4d). This suggests that these proteins and their ligands have developed a very stable complex system. (3,5-Dihydroxyoxolan-2-yl)methyl dihydrogen phosphate complex showed minor fluctuation up to residue number 200.

Solvent accessible surface area (SASA)

SASA's major trait is that it keeps macromolecules stable and folded. The SASA values for the mutant and wild forms have been determined. The average SASA values for Adenosine, 2'-deoxy-/ E2 ubiquitin ligase TRAF3IP2 complex were 50 to 150 while 2-Amino-1-hydroxyoctadecan-3-one/ E2 ubiquitin ligase TRAF3IP2 showed stability between 30 to 60 residue number 152, respectively, suggesting that throughout the simulation phase, no variations in the area covered by all the systems were noticed (Figures 4e and 4f). (3,5-Dihydroxyoxolan-2-yl)methyl dihydrogen phosphate showed stability between the residues number 150 to 200.

Conventional drug designing is a time-consuming process, a long time is needed to invent a novel drug by this approach. Computer-based drug designing methods change the drug-designing methods dramatically and made the drug in a very short duration with the least side effects. In-silico approaches are not only less time-consuming but also cost-effective as compared to conventional drug design [29, 30]. With the advancements in the field of Bioinformatics, different tools and algorithms have been advised that identify and discover the drug targets [31]. Molecular docking is one of the popular and modern in silico approaches in the field of medical research that led to the invention of new compounds by observing the pattern of the small molecule’s binding to the target and hence discovering the new drug that combats the deadly diseases [32, 33].

The present study illustrated that by identifying the active sites of the target proteins, we can inhibit their expression. Some compounds have considerable interactions with these three proteins involved in osteoporosis. The molecular properties and drug-likeliness of the selected complexes were estimated according to the "Lipinski Rule of Five". This rule states that the molecular weight of the compound must be less than 500 Daltons, less than 5 Hydrogen bond donors, no more than 10 Hydrogen bond acceptors, and an AlogP value of fewer than 5. All compounds Adenosine, 2'-deoxy-, 2-Amino-1-hydroxyoctadecan-3-one, (3,5-Dihydroxyoxolan-2-yl)methyl dihydrogen phosphate and 2-(3,7,11,15,19,23,27,31-Octamethyldotriaconta-2,6,10,14,18,22,26,30 -octaenyl)naphthalene-1,4-dione fulfill the Lipinski's Rule of Five and show no violation. Selected compounds have low scoring values and have RMSD values less than 3. ADMET analysis is a challenging process in drug discovery. This is achieved through the SwissADME database and showed that selected compounds have good pharmacokinetic properties. The compounds were then evaluated for blood-brain barrier (BBB) penetration, HIA (Human Intestinal Absorption), and AMES monitoring. Predicting the ADMET characteristics is an important predictor of the drug candidate's behavior, toxicity level, and fate in the human body [34, 35]. It indicates the possibility of the candidate's capacity to enter the intestinal absorption, metabolism, blood-brain barrier, subcellular localization, and, most importantly, the potential for injury to the body. Many drugs that do not go through the process of drug development are due to toxicity and poor pharmacokinetic properties. Identifying active lead compounds at early drug discovery is facilitated by producing high-performance and quick ADMET profiling assays.  To cross-check the stability of performed analysis top ligands along with the targeted receptor having the least docking score were chosen and the MD simulation analysis was performed. To further analyze the stability of the complexes. RMSD, RMSF, and SASA values were calculated and they showed that the compounds showed overall good stability with the receptor protein. So it is concluded that the extracted phytochemicals can inhibit the activity of the proteins involved in osteoporosis by targeting their binding pockets and hence they can be used as effective drug candidates against the disease. The efficacy of the reported candidates was cross checked by comparing the results with already available drug

CONCLUSION

This study aims to identify the novel compounds that can be used as effective drugs against this highly over-expressed gene of osteoporosis. A Revolution in the field of Bioinformatics has advised different tools and algorithms that lead to the drug target identification. MOE is the latest tool in the field of bioinformatics for drug discovery. Proteins of these genes were isolated and docked with the ready to dock library of compounds using MOE. The resulted compounds Adenosine, 2'-deoxy-, 2-Amino-1-hydroxyoctadecan-3-one, (3,5-Dihydroxyoxolan-2-yl)methyl dihydrogen phosphate and  2-(3,7,11,15,19,23,27,31-Octamethyldotriaconta-2,6,10,14,18,22,26,30 -octaenyl)naphthalene-1,4-dione were analyzed using all the parameters like Lipinski rule of five, ADMET assay etc. ADMET analysis brings forth its absorption in the intestine, distribution, the tendency of a compound to pass through the blood-brain barrier, metabolism, and the toxicity and damage it can produce in the body. After using different bioinformatics tools, it is observed that the selected compounds have the potential for drug-likeness properties. Adenosine, 2'-deoxy- and 2-Amino-1-hydroxyoctadecan-3-one were further validated by performing energy calculations and MD simulation analysis. All of the analyses showed stability throughout the 20 ns period. The current study results might help invent and design the novel compounds that have the better inhibitory activity against the overexpressed proteins of up-regulated genes. It is greatly advisable to examine the effectiveness of the phytochemicals that are used in this research to study for additional studies.

ACKNOWLEDGMENTS : None

CONFLICT OF INTEREST : None

FINANCIAL SUPPORT : None

ETHICS STATEMENT : None