Discovery of Some Heterocyclic Molecules as Bone Morphogenetic Protein 2 (BMP-2)-Inducible Kinase Inhibitors: Virtual Screening, ADME Properties, and Molecular Docking Simulations

Bone morphogenetic proteins (BMPs) are growth factors that have a vital role in the production of bone, cartilage, ligaments, and tendons. Tumors’ upregulation of bone morphogenetic proteins (BMPs) and their receptors are key features of cancer progression. Regulation of the BMP kinase system is a new promising strategy for the development of anti-cancer drugs. In this work, based on a careful literature study, a library of benzothiophene and benzofuran derivatives was subjected to different computational techniques to study the effect of chemical structure changes on the ability of these two scaffolds to target BMP-2 inducible kinase, and to reach promising candidates with proposed activity against BMP-2 inducible kinase. The results of screening against Lipinski’s and Veber’s Rules produced twenty-one outside eighty-four compounds having drug-like molecular nature. Computational ADMET studies favored ten compounds (11, 26, 27, 29, 30, 31, 34, 35, 65, and 72) with good pharmacokinetic profile. Computational toxicity studies excluded compound 34 to elect nine compounds for molecular docking studies which displayed eight compounds (26, 27, 29, 30, 31, 35, 65, and 72) as promising BMP-2 inducible kinase inhibitors. The nine fascinating compounds will be subjected to extensive screening against serine/threonine kinases to explore their potential against these critical proteins. These promising candidates based on benzothiophene and benzofuran scaffolds deserve further clinical investigation as BMP-2 kinase inhibitors for the treatment of cancer.


Introduction
Bone morphogenetic proteins (BMPs) are a type of transforming growth factor-β (TGF-β) superfamily that has multiple functions. More than a quarter of a million BMP members have been identified. Cell proliferation, survival, differentiation, and apoptosis are all regulated by BMPs. It is important to produce bone, cartilage, ligaments, and tendons [1].
When BMPs attach to their cell surface receptors on mesenchymal cells, the BMP signaling cascade is initiated, and signals are transmitted to the cell nucleus via particular proteins. The mesenchymal cell becomes a chondrocyte or an osteoblast as a result of the expression of genes that lead to the production of macromolecules involved in cartilage and bone development [2]. Following the implantation of this protein component of the bone matrix, a complex series of cellular events occurred, including mesenchymal cell infiltration, cartilage formation, vascularization, bone formation, and finally remodeling of new bone tissue, as well as population by hematopoietic bone marrow elements [3]. BMPs have been found to directly differentiate cells into the osteoblast phenotype, in addition to chondrocyte lineage differentiation [4]. Bone morphogenetic protein-2 (BMP-2) is the only osteoinductive growth factor approved by the FDA for use as a bone graft alternative.
BMP signaling is associated with cancer-related cellular phenotypes. BMPs can increase cell migration and invasiveness by arresting the cell cycle of many different cell types in the early G1 phase [5]. BMPs have been shown to play a function in a wide range of cancers and other malignancies. BMPs can either suppress or promote tumorigenesis, with the majority of cases favoring metastasis [6]. Intestinal tumorigenesis can be triggered by disruption of the tightly controlled homeostatic BMP signaling gradients [7]. The most important factor of hereditary risk and predisposition for sporadic colorectal cancer (CRC) susceptibility, according to recent research, is variation in the BMP pathway [8]. As rhBMP-2 has been observed to be associated with a higher risk of developing new cancer than vertebral bone graft, its use in spine surgery has been the subject of significant discussion [9]. BMP-2-inducible kinase belongs to the Numb-associated kinase (NAK) family of serine/threonine kinases [10]. BMP-2 is crucial to the occurrence and progression of colon cancer, prostatic carcinoma, and lung cancer [11,12]. Furthermore, BMP-2 is essential for prenatal and postnatal mammary gland development, but it has also been identified in breast cancer cells [13]. Additionally, BMP-2 increases gastric cancer cell motility and invasion by activating PI-3 kinase/Akt; blocking this pathway may limit BMP-2-mediated metastasis [14]. Various heterocyclic molecules are found in a variety of drugs and have become a key research foundation in medicinal chemistry. This is mostly owing to the adaptability and specific physicochemical properties of heterocyclic molecules. Benzofuran is one of the identified heterocyclic compounds [15]. Benzothiophene and Benzofuran scaffold is one of the privileged frameworks in drug development, as this core displays diverse biological activities allowing them to function as anti-convulsant, anti-cancer, antidiabetic, anti-tubercular, anti-oxidant, anti-inflammatory, anti-microbial, and many other agents [16,17].
Many literature studies revealed that substituted benzothiophene and benzofuran derivatives were a new class of small molecules that act as potential anabolic agents targeting BMP-2 [18,19]. In particular, compounds I, II III, and IV ( Figure 1) enhanced BMP-2 expression in vitro and stimulated bone formation and trabecular connectivity restoration in vivo [18]. In contrast, benzothiophene and benzofuran derivatives can inhibit several protein kinases and act as anticancer agents such as compounds V-VIII ( Figure 2) [16,[20][21][22].
In this work, we will try to find new inhibitors for BMP-2 that may be useful in the treatment of several types of cancers. We decided to investigate some reported benzothiophene and benzofuran derivatives [23,24] for their potential as BMP-2 inducible kinase inhibitors. The selected library (Figure S1 Supplementary data) was subjected to different computational techniques such as ADMET, toxicity, and docking studies to reach promising candidates with proposed activity against BMP-2-inducible kinase. Finally, the promising compounds will undergo additional screening against kinases related to serine/threonine kinases (CDK2, Pim1, cell division protein kinase 2, casein kinase II, and  eukaryotic translation initiation factor 2-alpha kinase 3).  In this work, we will try to find new inhibitors for BMP-2 that may be useful in the treatment of several types of cancers. We decided to investigate some reported benzothiophene and benzofuran derivatives [23,24] for their potential as BMP-2 inducible kinase inhibitors. The selected library (Figure S1 Supplementary data) was subjected to different computational techniques such as ADMET, toxicity, and docking studies to reach promising candidates with proposed activity against BMP-2-inducible kinase. Finally, the promising compounds will undergo additional screening against kinases related to serine/threonine kinases (CDK2, Pim1, cell division protein kinase 2, casein kinase II, and eukaryotic translation initiation factor 2-alpha kinase 3).  In this work, we will try to find new inhibitors for BMP-2 that may be useful in the treatment of several types of cancers. We decided to investigate some reported benzothiophene and benzofuran derivatives [23,24] for their potential as BMP-2 inducible kinase inhibitors. The selected library (Figure S1 Supplementary data) was subjected to different computational techniques such as ADMET, toxicity, and docking studies to reach promising candidates with proposed activity against BMP-2-inducible kinase. Finally, the promising compounds will undergo additional screening against kinases related to serine/threonine kinases (CDK2, Pim1, cell division protein kinase 2, casein kinase II, and eukaryotic translation initiation factor 2-alpha kinase 3).

Virtual Screening
In the current study, an in silico computational study was conducted to determine the number of rotatable bonds, topology polar surface area (TPSA), and other physicochemical properties for the tested candidates according to the directions of Veber's and Lipinski's Rules of five [25].
Lipinski proposed that the absorption of an orally administered compound is more likely to be better if the molecule satisfies at least three out of four of the following rules: (1) hydrogen bond donors (HBD) ≤ 5; (2) hydrogen bond acceptors (HBA) ≤ 10; (3) molecular weight < 500; (4) a coefficient of partition between octanol and water logP < 5. Compounds contravening more than one of these rules could not have good bioavailability. Moreover, reduced molecular flexibility, as measured by the number of rotatable bonds, and low polar surface area are found to be important predictors of good oral bioavailability [26]. In this regard, Veber's Rule says that a compound with 10 or fewer rotatable bonds (RTB) and a polar surface area (TPSA) no greater than 140 Å2 should present good oral bioavailability [26].
The results presented in Table 1 showed that twenty-one outside eighty-four tested compounds (11,(26)(27)(28)(29)(30)(31)(32)(33)(34)(35)40, 61-63, 65, and 68-72) showed no contravention of Lipinski's and Veber's Rules, respectively, and hence display a drug-like molecular nature.  The obtained values match with the fourth parameter of Lipinski's Rule, which suggests that good drug candidates should have logP values less than 5. According to Veber's Rule, number of rotatable bonds was calculated, and this is an important parameter to measure the molecular flexibility and oral bioavailability of the drug candidates. The results revealed that compounds 11, 26-35, 40, 61-63, 65, and 68-72 displayed an acceptable number of rotatable bonds ranging from 2-10 that meet the criteria of Veber's Rule. Moreover, the number of TPSA (a physicochemical property describing the polarity of molecules) of such compounds is within the acceptable values (ranging from 74.77 to 123.91 Å 2 ) of less than 140 Å 2 . On the other hand, the rest of the tested compounds did not meet the criteria of Lipinski's and Veber's Rules, respectively. A closer look at the data presented in Table S1 Molecules 2022, 27, 5571
The results revealed that ten out of twenty-one compounds showed good ADMET profiles and drug-likeness properties. In detail, compounds 11, 26, 27, 29, 30, 31, 34, 35, 65, and 72 exhibited medium to very low BBB penetration levels ( Figure 3); therefore, such compounds were anticipated to be safe for the CNS. In addition, it was found that all of the ten compounds had moderate to good absorption behavior. Moreover, the solubility level of the ten compounds is expected to be better than or even similar to that of the reference drug that showed a low solubility level. For cytochrome P450 2D6 (CYP2D6) inhibition, all the examined members were predicted as non-inhibitors. Finally, all compounds were expected to bind to plasma protein by more than 90% (Table 2).
All the tested compounds showed TD 50 values better than IDK, and the values ranged from 6.036 to 431.948 g/kg body. Furthermore, compounds 27, 29, 30, and 31 showed rat maximum tolerated dose values higher than that of IDK, and the values range is from 0.190 to 0.210 g/kg body weight. Conversely, the rest of the compounds showed maximum tolerated doses ranging from 0.035 to 0.115 g/kg that were lower than that of IDK. Compounds 34, 35, and 65 revealed rat oral LD 50 values of 2.364, 5.598, and 5.308, respectively, which were higher than that of IDK (2.626 g/kg body weight). On the other hand, the rest of the compounds showed oral LD 50 values ranging from 0.230 to 2.223 g/kg body weight, which were lower than that of IDK. Additionally, all compounds were estimated to be non-toxic against the developmental toxicity potential model except compound 34. Moreover, compounds 26, 27, 34, and 72 showed LOAEL values ranging from 0.053 to 0.088 g/kg body weight whereas IDK exhibited 0.049 g/kg body weight. Finally, all the elected compounds showed no irritancy in the skin irritancy model as shown in Table S2 (Supplementary Data).

Molecular Docking Studies
Molecular docking was performed using Molecular Operating Environment (MOE) [35] to recognize the binding modes and interactions with the crucial amino acids [36][37][38][39][40][41]. Nine compounds (11, 26, 27, 29, 30, 31, 35, 65, and 72) displayed good computational ADMET and toxicity profiles. They were subjected to further docking studies. Such studies were carried out for the tested candidates to inspect their binding free energies (∆G) and binding modes against BMP-2-inducible kinase (PDB ID: 5I3R) using IDK as a reference (Table 3). At first, the docking procedure was validated and the RMSD value was 1.05, which indicated the validity of the docking process ( Figure S2) (See Supplementary data).
IDK orientation with the amino acids of the pocket has been studied and displayed in Figure 4. The proposed binding pattern revealed an affinity value of −22.72 Kcal/mol. The interaction between IDK and BMP-2-inducible kinase is stabilized through five hydrogen bonds and a series of hydrophobic interactions. The N-(1H-indazol-3-yl)cyclopropanecarboxamide moiety forms three hydrogen bonds with the two crucial amino acids, Cys133 and Glu131, and nine hydrophobic interactions with Cys133, Leu57, Leu187, Val65, and Ala77. Moreover, 1-cyclopropyl-N-phenylmethanesulfonamide moiety was buried in the active pocket of BMP-2-inducible kinase through the formation of two hydrogen bonds with Gln137 and Asn185 and four hydrophobic interactions with Val65, Ala58, and Lys79 ( Figure 4).
Compound 26 (affinity value of −22.90 Kcal/mol) combined with the receptor protein in a manner similar to that of IDK. The 3-methoxybenzo[b]thiophene moiety formed two hydrogen bonds with Cys197 and Gln137, in addition to seven hydrophobic interactions with Val65, Cys197, Ala58, Ala77, and Leu187, and two pi-sulfur interactions with Met130. Furthermore, the 2-carboxylic acid moiety formed two hydrogen bonds with the two essential amino acids, Cys133 and Glu131 ( Figure 5). drogen bonds and a series of hydrophobic interactions. The N-(1H-indazol-3-yl)cyclopropanecarboxamide moiety forms three hydrogen bonds with the two crucial amino acids, Cys133 and Glu131, and nine hydrophobic interactions with Cys133, Leu57, Leu187, Val65, and Ala77. Moreover, 1-cyclopropyl-N-phenylmethanesulfonamide moiety was buried in the active pocket of BMP-2-inducible kinase through the formation of two hydrogen bonds with Gln137 and Asn185 and four hydrophobic interactions with Val65, Ala58, and Lys79 (Figure 4).  Docking simulation of compound 27 revealed that it has a good fitting into the enzyme active site with a docking score of −23.34 Kcal/mol. The benzo[b]thiophene-2-carboxylic acid moiety formed two hydrogen bonds with Cys133 and Cys197. In addition, it formed seven pi-alkyl and pi-sigma interactions with Ala77, Val65, Cys197, Ala58, and Leu187. Furthermore, the 2-amino-2-oxoethoxy moiety formed two hydrogen bonds with Gln137 ( Figure 6).
Interestingly, the docking result of compound 31 (affinity value of −22.22 Kcal/mol) is almost like that of IDK. The 3-ethoxybenzo[b]thiophene-2-carboxylic acid moiety formed three hydrogen bonds with Cys133 and Glu131. In addition, it is involved in binding with the receptor through the formation of seven hydrophobic interactions with Ala77, Ala58, Cys197, Leu187, and Val65, in addition to two pi-sulfur interactions with Met130. On the other hand, the 2-nitrobenzene moiety formed one pi-pi interaction with Leu57 ( Figure 7).

Further Investigation against Protein Kinases
The promising nine compounds will be subjected to further screening against kinases in general to explore their potential against these crucial proteins. Compound 11 was revealed to target 48 types of kinases (Supplementary Data) with a binding affinity ranging from −4.694 to −10.033 Kcal/mol. Compound 11 showed a ligand similarity score of 0.428 with Cell division protein kinase 2 ligand (A27). Moreover, it showed a binding similarity score of more than 66% with a docking score equal to −9.463 Kcal/mol. Figure 10 illustrates the mode and interactions of compound 11 with Cell division protein kinase 2. Compound 26 (figure 11) revealed good potential to target three types of kinases, Casein kinase II subunit alpha, Proto-oncogene serine/threonine-protein kinase pim-1, and Serine/threonine-protein kinase pim-1 with a ligand similarity score of more than 40% and a binding affinity range from −6.4 to −7.00 Kcal/mol. Compound 27 showed the ability to target seven types of kinases with a ligand similarity score of more than 40% and a binding affinity range from −6.4 to −6.7 Kcal/mol. The best binding affinity was assigned for Serine/threonine-protein kinase pim-1. Furthermore, compound 27 revealed a good potential to target seven kinases with ligand similarity of more than 40% and a binding affinity range from −6.7 to 8.7 Kcal/mol, and it has shown a ligand and binding similarity mode of more than 40% and a docking score energy equal to −8.7 Kcal/mol with Cyclin-dependent kinase 2. Compound 30 showed potential against nine types of protein kinases, and the best proposed potential was assigned for cell division protein kinase 2 with a binding

Further Investigation against Protein Kinases
The promising nine compounds will be subjected to further screening against kinases in general to explore their potential against these crucial proteins. Compound 11 was revealed to target 48 types of kinases (Supplementary Data) with a binding affinity ranging from −4.694 to −10.033 Kcal/mol. Compound 11 showed a ligand similarity score of 0.428 with Cell division protein kinase 2 ligand (A27). Moreover, it showed a binding similarity score of more than 66% with a docking score equal to −9.463 Kcal/mol. Figure 10 illustrates the mode and interactions of compound 11 with Cell division protein kinase 2. Compound 26 ( Figure 11) revealed good potential to target three types of kinases, Casein kinase II subunit alpha, Proto-oncogene serine/threonine-protein kinase pim-1, and Serine/threonine-protein kinase pim-1 with a ligand similarity score of more than 40% and a binding affinity range from −6.4 to −7.00 Kcal/mol. Compound 27 showed the ability to target seven types of kinases with a ligand similarity score of more than 40% and a binding affinity range from −6.4 to −6.7 Kcal/mol. The best binding affinity was assigned for Serine/threonine-protein kinase pim-1. Furthermore, compound 27 revealed a good potential to target seven kinases with ligand similarity of more than 40% and a binding affinity range from −6.7 to 8.7 Kcal/mol, and it has shown a ligand and binding similarity mode of more than 40% and a docking score energy equal to −8.7 Kcal/mol with Cyclin-dependent kinase 2. Compound 30 showed potential against nine types of protein kinases, and the best proposed potential was assigned for cell division protein kinase 2 with a binding affinity of −8.7 Kcal/mol. Figure 12      Compound 31 showed binding affinity equal to −8.53 against Cyclin-dependent kinase 2 and compound 35 showed more binding affinity towards protein kinase named Eukaryotic translation initiation factor 2-alpha kinase, and its binding affinity score is −11.812 Kcal/mol. Figure 13 show 2D binding mode of compound 35 with Eukaryotic translation initiation factor 2-alpha kinase-3. Compounds 65 and 72 revealed no potential to target kinases through the ligTMap tool.  Compound 31 showed binding affinity equal to −8.53 against Cyclin-dependent kinase 2 and compound 35 showed more binding affinity towards protein kinase named Eukaryotic translation initiation factor 2-alpha kinase, and its binding affinity score is −11.812 Kcal/mol. Figure 13 show 2D binding mode of compound 35 with Eukaryotic translation initiation factor 2-alpha kinase-3. Compounds 65 and 72 revealed no potential to target kinases through the ligTMap tool. Compound 31 showed binding affinity equal to −8.53 against Cyclin-dependent kinase 2 and compound 35 showed more binding affinity towards protein kinase named Eukaryotic translation initiation factor 2-alpha kinase, and its binding affinity score is −11.812 Kcal/mol. Figure 13 show 2D binding mode of compound 35 with Eukaryotic translation initiation factor 2-alpha kinase-3. Compounds 65 and 72 revealed no potential to target kinases through the ligTMap tool.

Screening against Lipinski's and Veber's Rules
Lipinski and Veber descriptors were determined using Discovery studio 4.0. At first, the CHARMM force field was applied then the tested compounds were prepared and minimized according to the preparation of small molecule protocol. Then, the Lipinski and Veber descriptors protocol was applied to carry out these studies.

Computational ADMET Studies
ADMET descriptors (absorption, distribution, metabolism, excretion and toxicity) of the compounds were determined using Discovery studio 4.0. At first, the CHARMM force field was applied then the tested compounds were prepared and minimized according to the preparation of small molecule protocol. Then, ADMET descriptors protocol was applied to carry out these studies [42][43][44][45][46][47].

Computational Toxicity Studies
The toxicity parameters of the synthesized compounds were calculated using Discovery studio 4.0. IDK was used as a reference drug. At first, the CHARMM force field was applied then the compounds were prepared and minimized according to the preparation of small molecule protocol. Then, different parameters were calculated from the toxicity prediction (extensible) protocol [29,48].

Docking Studies
Crystallographic structure of tubulin [PDB ID: 5I3R, resolution: 2.40 Å] was retrieved from Protein Data Bank (http://www.pdb.org) and considered a target for docking simulation. The docking analysis was performed using MOE2014 software (Montreal, QC, Canada) to evaluate the free energy and binding modes of the synthesized compounds against BMP-2-inducible kinase. At first, the crystal structure of the target was prepared by removing water molecules and retaining the essential chain and the co-crystallized ligand (IDK). Then, the protein structure was protonated, and the hydrogen atoms were hidden. Next, the energy was minimized, and the binding pocket of the protein was defined.
The 2D structures of the tested compounds and reference ligand (IDK) were sketched using ChemBioDraw Ultra 14.0 and saved in MDL-SD format. Then, the saved files were opened using MOE and 3D structures were protonated. Next, energy minimization was applied. Before the docking process, validation of the docking protocol was carried out by running the simulation only using the co-crystallized ligand (IDK), which showed a low RMSD value. The molecular docking of the tested compounds was performed using a default protocol against the target receptor. In each case, 30 docked structures were generated using genetic algorithm searches, ASE was used for scoring, and forcefield (MMFF94X) for refinement. The ASE scoring function estimates the free energy of binding of the ligand from a given pose. The functional form is a sum of terms: The output from MOE was further analyzed and visualized using Discovery Studio 4.0 software [49][50][51][52][53][54][55].