Molecular docking based screening of Noggin inhibitors

Noggin (NOG) a BMP (bone morphogenetic protein) antagonist plays a key role in preferentially driving a subset of breast cancer cells towards the bone and causing osteolytic lesions leading to severe pain and discomfort in the patients. Owing to its role in bone metastasis, NOG could be promising molecular target in bone metastasis and that identifying small molecule inhibitors could aid in the treatment. Towards identifying cognate inhibitors of NOG, structure based virtual screen was employed. A total of 8.5 million ligands from e-molecule database were screened at a novel binding site on NOG identified by the Sitemap tool, employing GLIDE algorithm. Potential eight molecules were selected based on the Glide score, binding mode and H-bond interactions. Free energy of binding was calculated using Molecular mechanics based MMGBSA and the obtained energy was used in the prioritizing the compounds with the similar structures and glide score. Further, the compounds were evaluated for their druggability employing physico-chemical property analysis. Our study helped in identifying novel potential NOG inhibitors that can further be validated using in-vivo and in-vitro studies and these molecules can also be employed as tool compounds to study the functions of BMP.


Background:
Early diagnosis of breast cancer is pivotal in the maximizing the survival rates of the cancer patients. Often, breast cancers are detected only after they are metastasized. One of the major metastatic sites of the breast cancer is the bone [1]. Bone metastasis leads to pathological fractures, life threatening hypercalcemia, spinal cord compression, severe pain and morbidity.
Understanding, the underlying molecular mechanisms in bone metastasis helps in identifying plausible novel targets, which could ameliorate pain and reduce morbidity.
Bone tissue is made up of osteoblasts, osteoclasts and osteocytes. Osteoblasts are involved in the bone formation, while osteoclasts in the re-sorption of the bone. RANKL (Receptor activated NF kappaBLigand) is a member of the tumor necrosis factor cytokine family and is responsible for osteoclast differentiation and activation. OPG (Osteoprotegrin) is an osteoblast-secreted decoy receptor that functions as a negative regulator of bone re-sorption. Always equilibrium is maintained between the RANKL and OPG. Shift of this equilibrium towards the RANKL results in lesions that destruct the bone conversely, shift towards OPG results in bone formation, which is brought about by the family of growth factors called Bone morphogenetic proteins (BMP) [2,3]. Tumor  NOG is a secreted glycosylated homodimer and acts by directly binding to the BMP and preventing BMPs from binding to their receptors. NOG is preferentially expressed in the breast cancer cells that metastatize to the bone. It is involved in the numerous developmental processes. Binding of NOG to BMPs shifts the equilibrium between the RANKL and OPL towards RANKL there by resulting inosteolytic lesions [6].
Recent evidence suggests that NOG plays a significant role in the tumor growth and progression. Keratin 14-driven NOG over expression in mice results in development of skin tumors [7]. The osteolytic lesions in bones xenografted with the PC3 (human prostate cancer cell line) cells showed increased osteoclast activity and reduced osteoblast activity. Interestingly, when NOG-silenced PC3 cells were used repair activity was seen in lesions emphasizing the role of NOG in prostate cancer [8]. Expression of NOG in breast cancer cells provides them with bone colonization capabilities and also increased osteoclast activity and when NOG was silenced the osteoclast activity was reduced [9]. From these results we hypothesized that NOG inhibition could help in reducing bone metastatic cancer progression thereby alleviating pain in the metastatic bone lesions.
Previous studies by Karen etal., identified flavonoids that activate the BMP signaling pathway by inhibiting NOG [10]. Here for the first time we intended to identify small molecule inhibitors of NOG using structure based virtual screening that would possibly increase the available BMP levels, thereby may aid in restoring the bone damage and thus inhibit bone metastatic cancer progression. Alternatively, some of these molecules can be used, as tool compounds that would help to further understand the functions of NOG and BMPs in the context of various cancers. In order to achieve the above-mentioned objective we employed high throughput SBVS of small molecules.

Methodology:
Protein preparation: Structure of the NOG was retrieved from PDB with the identification number 1M4U [11]. Loops missing in the PDB structure were modeled using SwissModel (https://swissmodel.expasy.org) [12]. To ensure correct starting structures initial structure of the protein was refined and subjected to energy minimization. The 3D model of the protein was prepared using the Protein Preparation Wizard in Maestro [13]. Protein was prepared by adding the hydrogen atoms, optimizing hydrogen bonds, removing atomic clashes, adding formal charges to the hetero groups and then optimizing at neutral pH. Finally the structure was minimized using optimized potential for liquid simulations (OPLS-2005) force field.

Binding site analysis:
Binding site was generated using the SiteMap tool from schrodinger [14]. SiteMap characterizes sites through a combination of properties calculated at each site point which include size of the site as measured by number of site points, the degree of enclosure by protein, degree of exposure to solvent, tightness with which sitepoints interact with the protein, hydrophobic and hydrophilic character of the site as well as the balance between these terms and degree to which a ligand and has possibility of donating or accepting hydrogen bonds. An overall SiteScore and druggability score are then used in selection of potential binding sites.

Ligand and Library preparation:
The Ligand and molecules from e-molecule database are retrieved in structure data file (SDF) format. The molecules were subjected to ligand preparation using Ligand Prep module of Schrödinger suite (Schrödinger). The ligand preparation process of molecules involves preserving the definite chiralities, to generate minimum five low-energy stereoisomers per ligand, using default conditions at pH 7.0 ± 2.0. The resulting ligands are subjected to high throughput virtual screening using the GLIDE (Grid based Ligand and Docking with Energetics) module of Schrödinger suite.

Docking:
Receptor grid was generated using receptor grid generation in the Glide application (Glide, version) of Maestro. Once receptor grid was generated, ligands were docked to the protein-using Glide docking protocol. The ligands were docked by using a three tire docking which starts with "High throughput Virtual Screening" (HTVS) followed by "Standard Precision" (SP) and then by "Extraprecision" mode (XP). The docked conformers were evaluated using Glide (G) Score [15,16].

MMGBA:
The docked complexes were subjected to Molecular Mechanics/Generalized Born Surface Area (MM-GBSA) analysis for predicting the binding energy by prime approach. The binding energies, obtained through MM-GBSA OPLS-2005 are considered much more accurate than the XP GScore [17]. ADME property analysis: ADME properties of selected ligands were analyzed using QikProp tool of Schrodinger suite. The tool predicts physiochemical properties with a detailed analysis of: (i) Molecular Weight (ii) partition coefficient (iv) hydrogen bond donors (v) hydrogen bond acceptors (vi) number of rotatable bonds

Results: Protein preparation and binding site analysis:
NOG small molecule antagonist could be used locally to promote bone formation. In order to identify suitable binding pocket for virtual screening of small molecule inhibitors; binding site analysis was performed using SiteMap tool of Schrodinger software. SiteMap identified five sites based on site score that includes size, volume, amino acid exposure, enclosure, contact, hydrophobicity, hydrophilicity and donor/acceptor ratio. The site having site score and druggability score of 0.95 and 0.98 respectively was used for further screening of the novel compounds using docking protocol. The sites with site score of 1 and above can be a suitable site for the ligand binding. In addition druggability score of >1 is

Small molecule Database preparation:
The eMolecule database containing 7303475 compounds was used and filtered using Rapid Elimination Of Swill (REOS) which is a hybrid method that combines a set of functional group filters with some simple counting schemes with more than 200 rules [19]. A total of 5744923 compounds were obtained after REOS filtering and subsequently passed through the Pan-Assay Interference Compounds (PAINS) to remove non-specific compounds [20]. This filter further reduced the number to 5734778 compounds and these were clustered using molprint2D finger print with 64 bit and with Leader-follower clustering method using Canvas tool of Schrodinger software to obtain diverse compounds. Finally compounds were prepared using Ligand Preparation tool with pH 7 .0 (+/-1.0). Ligand preparation generated 8537456 compounds (maegz format) with proper ionization and tautomeric states were generated and are used for virtual screening.

Virtual Screening and Identification of small molecule inhibitors:
To identify small molecule inhibitors we have performed SBVS using prepared 8.5 million druggable compounds of eMolecule database against the selected site. Before running virtual screening a grid was generated on protein binding site using selected sitemap site to calculate potentials of binding site that are used for docking and scoring the compounds. The docking based screening was performed multi-tiered screening protocol, starting with HTVS followed by SP and XP methods [16]. The High Throughput Virtual Screening (HTVS) mode is designed to screen large libraries quickly with rough scoring functions, hence to screen 8.5 million compounds we started with this method. The top ranked hits (top 20%) were passed through Standard Precision (SP) mode, which is ten times slower and more precise than HTVS. The SP method is more exhaustive in conformational sampling and more precise than HTVS method with the expense of time. About 20,000 compounds obtained from SP method (best 50% of the compounds) were further evaluated with even more precise and more computationally intensive Extra Precision (XP) method. About 1000 compounds obtained from XP method were shortlisted based on docking score that are -9.0 and above. The high glide score indicated a high binding affinity towards the target. We checked for the following interactions, hydrogen bonds, salt bridges, halogen bonds, aromatic bonds, pi-cation and also pi-pi interactions all of which contribute towards the stability of the protein-ligand complexes. All the compounds

Molecular Mechanics Genralized Born Surface Area (MMGBSA):
The ΔG bind between NOG and each ligand, respectively, were calculated using the MMGBSA and are shown in Figure 1. Ligand 7 showed the maximum ΔG bind and the ligand 8 showed the minimum ΔG bind.

ADMET Analysis:
In order to understand whether these compounds could further be utilized as drug molecules we performed ADMET analysis on the calculated ADMET properties. The analysis shows that all the selected molecules follow Lipinski's rule of five and number of rotatable bonds <10 with an exception of ligand 7, molecular weight and other parameters do not follow Lipinski's rule. Similarly ligand 8 is also predicted to have low oral absorption. The results are shown in Table 1.

Discussion:
Based on the similarity of the structure, molecules were grouped into five sets. Ligand 1 and Ligand 2 in set 1 has a substituted thiazole attached to ethylamide linker, which in turn is connected to pyridine ring. Amide -NH in the ligand 1 forms a NH-O interaction with backbone carbonyl atom of Pro 223. The pyridone carbonyl on ligand 1 and 2 interact with the backbone -NH of Ile 225 forming a NH-O interaction. Amine group in the pyridone on the ligand 1 and 2 forms NH-O bond with the backbone carbonyl atom of the Gly 176. Additionally, pyridone carbonyl on ligand 1 also forms a NH-O with the Gly176 resulting in a bifurcated H-bond. Since, each strong H-bond accounts to about 15 kcal of enthalpic stabilization bifurcated Hbond interactions are expected to contribute more stabilization in comparison with a single H-bond interaction. Ligand 1 also forms CH-O interaction with the backbone carbonyl atom of Pro 223. Next calculating the ΔG bind values helps in the prioritizing the potential compounds. Ligand 2 with a ΔG bind value of -46kcal/mol and ligand 1 with a ΔG bind value of -39 kcal/mol is throwing a different picture in contrary to the earlier mentioned contribution of bifurcated H-bonds in ligand1. In this context, if molecules are short-listed just based on docking scores without binding energy might show a different order of priority in comparison when we judge them with binding free energy contribution. Therefore we prioritized ligand 2 in set 1 vis-á-vis ligand 1 from the same set. ring formed NH-O interaction with the backbone carbonyl group of Gly176. The suphonyl oxygen forms hydrogen bond with the backbone amine group of Asp 45. A good glide score and the ΔG bind of -9.02 and -43.280 kcal/mol respectively and also abides to the Lipinski rule of 5, indicating that it could be a potential drug candidate. The ligand 7 in the set 4 had good glide score of -9.038 and binding energy of -80.995 kcal /mol, shows strong hydrogen bond interactions but owing to its poor physicochemical properties cannot be considered as a potential drug candidate. Ligand 8 although had a good glide score of -9.73 and binding energy was very high -26.20kcal /mol. Therefore ligand 8 hence is not a good candidate. In the light of the docking scores, hydrogen bond, binding energy and physicochemical properties we conclude that ligand 2, ligand 6 and ligand 1 to be potential Noggin inhibitors.

Conclusion:
In summary, we used structure based virtual screening to identify novel small molecule inhibitors of NOG. We first ranked the compounds based on their glide score and binding mode. The compounds with the glide score greater than -9.0 were further analyzed for the hydrogen bond interactions. Our analysis indicated that the all the ligands interacted with the backbone carbonyl oxygen of the Gly 176 of β2 region. Pro223 and Ile225 are the key residues involved in the binding. Further, we utilized binding free energy for the selection of potential inhibitors. ADMET analysis was used to understand the physic chemical properties and drug likeliness of the compounds. Three compounds, Ligand 2, Ligand 6 and ligand 1 are found to be potential NOG inhibitors that had a good glide score, binding energy and had favorable physiochemical properties so as to be considered as drugs. These molecules can be used in the prevention and as adjuvant therapy in bone metastasis. Additionally, these can also be used, as tool compounds to further understand the function of the NOG and BMP.