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

 

Bandar Hamad Aloufi1*

1Department 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.

 

Figure 1. Flow Chart a Graphical representation of methodology.

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

Compounds ID

Compounds Name

2D Structures

Binding Affinity

RMSD

Interacting Residues

636

Adenosine, 2'-deoxy-

-14.32 kcal/mol-1

1.95

Lys A72

Phe A75

Arg B40

Leu B31

Arg A74

Tyr A152

631

2-Amino-1-hydroxyoctadecan-3-one

-12.90 kcal/mol-1

1.94

IIe A91

Thr A77

Lys D28

Phe A75

 

635

(3,5-Dihydroxyoxolan-2-yl)methyl dihydrogen phosphate

-12.80 kcal/mol-1

1.64

IIe A91

Arg B40

Leu B31

Arg A74

 

633

2-(3,7,11,15,19,23,27,31-Octamethyldotriaconta-2,6,10,14,18,22,26,30-octaenyl)naphthalene-1,4-dione

-12.66 kcal/mol-1

1.49

Lys A72

IIe A91

Lys D28

 

Standard Drug

68740

Zoledronic acid

 

-7.05 kcal/mol-1

2.09

Arg A74,

Arg B40

Leu B31

 

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.

a)

b)

Figure 2. Three-dimensional representation of molecular docking analysis between top two compounds with interacting residues. a) Adenosine, 2'-deoxy-/ E2 ubiquitin ligase TRAF3IP2 b) 3,5-Dihydroxyoxolan-2-yl)methyl dihydrogen phosphate/ E2 ubiquitin ligase TRAF3IP2.

 (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.