Skip to main content

ORIGINAL RESEARCH article

Front. Pharmacol., 21 April 2022
Sec. Experimental Pharmacology and Drug Discovery
This article is part of the Research Topic Chemoinformatics Approaches to Structure- and Ligand-Based Drug Design, Volume II View all 20 articles

Identification of Potent and Selective JAK1 Lead Compounds Through Ligand-Based Drug Design Approaches

  • 1Computational Biology Lab, Department of Genetic Engineering, School of Bioengineering, SRM Institute of Science and Technology, SRM Nagar, Kattankulathur, India
  • 2Department of Clinical Immunology, Jawaharlal Institute of Post-Graduate Medical Education and Research, Pondicherry, India
  • 3Department of Chemistry and Department of Carbon Materials, Chosun University, Gwangju, South Korea

JAK1 plays a significant role in the intracellular signaling by interacting with cytokine receptors in different types of cells and is linked to the pathogenesis of various cancers and in the pathology of the immune system. In this study, ligand-based pharmacophore modeling combined with virtual screening and molecular docking methods was incorporated to identify the potent and selective lead compounds for JAK1. Initially, the ligand-based pharmacophore models were generated using a set of 52 JAK1 inhibitors named C-2 methyl/hydroxyethyl imidazopyrrolopyridines derivatives. Twenty-seven pharmacophore models with five and six pharmacophore features were generated and validated using potency and selectivity validation methods. During potency validation, the Guner-Henry score was calculated to check the accuracy of the generated models, whereas in selectivity validation, the pharmacophore models that are capable of identifying selective JAK1 inhibitors were evaluated. Based on the validation results, the best pharmacophore models ADHRRR, DDHRRR, DDRRR, DPRRR, DHRRR, ADRRR, DDHRR, and ADPRR were selected and taken for virtual screening against the Maybridge, Asinex, Chemdiv, Enamine, Lifechemicals, and Zinc database to identify the new molecules with novel scaffold that can bind to JAK1. A total of 4,265 hits were identified from screening and checked for acceptable drug-like properties. A total of 2,856 hits were selected after ADME predictions and taken for Glide molecular docking to assess the accurate binding modes of the lead candidates. Ninety molecules were shortlisted based on binding energy and H-bond interactions with the important residues of JAK1. The docking results were authenticated by calculating binding free energy for protein–ligand complexes using the MM-GBSA calculation and induced fit docking methods. Subsequently, the cross-docking approach was carried out to recognize the selective JAK1 lead compounds. Finally, top five lead compounds that were potent and selective against JAK1 were selected and validated using molecular dynamics simulation. Besides, the density functional theory study was also carried out for the selected leads. Through various computational studies, we observed good potency and selectivity of these lead compounds when compared with the drug ruxolitinib. Compounds such as T5923555 and T5923531 were found to be the best and can be further validated using in vitro and in vivo methods.

Introduction

Janus kinase 1 (JAK1) is the most widely employed JAK, according to biochemical and genetic research, since it is involved in the signaling of the gamma common (γc), beta common (βc), gp130, type I and type II interferon, IL-6, and IL-10 cytokine subfamilies (Kulagowski et al., 2012). JAK1 comprises seven homology domains (JH1–JH7) (Harpur et al., 1992). The C-terminal kinase module (JH1) is the protein’s physiologically active catalytic domain. The JH2 domain is a catalytically inactive pseudokinase domain that has been found to interact with the JH1 domain and control its activity (Saharinen and Silvennoinen 2002). Two Src homology 2 (SH2) domains (JH3 and JH4) precede the FERM domain (JH5–JH7) at the N-terminus. The JH1 domain has an ATP-binding site, which has been the target of a number of small-molecule inhibitors. All four members of JAK have a highly conserved kinase domain, particularly at the ATP-binding region, which complicates the development of particular inhibitors (Caspers et al., 2016). The active sites of JAKs comprise multiple subdomains that include the β-glycine loop, the catalytic loop, and activation loops (Taldaev et al., 2022). The amino acid present in and around the hinge region serves critical functions in the integrity of kinase activity control. Furthermore, since this area is adjacent to the ATP-binding site in the catalytic cleft, it is reasonable to believe that the mutations in this region might promote constitutive activation of the kinase (Gorden et al., 2010; Haan et al., 2010).

All STAT proteins (STAT1–STAT6) that are ubiquitously expressed in all the tissues may be phosphorylated by JAK1 enzyme (Gruber et al., 2020). JAK1 has been shown in mouse knockout experiments to have a critical function in signal transduction (Itteboina et al., 2017). According to earlier research, JAK1 is ascendant over JAK3, and in the absence of JAK1, JAK3 is unable to activate STATs (Haan et al., 2011). Furthermore, recent studies have shown that JAK1 rather than JAK3 kinase is the primary driver of the immune-relevant cytokine activity (Menet et al., 2015). JAK1 is involved in various types of cancer. Activation of JAK1 kinase by IL-6 family cytokines appeared to be the mechanism for constitutive STAT3 activation in human ovarian cancer cells (Wen et al., 2014). In gastric cancer, by activating the JAK1/STAT3 pathway, the upregulation of HOXA10 gene increased cell proliferation, cloning formation, and tumorigenesis and lowered cell apoptosis (Chen et al., 2019). In lung adenocarcinoma patients, the overall survival time was substantially reduced in patients with EGFR-amplified tumors expressing greater levels of phosphorylated JAK1 compared with individuals with tumors without one or both of these traits. Additionally, JAK inhibition was demonstrated to limit the development of human lung adenocarcinoma with a K-RAS mutation (Xie et al., 2021). AML and breast cancer patients have exhibited several STAT5-activating JAK1 mutations (Hornakova et al., 2011). Moreover, in ER-negative breast cancer cell lines, the upregulation of phosphorylated JAK1 expression was observed (Yeh et al., 2007).

According to clinical and experimental investigations, rheumatoid arthritis synovial response may be influenced by the JAK1-mediated cytokine (IFN and IL-6) signaling. As a result, inhibiting JAK1 is regarded as a significant therapeutic strategy for the successful treatment of rheumatoid arthritis (Keretsu et al., 2021b). Recently, it has been discovered that inhibiting JAK1 selectively may be an effective therapy option for patients suffering from autoimmune and hematological illnesses because of the role that altered JAK1 signaling plays in these conditions (Kleppe et al., 2017). Moreover, JAK1 expression in cancer cells allows individual cells to contract, perhaps enabling them to transcend their tumor and spread to other areas of the body (Nordqvist 2011). Mutations in JAK1 are less common than in T-ALL patients with B-ALL or leukemia of the myeloid origin. In two AML patients, a JAK1 mutation V623A was found, emphasizing the capacity of constitutively active JAK1 to induce a variety of leukemias (Xiang et al., 2008; Raivola et al., 2021).

JAK inhibitors, which have been authorized for the treatment of cancer and autoimmune illnesses, have provided the first insight on the importance of JAK1 in NK cell biology (Schwartz et al., 2017). Ruxolitinib, JAK1/JAK2 inhibitor, has lowered the number of NK cells and hampered maturation and function in both mice and human patients (Schönberg et al., 2015; Bottos et al., 2016). Ruxolitinib’s influence on NK cell development has been linked to JAK2 as well; therefore, it is not clear which of the two kinases is accountable for the reported results (Bottos et al., 2016; Kim et al., 2017). The fascinating finding by Sohn et al. (2013) has highlighted the importance of JAK1 inhibitor on the IL-6, IL-22, and INF-pathways. JAK1 inhibitors including ruxolitinib, tofacitinib, filgotinib, peficitinib, and numerous additional second-generation inhibitors are now under investigation for the treatment of inflammatory and autoimmune illnesses. Because of limited potency, non-targeting, and off-target effects (Keretsu et al., 2021a), new JAK1 inhibitors with high potency and selectivity are urgently needed.

Pharmacophore models are widely employed to quantitatively explore common chemical characteristics among a considerable number of structures with great diversity (Taha et al., 2008; Xie et al., 2009). It is one of the widely used approaches to search for chemical databases and identify novel scaffolds for various targets (Wang H. et al., 2008; Lu et al., 2007). To discover the potent hits, the ligand-based and structure-based pharmacophore models can be used. In this study, the ligand-based pharmacophore models were generated using the 52 JAK1 inhibitors reported by Zak et al. It elucidates the spatial arrangement of structural features of various potent and structurally diverse inhibitors crucial for biological recognition. One efficacious approach toward the discovery and development of the drugs is the virtual screening of molecular libraries (Stahl et al., 2006). Virtual screening helps to identify the potential lead molecules and reduces the time and cost of the drug discovery process (Reddy et al., 2007). Thus, pharmacophore-based virtual screening was implemented. In many research works, it was proposed that the combination of pharmacophore modeling and molecular docking is a successful method to discover the novel and potent lead compounds (Sakkiah et al., 2009; Sakkiah et al., 2010; Sakkiah et al., 2011). Hence, the results of pharmacophore-based virtual screening were taken for molecular docking.

Docking results were used to predict the binding orientations of the hits as well as the filter to select the hits. The molecular docking results were validated by calculating the free energy of binding using the molecular mechanics-generalized born surface area (MM-GBSA) method for the protein–ligand complexes (Friesner et al., 2006). Furthermore, induced fit docking (IFD) was carried out to get additional understanding about the structure and flexibility of these hits into the binding site since IFD has been reported to be a powerful method to account for both receptor and ligand flexibility (Zhong et al., 2009). Subsequently, the cross-docking method was used to identify the selective hits by docking every hit to every receptor. By examining the results, the top five hits were selected and taken for molecular dynamic simulation and density functional theory (DFT) study. To identify the potency and selectivity of the leads, a drug molecule named ruxolitinib was included in the study. The results of selected lead compounds and the drug were compared and analyzed.

Materials and Methods

Dataset Selection

For ligand-based pharmacophore modeling, a set of 52 JAK1 inhibitors (C-2 methyl/hydroxyethyl imidazopyrrolopyridines derivatives) reported by Zak et al. (2012) and Zak et al. (2013) were selected because of their diverse biological activity. The Ki values of these inhibitors (0.1–150 nM) were derived using biochemical and cell-based assays. These inhibitors have shown higher selectivity toward JAK1 over JAK2. The experimental Ki values were converted into pKi values that are simply the negative log of the Ki value. The chemical structures and biological activities of all molecules are given in Table 1.

TABLE 1
www.frontiersin.org

TABLE 1. The chemical structures and the biological activity of JAK1 inhibitors.

Pharmacophore Model Generation

Phase 4.3, a high-performance program module of Schrödinger 2015, was used to generate the ligand-based pharmacophore models (Dixon et al., 2006). It uses a fine-grained conformational sampling method to predict a hypothesis consisting of the common pharmacophore features. The Ligprep module was used to clean, minimize, and generate conformations of all compounds. Based on the diversity of the chemical structure and its biological activity, the quantitative pharmacophore models were generated using the Develop Pharmacophore Model option. On the basis of biological activity distribution (pKi values), the activity threshold value was set and the inhibitors were divided into actives, inactives, and moderately actives. In this study, both five and six featured pharmacophore hypotheses were generated by defining the minimum and maximum numbers of sites to five and six. The pharmacophore models were developed possessing different combinations of hydrogen-bond acceptor (A), hydrogen-bond donor (D), aromatic ring (R), hydrophobic group (H), positively ionizable (P), and negatively ionizable (N) groups. The resulting hypotheses were scored and ranked on the basis of scoring parameters. The scoring algorithm includes the alignment of site points and vectors, number of ligands matched, volume overlap, relative conformational energy, selectivity, and activity. The difference between the survival score and the survival inactive score notifies the ability of the hypotheses to correctly distinguish between actives and inactives.

Pharmacophore Model Validation

Since the pharmacophore model is just a theoretical model, it is necessary to analyze whether or not the generated model is able to predict the active compounds. Thus, two approaches, namely, potency validation and selectivity validation, were performed to measure the accuracy of pharmacophores in selecting the active compounds.

Potency Validation

Potency validation was carried out to test whether the pharmacophore model is good enough to pick a greater number of active molecules. This was achieved by screening the database consisting of both active molecules and decoys. Active molecules are the known inhibitors of JAK with higher biological activities, whereas decoys are the molecule that does not have any activity toward JAK and it was downloaded from DUD-E (a Database of Useful Decoys-Enhanced) database (Mysinger et al., 2012). DUD-E datasets were used only after removing the biasness through docking. Based on the number of actives and decoys retrieved by the pharmacophore models, statistical parameters such as Guner-Henry (GH) score, %A, %Y, and E score were calculated using the following formula:

GH score= (Ha (3A + Ht)4HtA)(1 Ht  HaD  A);
%A=HaA100;%Y=HaHt100;E=Ha/HtA/D,

where Ha is the number of actives in the hits list, Ht is the number of hits retrieved, A is the number of active compounds in the database, D is the number of compounds in the database, %A is the percentage of known active compounds obtained from the database, %Y is the percentage of known actives in the hits list, and E is the enrichment of the concentration of actives by the model relative to random screening without a pharmacophoric approach. GH score ranges from 0 to 1, which indicates a null model and an ideal model, respectively. GH score >0.6 indicates the acceptable quality of the pharmacophore model and is useful in differentiating the known active molecules from inactives and suitable for retrieving active JAK1 inhibitors (Sathe et al., 2014; Li et al., 2015).

Selectivity Validation

Selectivity validation was performed to check which pharmacophore models are more selective in choosing high number of JAK1 molecules. Selectivity validation was carried out in two ways. First, a database comprising 30 JAK1, 30 JAK2, 30 JAK3, and 30 TYK2 molecules (Yang et al., 2007; Wang et al., 2009; Pissot-Soldermann et al., 2010; Ioannidis et al., 2011; Kulagowski et al., 2012) was created and used for validation. The ability of pharmacophore models to differentiate the selective JAK inhibitors was evaluated using virtual screening workflow on a manually curated database. Second, to further confirm the selectivity of the selected models, the available 288 JAK1 (Kulagowski et al., 2012; Labadie et al., 2012; Hurley et al., 2013; Labadie et al., 2013), 627 JAK2 (Lucet et al., 2006; Wang et al., 2009; Pissot-Soldermann et al., 2010; Harikrishnan et al., 2011; Ioannidis et al., 2011; Schenkel et al., 2011; Dugan et al., 2012; Forsyth et al., 2012; Lynch et al., 2013; Vazquez et al., 2018), and 431 JAK3 (Chrencik et al., 2010; Thoma et al., 2011; Jaime-Figueroa et al., 2013; Lynch et al., 2013; Soth et al., 2013; De Vicente et al., 2014; Duan et al., 2014) inhibitors from diverse research papers that mention either IC50 values or Ki values of these inhibitors were taken for validation.

Pharmacophore-Based Virtual Screening

Virtual screening is the process where the complete databases are used to identify the molecules in the database which are most likely to bind to a drug target (Vyas et al., 2008). In this study, pharmacophore-based virtual screening was carried out using the “find matches to hypothesis” option available in the phase module which efficiently search for pharmacophore matches from the database of fixed conformers. The pharmacophore-based virtual screening was performed against Maybridge (53,000) (www.maybridge.com), Lifechemicals (12, 92, 000) (https://lifechemicals.com/), Enamine (24,91,318) (https://enamine.net/), Chemdiv (15,00,000) (https://www.chemdiv.com/), Asinex (398,022) (https://www.asinex.com/), and Zinc chemical and Zinc natural databases (https://zinc.docking.org/) (44,92,226) (Irwin and Shoichet 2005; Irwin et al., 2012; Sterling and Irwin 2015) to identify the new molecules with novel scaffolds. After screening, fitness score that is a measure of how well the hypothesis matched to the aligned ligand conformers based on RMSD site matching, volume terms, and vector alignments was used to filter the molecules.

Absorption, Distribution, Metabolism, and Excretion Prediction

After virtual screening, the molecular descriptors and pharmaceutically applicable properties of the hits were calculated using Qikprop 4.4. Qikprop generates the physicochemical properties for a compound to find whether the compound follows drug likeliness properties. Lipinski’s rule characterizes the important molecular properties of drug, including absorption, distribution, metabolism, and excretion (ADME) that is essential for a drug’s pharmacokinetics in the human body (Jorgensen and Duffy 2002). Parameters that determine the ADME of the molecules were Molweight (Molecular weight), QPlogPo/w (partition coefficient), QPlogS (water solubility), percentage of human oral absorption, and intestinal absorption parameters such as Caco-2 and MDCK permeability. The compounds are expected to be active in humans only if the molecule passes through Lipinski’s rule of five. Therefore, the compounds retrieved after filtration were subjected to ADME prediction and its physicochemical properties were analyzed.

Molecular Docking

Molecular docking predicts the binding mode and interaction of the small molecule to the protein. It distinguishes the behavior of small molecules in the binding site of target protein and explicates its fundamental biochemical processes (Gschwend et al., 1996; Lipinski 2000). The binding conformations of the hits inside the JAK1 ATP-binding site were investigated using Grid-based Ligand Docking with Energetics (Glide 6.7) module. The ATP-binding site of JAK1 comprises Leu881, Glu883, Val889, Ala906, Met956, Glu957, Phe958, Leu959, Gly962, Ser963, Glu966, Arg1007, Asn1008, Leu1010, Gly1020, and Asp1021 residues. Before docking, protein preparation wizard was used to prepare protein structure (3EYG) (Williams et al., 2009) applying the default parameters that include adding hydrogens, filling missing atoms and residues using PRIME, assigning correct bond orders, and hydrogen-bond optimization and minimization. In the Receptor Grid Generation panel, the center of the gird box was defined on the centroid of the co-crystallized ligand (MI1), and the volume in the active-site region of the receptor was calculated by default settings (van der Waals radius scaling factor 1.0 and partial charge cutoff 0.25). Molecular docking was performed using both the Standard Precision (SP) and Extra Precision (XP) docking modes in which the receptor was held rigid and the ligand was free to move (Jain 2003; Halgren et al., 2004). Glide score is a combination of hydrophobic, hydrophilic, van der Waals energy, metal binding groups, freezing rotatable bonds, and polar interactions with the receptor. Comparing Glide SP and XP score, Glide SP score is a softer and more forgiving function whereas Glide XP score is a harder function and adept at reducing the false positives. Therefore, Glide XP score was considered for the selection of hits and further analysis.

MM-GBSA Calculations

The binding free-energy calculations procured via the MM-GBSA method are more precise and consistent than the glide XP score and improve the ranking of potential leads (Lyne et al., 2006; Das et al., 2009; Yang et al., 2009). Therefore, the binding free energy (∆G bind) of the protein–ligand complexes was calculated using the Prime MM-GBSA module implemented in Schrödinger 2015. The Prime MM-GBSA Module incorporates the OPLS3 force field and the VSGB dissolvable model to look through calculations (Li et al., 2011). The energy difference between the free and complex states of protein and ligand was calculated. The energy components such as covalent binding energy, van der Waals energy, generalized born electrostatic solvation energy, Coulomb energy, total energy, and H-bond correction were retrieved from the calculations.

Induced Fit Docking

In the docking protocol, to retain the flexibility of the receptor, a mixed molecular docking protocol called induced fit docking (IFD) developed by Schrödinger 2015 was employed (Wang H. Y. et al., 2008). IFD uses the refinement module in Prime to account for the receptor flexibility and Glide to account for the ligand flexibility (Jacobson et al., 2004). Protein preparation wizard and the Ligprep module were used for protein and ligand preparation, respectively. Grid was generated on the ATP-binding site amino acid residues based on the co-crystallized ligand. The ATP-binding site residues and their flexibility were considered for the IFD protocol. IFD was carried out with default parameters, and 20 conformational poses were calculated for each ligand. IFD scores (IFD score = 1.0 Glide_Gscore +0.05 Prime_Energy) were calculated based upon the total energy of the system and the protein–ligand interaction energy and used to rank the IFD poses (Luo et al., 2014). The electrostatic interactions formed between the receptor and the ligand were calculated by the docking scores under “Electro,” and hydrophilic interactions under “Lipophilic Evdw” mention the lipophilicity component acquired from the hydrophobic grid. IFD poses were ranked based on the scores, and the best pose was chosen for each hit.

Cross Docking

Cross docking is the process of taking a series of complexes of ligand–receptor pairs and docking every ligand to every receptor. This is used to study the specificity of the ligands and the receptors and, thus, yield valuable report regarding the effects of ligand upon binding. Protein preparation wizard and LigPrep were used to prepare proteins and the shortlisted hits, respectively. Grid was generated on the ATP-binding site residues. The hits shortlisted from the molecular docking study were docked against JAK1, JAK2, and JAK3 using the Glide XP module to identify the selective lead compounds.

Molecular Dynamics Simulation

Docking results could be the instantaneous state and were not considered decisive because binding of the inhibitor to a protein in an in vivo state is a dynamic process. For advanced studies, the stable binding mode of the ligand is more reliable. Hence, to explore the detailed binding modes and compare the stability and molecular interactions of the docked lead complexes, molecular dynamics simulation was carried out for 100ns using GROningen MAchine for Chemical Simulations (GROMACS version 2016.3 installed in Centos 7.3) software (Abraham et al., 2015). GROMACS works according to Newton’s laws of motions and simulates the behavior of bio-molecules such as nucleic acids, proteins, lipids, ligands, ions, and water. The coordinates for MD simulations have been achieved from the docking results. The PRODRG server (http://davapc1.bioch.dundee.ac.uk/cgi-bin/prodrg) was used to calculate the ligand parameters in the framework of GROMOS96 54a7 force field. The SPC water model was used as a solvent during simulation. To achieve the stability of the simulated system, the potential energy, temperature, and pressure were monitored during the simulations. The temperature and pressure of the system were equilibrated (from ps to ns) till they reach 300 K and 1.05 bar, respectively. The stability of the secondary structure elements and conformational changes of the simulated complexes were evaluated by root mean square deviations (RMSDs), root mean square fluctuation (RMSF), radius of gyration (Rg), and solvent-accessible surface area (SASA) values obtained from MD trajectories. The molecular dynamics study was performed using High-Performance Computing server (Intel Xeon 14 core processor with 28 threads and 2.40 GHz processor speed).

MM-PBSA Calculation

The molecular mechanics energies combined with the Poisson–Boltzmann and surface area continuum solvation (MM/PBSA) method have been applied to predict binding free energies and to evaluate the relative stabilities of different bimolecular structures. The MM/PBSA calculations were performed for the simulated systems using g_mmpbsa, a GROMACS Tool for High-Throughput MM-PBSA Calculations (Kumari et al., 2014). Combined with molecular dynamics (MD) simulations, MM-PBSA can also incorporate conformational fluctuations and entropic contributions to the binding energy (Homeyer and Gohlke 2012).

Density Functional Theory Study

The density functional theory (DFT) study was carried out to observe the chemical behavior of the lead compounds using the electron density-relevant concepts (Zhao et al., 2011). Also, it provides a quantum-level understanding of the molecules and assists in building the relationship between the electronic properties and the biological activity of the molecule (Nagarajan et al., 2018). Molecular descriptors such as total energy, highest occupied molecular orbital (HOMO), lowest unoccupied molecular orbital (LUMO), band energy gap (ΔE), molecular dipole moment, absolute hardness (η), global softness (σ), chemical potential (μ), electronegativity (χ), and electrophilicity index (ω) were studied for the selected lead compounds using Gaussian 16 software. Initially, the molecules were optimized using the B3LYP function with a 6-31G(d) basis set to calculate their molecular properties such as total energy and molecule dipole moment (Becke 1998). The dipole moment relates to the electro-chemical reactivity of the compounds. The electron donating and accepting ability of the molecules HOMO energy (EHOMO) and LUMO energy (ELUMO), respectively, were calculated.

Results and Discussion

Pharmacophore Model Generation

In the phase module, ligand-based pharmacophore model generation was carried out utilizing 52 JAK1 inhibitors named C-2 methyl/hydroxyethyl imidazopyrrolopyridine derivatives (Table 1) along with their activity values. Ten molecules whose pKi > 8.9 were taken as actives, twelve molecules whose pKi < 8.1 were taken as inactives, and the remaining thirty-one molecules were considered to be intermediates. Twenty-seven different pharmacophore hypotheses (six with six featured pharmacophores and twenty-one with five featured pharmacophores) were generated and put through the stringent scoring function. The generated pharmacophores were ranked by aligning them with the active ligands, and the statistical data obtained after scoring are tabulated in Table 2. Besides the survival active score, survival inactive score, and post-hoc score, fitness score was considered to measure the quality of the pharmacophores. The fitness score was calculated between the pharmacophores and the highly active (compound 51) and highly inactive (compound 7) compounds in the dataset. For all pharmacophores, the fitness score was higher with the highly active compound compared with the inactive compound. Subsequently, the pharmacophores were evaluated using different validation methods.

TABLE 2
www.frontiersin.org

TABLE 2. The summary of statistical data obtained for the pharmacophore hypotheses.

Pharmacophore Model Validation

Potency Validation

For potency validation, a database containing 40 JAK1 actives and 1,000 decoys was created. The generated pharmacophore models were allowed to screen this database to calculate the GH score. It was observed that six featured hypotheses have picked less number of decoys compared with five featured hypotheses. The ADHRRR, DDHRRR, and DPRRR hypotheses were more potent because they do not pick any decoys. DDRRR, ADPRR, and AHRRR have picked very less number of decoys. DDRRR, ADHRR, DHRRR, AADHR, DDHRR, and DDRRR have picked more active molecules. The results of potency validation are tabulated in Table 3. Based on the number of actives and decoys retrieved by the hypotheses, the GH score was calculated. The hypotheses such as DDHRRR, ADHRRR, ADDHRR, DDRRR, DPRRR, DHRRR, ADRRR, AHRRR, DDHRR, and ADPRR have obtained the GH score >0.6 indicating the goodness of these hypotheses.

TABLE 3
www.frontiersin.org

TABLE 3. Pharmacophore validation results from potency validation.

Selectivity Validation

Initially, the selectivity validation was performed with a set of 30 JAK1, 30 JAK2, 30 JAK3, and 30 TYK2 molecules retrieved from different studies. DPRRR has picked only JAK1 molecules. The DHRRR, ADHRR, DDRRR, and DDHRR hypotheses have picked a high number of JAK1 molecules and very less JAK3 molecules. ADHRRR and DDHRRR have picked only few JAK1 and one JAK3 molecules. The results of selectivity validation are tabulated in Table 4. The fitness score for JAK1 inhibitors was greater than or equal to 1.5, whereas for other JAK inhibitors, the fitness score was <1.5 for most of the molecules indicating that the pharmacophore models were able to map well with the JAK1 inhibitors (Sathe et al., 2014; Babu et al., 2015).

TABLE 4
www.frontiersin.org

TABLE 4. Pharmacophore validation results from selectivity validation.

Six feature pharmacophore hypotheses were more potent but not highly selective to JAK1. Based on potency and selectivity validation results, the DDHRRR, ADHRRR, DDRRR, DPRRR, DHRRR, ADRRR, DDHRR, and ADPRR hypotheses were selected because they were successful in retrieving active compounds from the database. The representation of the selected JAK1 pharmacophore models showing the distances between the pharmacophoric sites is shown in Figures 1A–H. On mapping the selected pharmacophore models with highly active compound 51 and inactive compound 7, it was observed that the fitness score was >2.5 for the highly active compound mapping with all pharmacophore features whereas inactive compound 7 could map with either four or five pharmacophore features with low fitness score. The highest fitness score with compound 51 suggests screening using these models would pick the similar active compounds. From the results, we suggest the combination of two or three aromatic rings (R) and one or two donor atoms (D) with a hydrophobic (H) group is an important pharmacophoric feature for identifying the selective JAK1 inhibitors. The important pharmacophore features obtained were compared with the contribution maps obtained through the hologram-based fingerprint technique (Supplementary Figure S1). The contribution maps depict the imidazopyrrolopyridine ring which possesses one donor and three aromatic rings responsible for the intermediate contribution of the inhibitory activity. From the highly active compounds, we observed the cyano group attached to cyclohexanes (yellow) and the hydroxyethyl group attached to imidazopyrrolopyridines (green); a hydrophobic and an aromatic/donor group, respectively, are strongly responsible for the higher activity. Thus, these pharmacophore features are highly important for the inhibitory activity of JAK1.

FIGURE 1
www.frontiersin.org

FIGURE 1. The representation of selected pharmacophore models (A) ADHRRR, (B) DDHRRR, (C) DDRRR, (D) DPRRR, (E) DHRRR, (F) ADRRR, (G) DDHRR, and (H) ADPRR. Pharmacophore features are colored in light blue, brown, dark blue, brick red, and green contours representing the H-bond donor (D), H-bond acceptor (A), positives (P), aromatic ring (R), and hydrophobic (H) groups, respectively. The distances between the pharmacophore features (A˚) are given in pink dotted lines.

To confirm the selectivity of the selected pharmacophore models (DDHRRR, ADHRRR, DDRRR, DPRRR, DHRRR, ADRRR, DDHRR, and ADPRR), the second round of selectivity validation was carried out with a set of 288 JAK1, 627 JAK2, and 431 JAK3 inhibitors with diverse activity. It was observed (Table 4) that all the selected pharmacophore models were able to pick more number of JAK1 inhibitors compared with its subtypes. Hence, the selected pharmacophore models were capable of discriminating the JAK1 inhibitors and appropriate for retrieving the novel and selective JAK1 inhibitors.

Pharmacophore-Based Virtual Screening

The selected pharmacophore models were screened against Maybridge, Lifechemicals, Enamine, Chemdiv, Asinex, and Zinc (chemical and natural) databases for the identification of new hits. The identified hits contain the structural features that overlap with the selected pharmacophore models. The hits obtained were ranked and filtered based on the fitness score. The fitness score was set to >1.5 for the Maybridge, Asinex, Chemdiv, Lifechemicals, and Enamine databases, whereas for the Zinc database, the fitness score was set to >2 because of high number of molecules retrieved from the Zinc database. As a result of screening and filtration, 4,265 compounds were retrieved. The total numbers of hits retrieved from different databases are tabulated in Table 5.

TABLE 5
www.frontiersin.org

TABLE 5. Number of hits obtained from pharmacophore-based virtual screening.

The potentiality of the pharmacophore models was validated using receiver operating curves (ROCs) utilizing the screened molecules (Hevener et al., 2009); 10 compounds identified from the pharmacophore-based virtual screening were seeded with 500 decoys. Enrichment was estimated based on how well the compounds were fetched. After ranking the decoy set and docked compounds by the Glide score, the enrichment was calculated using the ROC plot that provides the report on sensitivity and specificity. The ROC plot inferred that Glide XP ranked seven compounds in top 10% with the ROC value as 0.93 and the AUC value as 0.92. 80% of the true positives were fetched in top 20% which indicates its capability of retrieving the active compounds. The gentle increase in the ROC curve (Supplementary Figure S2) was noticed in the beginning, which implies that number of true positives was sacrificed to reduce the amount of false positives.

Absorption, Distribution, Metabolism, and Excretion Prediction

Compounds that pass Lipinski’s rule of five and other ADME properties of the drug are expected to be active in humans. Properties such as molecular weight, H-bond donors, H-bond acceptors, log p, van der Waals surface, aqueous solubility, blockage of HERG K+ channels, apparent Caco-2 cell permeability, apparent MDCK cell permeability, brain/blood partition coefficient, skin permeability, binding to human serum albumin, and human oral absorption of the hits were studied. Finally, 2,856 compounds whose drug-like properties were in the acceptable range (according to qikprop recommended range) were selected and subsequently exposed to glide SP and XP docking protocols to remove both the false-positive and false-negative hits.

Molecular Docking

The molecular docking study was carried out using the Glide SP and XP modes to explore the binding mode and interaction of hits on the ATP-binding site. The crystal structure of JAK1 protein 3EYG in complex with MI1 (Williams et al., 2009) was used to perform molecular docking. The grid was developed on the centroid of co-crystallized ligand MI1 surrounding the ATP-binding site residues (Leu881, Glu883, Val889, Ala906, Met956, Glu957, Phe958, Leu959, Gly962, Ser963, Glu966, Arg1007, Asn1008, Leu1010, Gly1020, and Asp1021) of JAK1. Initially, the docking of MI1 into the ATP-binding site was performed to check the accuracy and reproducibility of the docking program. Subsequently, the highly active compound 51 and 2,856 hits were docked into the ATP-binding site. Considering the docking result of compound 51 (glide XP score -9.691), the glide XP threshold value was set to ≥ −9.60 to identify the novel hits. We observed that 90 molecules have exhibited glide score greater than the threshold and it was shortlisted (Supplementary Table S1). Among the JAK1 ATP-binding site residues, Leu959 and Glu957 that are present in the hinge region were found to be the most selective amino acid residues for the H-bond interaction and also crucial for selective inhibition of JAK1. Hence, the interactions with Leu959 and Glu957 were investigated for the hits. Compound 51 has shown H-bond interactions with Leu959, Glu957, and Leu881. The selected 90 hits have exhibited H-bond interaction with either Leu959 or Glu957 or both residues. Additionally, the Leu881, Glu883, Ser963, Glu966, and Arg1007 residues were involved in H-bond interaction with most of the hits. The hydrophobic interactions were formed mainly by the residues Leu881, Val889, Ala906, Val938, Met956, Phe958, Pro960, and Leu1010. The binding of compound 51 into the ATP-binding site is shown in Figure 2.

FIGURE 2
www.frontiersin.org

FIGURE 2. The binding of highly active compound 51 into the ATP-binding site of JAK1.

MM-GBSA Calculations

The highly ranked hits selected from glide docking were taken for MM-GBSA calculations to predict the binding energy of the protein–ligand complexes. The calculated free energy of binding (∆G bind) was lower than glide energy. It was observed that van der Waals (∆G_Bind_vdW) energy contributes more for the ligand binding, whereas covalent interaction (∆G_Bind_Covalent) and electrostatic salvation (∆G_Bind_Solv_GB) energy terms disfavor for the inhibitor binding.

Induced Fit Docking

IFD was performed on the highly ranked hits using the crystal structure of JAK1 (3EYG). It was observed that IFD also produces good IFD score and XP score comparable to glide XP score. The IFD score of JAK1 hits was greater than or equal to −590, and their corresponding XP score was greater than −8.00 which indicates the good binding ability of the hits. The observed H-bond and hydrophobic interaction with the IFD results was highly similar to glide results, indicating that these hits could bind and produce similar H-bond and hydrophobic interaction inside the binding site upon both receptor and ligand flexibility.

Cross Docking

Since an important objective of this work is to attain admissible levels of intra-family selectivity, the cross-docking approach was employed for the highly ranked hits. For cross docking, the crystal structure of JAK1 (3EYG), JAK2 (3KRR), and JAK3 (3ZEP) was used. Among 90 hits tested, the top five compounds (T6649932, ST088474, T5923555, T5923531, and T6763842) that have the highest docking score toward JAK1 (>-10.5) compared with JAK2 and JAK3 were selected and taken for further study.

Analysis of Selected Lead Compounds

The top five compounds that showed good potency and selectivity were selected and analyzed. To identify the potency and selectivity of the leads, a drug molecule named ruxolitinib was included in the study. Molecular docking, MM-GBSA calculations, IFD docking, and cross docking were performed for the drug and compared with the selected leads. Subsequently, the selected leads were taken for the molecular dynamics simulation study using the GROMACS and DFT calculations using Gaussian. The chemical structures of the selected lead compounds are shown in Figures 3A–E.

FIGURE 3
www.frontiersin.org

FIGURE 3. The chemical structure of selected lead compounds. (A) T6649932, (B) ST088474, (C) T5923555, (D) T5923531, and (E) T6763842.

Absorption, Distribution, Metabolism, and Excretion Properties

ADME properties are the key determinants for the successful development of new drugs. All the analyzed pharmacokinetic parameters of these lead compounds were found to be within the permissible range. The percentage of the human oral absorption was found to be greater than 50%. The partition coefficient and water solubility that are important for the assessment of absorption and distribution of drugs within the body ranged between −0.4 and 3.5 and −5.4 and −2.1, respectively. Compounds T6649932, T5923555, and T5923531 possessing good Caco-2 and MDCK permeability have good level of intestinal absorption. The drug-likeliness properties of the selected leads and ruxolitinib are given in Table 6.

TABLE 6
www.frontiersin.org

TABLE 6. The drug likeliness properties of the selected lead compounds and the drug.

Glide XP Docking Analysis

For the selected leads, the glide XP score was greater than −10.015, whereas for ruxolitinib, it was −9.282. The highest docking score was observed for T6763842 (−10.671). The major contribution of vdW interactions was observed which indicate that the vdW interaction favors the protein-ligand complex. The glide XP score, glide energy, glide evdw, and glide ecoul of the selected leads are given in Table 7. On analyzing the interaction, it was observed that the lead compounds showed conserved H-bond interactions with both the selective residues of JAK1 (Leu959 and Glu957) similar to the drug indicating its remarkable selectivity. Compounds such as T6649932, T5923555, and T5923531 formed another H-bond with Arg1007. Additionally, hydrophobic interactions were formed with the ATP-binding site residues Leu881, Val889, Ala906, Val938, Met956, Phe958, Pro960, and Leu1010. Figures 4A–F show the docked pose of lead compounds and drug inside JAK1 ATP-binding site.

TABLE 7
www.frontiersin.org

TABLE 7. Molecular docking results of the selected JAK1 lead compounds and the drug.

FIGURE 4
www.frontiersin.org

FIGURE 4. The representation of docked lead compounds and drug ((A) T6649932, (B) ST088474, (C) T5923555, (D) T5923531, (E) T6763842, and (F) ruxolitinib) present inside the ATP-binding site of JAK1 after molecular docking.

MM-GBSA Analysis

The binding free energy of the selected lead compounds ranges from −41.698 to −46.430 suggesting good binding affinity with JAK1. Many lead compounds have shown comparable free energy of binding with ruxolitinib. This provides the insight that these lead compounds have exhibited good specificity. Furthermore, the contribution of ∆G_Bind_vdW and ∆G_Bind Lipo components to the binding free energy was compared. The high binding free energy was majorly contributed by ∆G_Bind vdW than ∆G_Bind Lipo component. The predicted binding free energy of the selected leads and drug is tabulated in Table 8.

TABLE 8
www.frontiersin.org

TABLE 8. MM-GBSA results of the selected JAK1 lead compounds and the drug.

Induced Fit Docking Analysis

The accuracy of glide scoring function in identifying the leads was checked using the IFD method. The IFD score was greater than −594, and their corresponding docking score was greater than −8.5. The IFD scores of four lead compounds were higher compared to ruxolitinib (−595.395), whereas the docking scores of four lead compounds were little lower compared with ruxolitinib (−9.725). The highest IFD score and XP score were observed for T6763842 and T5923531. The Electro and Lipophilic Evdw scores of the lead compounds and drugs showed higher lipophilicity compared with electrostatic interactions which implies the important role of lipophilicity in inhibitory activity. The glide XP score, IFD score, lipophilic evdw, electro, and H-bond interaction of the selected leads are given in Table 9. The representation of docked lead compounds and drug present inside the ATP-binding site of JAK1 after induced fit docking is shown in Supplementary Figure S3. The IFD results also confirmed that the selected lead compounds have occupied the ATP-binding site of JAK1 irrespective of receptor flexibility.

TABLE 9
www.frontiersin.org

TABLE 9. Induced fit docking results of the selected JAK1 lead compounds and the drug.

Cross-Docking Analysis

The cross-docking results of selected leads indicate that their docking score was greater than −10.00 with JAK1, whereas with respect to other JAKs, their docking score ranges from −5.2 to −8.9. For ruxolitinib, docking scores were −9.178 with JAK1, −9.091 with JAK2, and −10.209 with JAK3. Therefore, the selected lead compounds have shown good selectivity in terms of docking score compared with ruxolitinib. The important components to determine the selectivity of the lead compounds were the electrostatic and hydrophilic components of the docking score. The drug and lead compounds showed higher lipophilicity compared with electrostatic interactions which implies that lipophilicity plays an important role in dictating the selectivity of these molecules. The cross-docking results of the drug and selected leads are given in Table 10. The cross-docking results of ruxolitinib showed higher binding affinity with JAK3 compared with JAK1 and JAK2 indicating its lesser selectivity, whereas the selected lead compounds showed greater affinity and selectivity toward JAK1 compared with other JAK subtypes. Hence, the cross-docking results indicate the selected lead compounds are more selective than the drug.

TABLE 10
www.frontiersin.org

TABLE 10. Cross-docking results of the selected JAK1 lead compounds and the drug.

Molecular Dynamics Simulation

RMSD Plot Analysis

RMSD relative to the respective initial conformations was monitored and analyzed to examine the stability and equilibration of all systems. The RMSD value for both the lead compounds and drug was calculated for the 100ns time scale using the apo form of JAK1 (3EYG) as reference. In Figure 5, it was observed that the RMSD values of both the drug and lead complexes were stable throughout the simulation. Furthermore, RMSD values for all protein backbone atoms attained convergence after 1 ns and maintained a plateau of 0.01 nm after the initial convergence. This suggests that all system simulations reached equilibrium and stabilization during the simulation. RMSD of both drug and lead complexes was found to be in the range of 0.02–0.03 nm indicating the similar stability. The observed smaller RMSD fluctuations for all compounds confirmed that the obtained binding conformations of these lead compounds and drug were highly reasonable.

FIGURE 5
www.frontiersin.org

FIGURE 5. The change in RMSD values of the backbone Cα atoms of JAK1 systems over a period of 100 ns after binding with the lead compounds and drug.

RMSF Plot Analysis

RMSF of backbone atoms were monitored to identify the strong binding interactions and exemplify the pliability of these lead complexes in the ATP-binding site of JAK1. The RMSF plot shown in Figure 6 indicates very minimal fluctuations were observed during the simulation except the terminal and loop regions of the protein. Most of the fluctuations were between 0.016 and 0.035 nm indicating the stability of the simulated system. Very minimal fluctuations were observed in the residues Pro912, His918, Glu946, Asn950, and Gly951 for all lead compounds and ruxolitinib. The JAK1 ATP-binding site residues that are crucial for binding and fixing the inhibitors have displayed insignificant fluctuations during the course of simulation. However, the most important and selective amino acid residues Leu959 and Glu957 that are important for inhibitor binding attained a quite stable behavior.

FIGURE 6
www.frontiersin.org

FIGURE 6. The change in RMSF values of JAK1 residues over a period of 100 ns after binding with the lead compounds and drug.

Rg Plot Analysis

The level of compactness in the structure of protein due to the presence or absence of ligands was calculated using radius of gyration (Rg) plot (Lobanov et al., 2008). It can be observed that all lead complexes and the drug showed consistently lower Rg values and exhibited a relatively similar nature of compactness in Figure 7. Thus, a relatively consistent Rg value indicates that a stably folded structure was observed throughout the MD simulation.

FIGURE 7
www.frontiersin.org

FIGURE 7. The change in Rg values over a period of 100 ns after binding with the lead compounds and drug.

Solvent Accessible Surface Area Plot Analysis

The solvent accessible surface area (SASA) calculation of the protein–ligand complexes was used for predicting the extent of the conformational changes that occurred during the interaction. The SASA plot shown in Figure 8 indicates that no significant changes in the protein structure were caused by these lead compounds and drug during simulation. Hence, the protein–ligand complexes are relatively stable throughout the simulation.

FIGURE 8
www.frontiersin.org

FIGURE 8. The change in SASA values over a period of 100 ns after binding with the lead compounds and drug.

Protein–Ligand Interaction Analysis

The most significant part in MD simulations is the analysis of protein–ligand interactions because it illustrates the changes in the binding mode of the ligands during simulations. Figure 9 shows the number of H-bond formations over the trajectory for lead compounds and the drug. The H-bonds were the principal binding forces between protein and ligand. The drug ruxolitinib has produced 2–4 H-bonds, whereas lead compounds have produced 0–2 H-bonds throughout the simulation. T5923555, T5923531, and T6763842 have produced 1–3 H-bonds, whereas ST088474 produced one H-bond with JAK1 all through the simulation. Ruxolitinib, T6763842, and T5923555 had retained two H-bonds, and T5923531 had retained one H-bond, whereas T6649932 does not have an H-bond at the end of the simulation. A strong hydrogen bond network was formed mainly by the residues Glu957 and Leu959. T5923555 and T5923531 retained hydrogen bonds with Leu959 and Glu957 at the end of simulation. Moreover, the ATP-binding site residues were almost hydrophobic, which can form strong nonpolar interactions with lead compounds. The detailed protein–ligand interaction residues before and after molecular dynamics simulation (Saddala and Adi 2018) were studied and are given in Table 11. T5923555 and T5923531 were found to be more stable and reliable before and after simulation, and their important interaction (Glu957 and Leu959) remains unchanged throughout the simulation.

FIGURE 9
www.frontiersin.org

FIGURE 9. The number of hydrogen bonds formed by lead compounds and drug over the simulation time.

TABLE 11
www.frontiersin.org

TABLE 11. The protein–ligand interaction analysis of the selected JAK1 lead compounds and the drug before, during, and after MD simulation.

The binding mode of the drug and lead compounds after simulation is represented in Figures 10A–F. It was inferred that the initial docked conformation and the final conformation of the lead compounds and drug lie in the same binding pocket (Supplementary Figure S4). Hence, the conformation of the lead compounds was stable inside the binding pocket which, in turn, validates the reliability of the docking results. Furthermore, these absolute results suggest that the identified lead compounds are highly selective and potent and they can be taken for in vitro and in vivo studies.

FIGURE 10
www.frontiersin.org

FIGURE 10. The representation of final conformation of the docked lead compounds and drug ((A) T6649932, (B) ST088474, (C) T5923555, (D) T5923531, (E) T6763842, and (F) ruxolitinib) present inside the ATP-binding site of JAK after molecular dynamics simulation.

MM-PBSA Calculation

The average binding energy of all the simulated complexes was calculated using the g_mmpbsa tool. The van der Waals energy, electrostatic energy, polar solvation energy, solvent-accessible surface area (SASA) energy, and binding energy were calculated and are tabulated in Table 12. T5923555 and T5923531 have shown good binding energy and van der Waals energy compared with other compounds.

TABLE 12
www.frontiersin.org

TABLE 12. MM-PBSA results obtained from the molecular dynamics trajectory for the selected JAK1 lead compounds and the drug.

Density Functional Theory Calculation

Molecular descriptors based on the electron density of the molecules were studied using Gaussian. Based on HOMO energy (EHOMO) and LUMO energy (ELUMO), descriptors such as ΔE, η, σ, μ, χ, and ω were calculated. The smaller energy gap (ΔE) for all lead compounds suggests that they can easily transit from HOMO to LUMO, which is important for the molecular reactivity. Since the decrease in electronegativity (χ) value is proportional to the increase in inhibitive efficiency (Zhan et al., 2003), these leads would have higher inhibitory activity because of their lower electronegativity value. The statistical values of the calculated molecular descriptors are tabulated in Table 13. The smaller energy gap, lower electronegativity, and higher dipole moment that are vital for the inhibitory effect of a molecule were observed which validates the better inhibitory activity for the selected lead compounds.

TABLE 13
www.frontiersin.org

TABLE 13. The statistical results of the DFT-based descriptors for the selected lead compounds and the drug.

Conclusion

Pharmacophore modeling, virtual screening, and molecular docking are the rational methods for the identification of novel leads with diverse chemical scaffold. Therefore, ligand-based pharmacophore modeling combined with virtual screening and docking was applied in this study to discover novel, potent, and selective virtual hits for JAK1 enzyme. Initially, the ligand-based pharmacophore models were generated and validated using the potency and selectivity validation methods. Eight pharmacophore models were selected and taken for pharmacophore-based virtual screening against six databases. The hits obtained from screening were filtered through ADME prediction and molecular docking. The binding free-energy calculation and induced fit docking methods were employed to validate the docking results. Subsequently, cross docking was carried out to identify the lead compounds that are more selective toward JAK1. Finally, the top five lead compounds were selected and taken for molecular dynamics and the DFT study. Among the five compounds, T5923555 and T5923531 were found to be the best leads and can be further validated using in vitro and in vivo methods.

Data Availability Statement

The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding authors.

Author Contributions

SB and TM have designed the experiment. SB has contributed to data acquisition, analysis, and interpretation. SB, SN, and SS validated the results. SB has drafted the manuscript. VN, HS, and TM contributed to the critical examination of the manuscript.

Funding

This research was supported by Start-Up Research Grant for Young Scientist (SB/YS/LS-128/2013), funded by the Science and Engineering Research Board (SERB), Department of Science and Technology (DST), Government of India.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors, and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Acknowledgments

This study was supported by the Computational Biology Lab, funded by the SERB Young Scientist grant (SB/YS/LS-128/2013). Author SB thanks CSIR, New Delhi, India, for providing Senior Research Fellowship (SRF). HS thanks the National Research Foundation of Korea (NRF) grant funded by the Korean Government (MSIT) (No. 2021R1F1A1062300) for the research support. This study was also supported by a research fund from Chosun University, 2021. The authors also express gratitude to StemOnc R&D Private Ltd., for the valuable revision of the manuscript. The authors thank the High-Performance Computing Centre, SRM Institute of Science and Technology, for providing the computational facility. This work is part of a PhD thesis submitted by SB to SRM Institute of Science and Technology.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fphar.2022.837369/full#supplementary-material

References

Abraham, M. J., Murtola, T., Schulz, R., Páll, S., Smith, J. C., Hess, B., et al. (2015). GROMACS: High Performance Molecular Simulations through Multi-Level Parallelism from Laptops to Supercomputers. SoftwareX 1-2, 19–25. doi:10.1016/j.softx.2015.06.001

CrossRef Full Text | Google Scholar

Babu, S., Kulkarni, S. A., Sohn, H., and Madhavan, T. (2015). Identification of Leads through In Silico Approaches Utilizing Benzylthio-1h-Benzo[d]imidazol-1-Yl Acetic Acid Derivatives: A Potent CRTh2 Antagonist. J. Mol. Struct. 1102, 25–41. doi:10.1016/j.molstruc.2015.08.031

CrossRef Full Text | Google Scholar

Becke, A. D. (1998). A New Inhomogeneity Parameter in Density-Functional Theory. J. Chem. Phys. 109 (6), 2092–2098. doi:10.1063/1.476722

CrossRef Full Text | Google Scholar

Bottos, A., Gotthardt, D., Gill, J. W., Gattelli, A., Frei, A., Tzankov, A., et al. (2016). Decreased NK-Cell Tumour Immunosurveillance Consequent to JAK Inhibition Enhances Metastasis in Breast Cancer Models. Nat. Commun. 7 (1), 12258–12312. doi:10.1038/ncomms12258

PubMed Abstract | CrossRef Full Text | Google Scholar

Caspers, N. L., Han, S., Rajamohan, F., Hoth, L. R., Geoghegan, K. F., Subashi, T. A., et al. (2016). Development of a High-Throughput crystal Structure-Determination Platform for JAK1 Using a Novel Metal-Chelator Soaking System. Acta Crystallogr. F Struct. Biol. Commun. 72 (11), 840–845. doi:10.1107/S2053230X16016356

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, W., Wu, G., Zhu, Y., Zhang, W., Zhang, H., Zhou, Y., et al. (2019). HOXA10 Deteriorates Gastric Cancer through Activating JAK1/STAT3 Signaling Pathway. Cancer Manag. Res. 11, 6625–6635. doi:10.2147/CMAR.S201342

PubMed Abstract | CrossRef Full Text | Google Scholar

Chrencik, J. E., Patny, A., Leung, I. K., Korniski, B., Emmons, T. L., Hall, T., et al. (2010). Structural and Thermodynamic Characterization of the TYK2 and JAK3 Kinase Domains in Complex with CP-690550 and CMP-6. J. Mol. Biol. 400 (3), 413–433. doi:10.1016/j.jmb.2010.05.020

PubMed Abstract | CrossRef Full Text | Google Scholar

Das, D., Koh, Y., Tojo, Y., GhoshMitsuya, A. K. H., and Mitsuya, H. (2009). Prediction of Potency of Protease Inhibitors Using Free Energy Simulations with Polarizable Quantum Mechanics-Based Ligand Charges and a Hybrid Water Model. J. Chem. Inf. Model. 49 (12), 2851–2862. doi:10.1021/ci900320p

PubMed Abstract | CrossRef Full Text | Google Scholar

De Vicente, J., Lemoine, R., Bartlett, M., Hermann, J. C., Hekmat-Nejad, M., Henningsen, R., et al. (2014). Scaffold Hopping towards Potent and Selective JAK3 Inhibitors: Discovery of Novel C-5 Substituted Pyrrolopyrazines. Bioorg. Med. Chem. Lett. 24 (21), 4969–4975. doi:10.1016/j.bmcl.2014.09.031

PubMed Abstract | CrossRef Full Text | Google Scholar

Dixon, S. L., Smondyrev, A. M., Knoll, E. H., Rao, S. N., Shaw, D. E., and Friesner, R. A. (2006). PHASE: a New Engine for Pharmacophore Perception, 3D QSAR Model Development, and 3D Database Screening: 1. Methodology and Preliminary Results. J. Comput. Aided. Mol. Des. 20 (10), 647–671. doi:10.1007/s10822-006-9087-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Duan, J. J., Lu, Z., Jiang, B., Yang, B. V., Doweyko, L. M., Nirschl, D. S., et al. (2014). Discovery of Pyrrolo[1,2-B]pyridazine-3-Carboxamides as Janus Kinase (JAK) Inhibitors. Bioorg. Med. Chem. Lett. 24 (24), 5721–5726. doi:10.1016/j.bmcl.2014.10.061

PubMed Abstract | CrossRef Full Text | Google Scholar

Dugan, B. J., Gingrich, D. E., Mesaros, E. F., Milkiewicz, K. L., Curry, M. A., Zulli, A. L., et al. (2012). A Selective, Orally Bioavailable 1,2,4-Triazolo[1,5-A]pyridine-Based Inhibitor of Janus Kinase 2 for Use in Anticancer Therapy: Discovery of CEP-33779. J. Med. Chem. 55 (11), 5243–5254. doi:10.1021/jm300248q

PubMed Abstract | CrossRef Full Text | Google Scholar

Forsyth, T., Kearney, P. C., Kim, B. G., Johnson, H. W., Aay, N., Arcalas, A., et al. (2012). SAR and In Vivo Evaluation of 4-Aryl-2-Aminoalkylpyrimidines as Potent and Selective Janus Kinase 2 (JAK2) Inhibitors. Bioorg. Med. Chem. Lett. 22 (24), 7653–7658. doi:10.1016/j.bmcl.2012.10.007

PubMed Abstract | CrossRef Full Text | Google Scholar

Friesner, R. A., Murphy, R. B., Repasky, M. P., Frye, L. L., Greenwood, J. R., Halgren, T. A., et al. (2006). Extra Precision glide: Docking and Scoring Incorporating a Model of Hydrophobic Enclosure for Protein-Ligand Complexes. J. Med. Chem. 49 (21), 6177–6196. doi:10.1021/jm051256o

PubMed Abstract | CrossRef Full Text | Google Scholar

Gordon, G. M., Lambert, Q. T., Daniel, K. G., and Reuther, G. W. (2010). Transforming JAK1 Mutations Exhibit Differential Signalling, FERM Domain Requirements and Growth Responses to Interferon-γ. Biochem. J. 432 (2), 255–265. doi:10.1042/BJ20100774

PubMed Abstract | CrossRef Full Text | Google Scholar

Gruber, C. N., Calis, J. J. A., Buta, S., Evrony, G., Martin, J. C., Uhl, S. A., et al. (2020). Complex Autoinflammatory Syndrome Unveils Fundamental Principles of JAK1 Kinase Transcriptional and Biochemical Function. Immunity 53 (3), 672–e11. doi:10.1016/j.immuni.2020.07.006

PubMed Abstract | CrossRef Full Text | Google Scholar

Gschwend, D. A., Good, A. C., and Kuntz, I. D. (1996). Molecular Docking towards Drug Discovery. J. Mol. Recognit 9 (2), 175–186. doi:10.1002/(sici)1099-1352(199603)9:2<175::aid-jmr260>3.0.co;2-d

PubMed Abstract | CrossRef Full Text | Google Scholar

Haan, C., Behrmann, I., and Haan, S. (2010). Perspectives for the Use of Structural Information and Chemical Genetics to Develop Inhibitors of Janus Kinases. J. Cel Mol. Med. 14, 504–527. doi:10.1111/j.1582-4934.2010.01018.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Haan, C., Rolvering, C., Raulf, F., Kapp, M., Drückes, P., Thoma, G., et al. (2011). Jak1 Has a Dominant Role over Jak3 in Signal Transduction through γc-containing Cytokine Receptors. Chem. Biol. 18 (3), 314–323. doi:10.1016/j.chembiol.2011.01.012

PubMed Abstract | CrossRef Full Text | Google Scholar

Halgren, T. A., Murphy, R. B., Friesner, R. A., Beard, H. S., Frye, L. L., Pollard, W. T., et al. (2004). Glide: a New Approach for Rapid, Accurate Docking and Scoring. 2. Enrichment Factors in Database Screening. J. Med. Chem. 47 (7), 1750–1759. doi:10.1021/jm030644s

PubMed Abstract | CrossRef Full Text | Google Scholar

Harikrishnan, L. S., Kamau, M. G., Wan, H., Inghrim, J. A., Zimmermann, K., Sang, X., et al. (2011). Pyrrolo[1,2-f]triazines as JAK2 Inhibitors: Achieving Potency and Selectivity for JAK2 over JAK3. Bioorg. Med. Chem. Lett. 21 (5), 1425–1428. doi:10.1016/j.bmcl.2011.01.022

PubMed Abstract | CrossRef Full Text | Google Scholar

Harpur, A. G., Andres, A. C., Ziemiecki, A., Aston, R. R., and Wilks, A. F. (1992). JAK2, a Third Member of the JAK Family of Protein Tyrosine Kinases. Oncogene 7 (7), 1347–1353.

PubMed Abstract | Google Scholar

Hevener, K. E., Zhao, W., Ball, D. M., Babaoglu, K., Qi, J., White, S. W., et al. (2009). Validation of Molecular Docking Programs for Virtual Screening against Dihydropteroate Synthase. J. Chem. Inf. Model. 49 (2), 444–460. doi:10.1021/ci800293n

PubMed Abstract | CrossRef Full Text | Google Scholar

Homeyer, N., and Gohlke, H. (2012). Free Energy Calculations by the Molecular Mechanics Poisson-Boltzmann Surface Area Method. Mol. Inform. 31 (2), 114–122. doi:10.1002/minf.201100135

PubMed Abstract | CrossRef Full Text | Google Scholar

Hornakova, T., Springuel, L., Devreux, J., Dusa, A., Constantinescu, S. N., Knoops, L., et al. (2011). Oncogenic JAK1 and JAK2-Activating Mutations Resistant to ATP-Competitive Inhibitors. Haematologica 96 (6), 845–853. doi:10.3324/haematol.2010.036350

PubMed Abstract | CrossRef Full Text | Google Scholar

Hurley, C. A., Blair, W. S., Bull, R. J., Chang, C., Crackett, P. H., Deshmukh, G., et al. (2013). Novel Triazolo-Pyrrolopyridines as Inhibitors of Janus Kinase 1. Bioorg. Med. Chem. Lett. 23 (12), 3592–3598. doi:10.1016/j.bmcl.2013.04.018

PubMed Abstract | CrossRef Full Text | Google Scholar

Ioannidis, S., Lamb, M. L., Wang, T., Almeida, L., Block, M. H., Davies, A. M., et al. (2011). Discovery of 5-Chloro-N2-[(1s)-1-(5-Fluoropyrimidin-2-Yl)ethyl]-N4-(5-Methyl-1h-Pyrazol-3-Yl)pyrimidine-2,4-Diamine (AZD1480) as a Novel Inhibitor of the Jak/Stat Pathway. J. Med. Chem. 54 (1), 262–276. doi:10.1021/jm1011319

PubMed Abstract | CrossRef Full Text | Google Scholar

Irwin, J. J., and Shoichet, B. K. (2005). ZINC--a Free Database of Commercially Available Compounds for Virtual Screening. J. Chem. Inf. Model. 45 (1), 177–182. doi:10.1021/ci049714+

PubMed Abstract | CrossRef Full Text | Google Scholar

Irwin, J. J., Sterling, T., Mysinger, M. M., Bolstad, E. S., and Coleman, R. G. (2012). ZINC: a Free Tool to Discover Chemistry for Biology. J. Chem. Inf. Model. 52 (7), 1757–1768. doi:10.1021/ci3001277

PubMed Abstract | CrossRef Full Text | Google Scholar

Itteboina, R., Ballu, S., Sivan, S. K., and Manga, V. (2017). Molecular Modeling-Driven Approach for Identification of Janus Kinase 1 Inhibitors through 3D-QSAR, Docking and Molecular Dynamics Simulations. J. Recept Signal. Transduct Res. 37 (5), 453–469. doi:10.1080/10799893.2017.1328442

PubMed Abstract | CrossRef Full Text | Google Scholar

Jacobson, M. P., Pincus, D. L., Rapp, C. S., Day, T. J., Honig, B., Shaw, D. E., et al. (2004). A Hierarchical Approach to All-Atom Protein Loop Prediction. Proteins 55 (2), 351–367. doi:10.1002/prot.10613

PubMed Abstract | CrossRef Full Text | Google Scholar

Jaime-Figueroa, S., De Vicente, J., Hermann, J., Jahangir, A., Jin, S., Kuglstatter, A., et al. (2013). Discovery of a Series of Novel 5H-Pyrrolo[2,3-B]pyrazine-2-Phenyl Ethers, as Potent JAK3 Kinase Inhibitors. Bioorg. Med. Chem. Lett. 23 (9), 2522–2526. doi:10.1016/j.bmcl.2013.03.015

PubMed Abstract | CrossRef Full Text | Google Scholar

Jain, A. N. (2003). Surflex: Fully Automatic Flexible Molecular Docking Using a Molecular Similarity-Based Search Engine. J. Med. Chem. 46 (4), 499–511. doi:10.1021/jm020406h

PubMed Abstract | CrossRef Full Text | Google Scholar

Jorgensen, W. L., and Duffy, E. M. (2002). Prediction of Drug Solubility from Structure. Adv. Drug Deliv. Rev. 54 (3), 355–366. doi:10.1016/s0169-409x(02)00008-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Keretsu, S., Bhujbal, S. P., and Cho, S. J. (2021b). Molecular Modeling Studies of Pyrrolo[2,3-D]pyrimidin-4-Amine Derivatives as JAK1 Inhibitors Based on 3D-QSAR, Molecular Docking, Molecular Dynamics (MD) and MM-PBSA Calculations. J. Biomol. Struct. Dyn. 39 (3), 753–765. doi:10.1080/07391102.2020.1714483

CrossRef Full Text | Google Scholar

Keretsu, S., Ghosh, S., and Cho, S. J. (2021a). Computer Aided Designing of Novel Pyrrolopyridine Derivatives as JAK1 Inhibitors. Sci. Rep. 11 (1), 23051–23112. doi:10.1038/s41598-021-0236410.1038/s41598-021-02364-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Kim, W. S., Kim, M. J., Kim, D. O., Byun, J. E., Huy, H., Song, H. Y., et al. (2017). Suppressor of Cytokine Signaling 2 Negatively Regulates NK Cell Differentiation by Inhibiting JAK2 Activity. Sci. Rep. 7 (1), 46153–46212. doi:10.1038/srep46153

PubMed Abstract | CrossRef Full Text | Google Scholar

Kleppe, M., Spitzer, M. H., Li, S., Hill, C. E., Dong, L., Papalexi, E., et al. (2017). Jak1 Integrates Cytokine Sensing to Regulate Hematopoietic Stem Cell Function and Stress Hematopoiesis. Cell stem cell 21 (4), 489–e7. doi:10.1016/j.stem.2017.08.011

PubMed Abstract | CrossRef Full Text | Google Scholar

Kulagowski, J. J., Blair, W., Bull, R. J., Chang, C., Deshmukh, G., Dyke, H. J., et al. (2012). Identification of Imidazo-Pyrrolopyridines as Novel and Potent JAK1 Inhibitors. J. Med. Chem. 55 (12), 5901–5921. doi:10.1021/jm300438j

PubMed Abstract | CrossRef Full Text | Google Scholar

Kumar Nagarajan, S., Babu, S., Sohn, H., Devaraju, P., Madhavan, T., and Madhavan, T. (2018). Toward a Better Understanding of the Interaction between Somatostatin Receptor 2 and its Ligands: a Structural Characterization Study Using Molecular Dynamics and Conceptual Density Functional Theory. J. Biomol. Struct. Dyn. 37 (12), 3081–3102. doi:10.1080/07391102.2018.1508368

PubMed Abstract | CrossRef Full Text | Google Scholar

Kumari, R., Kumar, R., and Lynn, A. (2014). g_mmpbsa-A GROMACS Tool for High-Throughput MM-PBSA Calculationsg_mmpbsa-Aa GROMACS Tool for High-Throughput MM-PBSA Calculations. J. Chem. Inf. Model. 54 (7), 1951–1962. doi:10.1021/ci500020m

PubMed Abstract | CrossRef Full Text | Google Scholar

Labadie, S., Barrett, K., Blair, W. S., Chang, C., Deshmukh, G., Eigenbrot, C., et al. (2013). Design and Evaluation of Novel 8-Oxo-Pyridopyrimidine Jak1/2 Inhibitors. Bioorg. Med. Chem. Lett. 23 (21), 5923–5930. doi:10.1016/j.bmcl.2013.08.082

PubMed Abstract | CrossRef Full Text | Google Scholar

Labadie, S., Dragovich, P. S., Barrett, K., Blair, W. S., Bergeron, P., Chang, C., et al. (2012). Structure-based Discovery of C-2 Substituted Imidazo-Pyrrolopyridine JAK1 Inhibitors with Improved Selectivity over JAK2. Bioorg. Med. Chem. Lett. 22 (24), 7627–7633. doi:10.1016/j.bmcl.2012.10.008

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, J., Abel, R., Zhu, K., Cao, Y., Zhao, S., and Friesner, R. A. (2011). The VSGB 2.0 Model: a Next Generation Energy Model for High Resolution Protein Structure Modeling. Proteins 79 (10), 2794–2812. doi:10.1002/prot.23106

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, R. J., Wang, Y. L., Wang, Q. H., Wang, J., and Cheng, M. S. (2015). In Silico Design of Human IMPDH Inhibitors Using Pharmacophore Mapping and Molecular Docking Approaches. Comput. Math. Methods Med. 2015, 418767. doi:10.1155/2015/418767

PubMed Abstract | CrossRef Full Text | Google Scholar

Lipinski, C. A. (2000). Drug-like Properties and the Causes of Poor Solubility and Poor Permeability. J. Pharmacol. Toxicol. Methods 44 (1), 235–249. doi:10.1016/S1056-8719(00)00107-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Lobanov, M. Iu., Bogatyreva, N. S., and Galzitskaia, O. V. (2008). Radius of Gyration Is Indicator of Compactness of Protein Structure. Mol. Biol. (Mosk) 42 (4), 701–706. doi:10.1134/S0026893308040195

PubMed Abstract | CrossRef Full Text | Google Scholar

Lu, A., Zhang, J., Yin, X., Luo, X., and Jiang, H. (2007). Farnesyltransferase Pharmacophore Model Derived from Diverse Classes of Inhibitors. Bioorg. Med. Chem. Lett. 17 (1), 243–249. doi:10.1016/j.bmcl.2006.09.055

PubMed Abstract | CrossRef Full Text | Google Scholar

Lucet, I. S., Fantino, E., Styles, M., Bamert, R., Patel, O., Broughton, S. E., et al. (2006). The Structural Basis of Janus Kinase 2 Inhibition by a Potent and Specific Pan-Janus Kinase Inhibitor. Blood 107 (1), 176–183. doi:10.1182/blood-2005-06-2413

PubMed Abstract | CrossRef Full Text | Google Scholar

Luo, H. J., Wang, J. Z., Huang, N. Y., Deng, W. Q., and Zou, K. (2014). Induced-fit Docking and Virtual Screening for 8-Hydroxy-3-Methoxy- 5H-Pyrido [2,1-c] Pyrazin-5-One Derivatives as Inducible Nitric Oxide Synthase Inhibitors. J. Chem. Pharm. Res. 6 (3), 1187–1194.

Google Scholar

Lynch, S. M., DeVicente, J., Hermann, J. C., Jaime-Figueroa, S., Jin, S., Kuglstatter, A., et al. (2013). Strategic Use of Conformational Bias and Structure Based Design to Identify Potent JAK3 Inhibitors with Improved Selectivity against the JAK Family and the Kinome. Bioorg. Med. Chem. Lett. 23 (9), 2793–2800. doi:10.1016/j.bmcl.2013.02.012

PubMed Abstract | CrossRef Full Text | Google Scholar

Lyne, P. D., Lamb, M. L., and Saeh, J. C. (2006). Accurate Prediction of the Relative Potencies of Members of a Series of Kinase Inhibitors Using Molecular Docking and MM-GBSA Scoring. J. Med. Chem. 49 (16), 4805–4808. doi:10.1021/jm060522a

PubMed Abstract | CrossRef Full Text | Google Scholar

Menet, C. J., Mammoliti, O., and López-Ramos, M. (2015). Progress toward JAK1-Selective Inhibitors. Future Med. Chem. 7 (2), 203–235. doi:10.4155/fmc.14.149

PubMed Abstract | CrossRef Full Text | Google Scholar

Mysinger, M. M., Carchia, M., Irwin, J. J., and Shoichet, B. K. (2012). Directory of Useful Decoys, Enhanced (DUD-E): Better Ligands and Decoys for Better Benchmarking. J. Med. Chem. 55 (14), 6582–6594. doi:10.1021/jm300687e

PubMed Abstract | CrossRef Full Text | Google Scholar

Nordqvist, C. (2011). Protein JAK Makes Cancer Cells Contract, So They Can Squeeze Out of a Tumor. Medical News Today.

Google Scholar

Pissot-Soldermann, C., Gerspacher, M., Furet, P., Gaul, C., Holzer, P., McCarthy, C., et al. (2010). Discovery and SAR of Potent, Orally Available 2,8-Diaryl-Quinoxalines as a New Class of JAK2 inhibitorsDiscovery and SAR of Potent, Orally Available 2, 8-Diaryl-Quinoxalines as a New Class of JAK2 Inhibitors. Bioorg. Med. Chem. Lett. 20 (8), 2609–2613. doi:10.1016/j.bmcl.2010.02.056

PubMed Abstract | CrossRef Full Text | Google Scholar

Raivola, J., Haikarainen, T., Abraham, B. G., and Silvennoinen, O. (2021). Janus Kinases in Leukemia. Cancers 13 (4), 800. doi:10.3390/cancers13040800

PubMed Abstract | CrossRef Full Text | Google Scholar

Reddy, A. S., Pati, S. P., Kumar, P. P., Pradeep, H. N., and Sastry, G. N. (2007). Virtual Screening in Drug Discovery -- a Computational Perspective. Curr. Protein Pept. Sci. 8 (4), 329–351. doi:10.2174/138920307781369427

PubMed Abstract | CrossRef Full Text | Google Scholar

Saddala, M. S., and Adi, P. J. (2018). Discovery of Small Molecules through Pharmacophore Modeling, Docking and Molecular Dynamics Simulation against Plasmodium Vivax Vivapain-3 (VP-3). Heliyon 4 (5), e00612. doi:10.1016/j.heliyon.2018.e00612

PubMed Abstract | CrossRef Full Text | Google Scholar

Saharinen, P., and Silvennoinen, O. (2002). The Pseudokinase Domain Is Required for Suppression of Basal Activity of Jak2 and Jak3 Tyrosine Kinases and for Cytokine-Inducible Activation of Signal Transduction. J. Biol. Chem. 277 (49), 47954–47963. doi:10.1074/jbc.M205156200

PubMed Abstract | CrossRef Full Text | Google Scholar

Sakkiah, S., Thangapandian, S., John, S., Kwon, Y. J., and Lee, K. W. (2010). 3D QSAR Pharmacophore Based Virtual Screening and Molecular Docking for Identification of Potential HSP90 Inhibitors. Eur. J. Med. Chem. 45 (6), 2132–2140. doi:10.1016/j.ejmech.2010.01.016

PubMed Abstract | CrossRef Full Text | Google Scholar

Sakkiah, S., Krishnamoorthy, N., Gajendrarao, P., Thangapandian, S., Lee, Y. O., Kim, S. M., et al. (2009). Pharmacophore Mapping and Virtual Screening for SIRT1 Activators. Bull. Korean Chem. Soc. 30 (5), 1152–1156. doi:10.5012/bkcs.2009.30.5.1152

CrossRef Full Text | Google Scholar

Sakkiah, S., Thangapandian, S., John, S., and Lee, K. W. (2011). Identification of Critical Chemical Features for Aurora Kinase-B Inhibitors Using Hip-Hop, Virtual Screening and Molecular Docking. J. Mol. Struct. 985 (1), 14–26. doi:10.1016/j.molstruc.2010.08.050

CrossRef Full Text | Google Scholar

Sathe, R. Y., Kulkarni, S. A., Sella, R. N., and Madhavan, T. (2014). Computational Identification of JAK2 Inhibitors: a Combined Pharmacophore Mapping and Molecular Docking Approach. Med. Chem. Res. 24, 1449–1467. doi:10.1007/s00044-014-1223-6

CrossRef Full Text | Google Scholar

Schenkel, L. B., Huang, X., Cheng, A., Deak, H. L., Doherty, E., Emkey, R., et al. (2011). Discovery of Potent and Highly Selective Thienopyridine Janus Kinase 2 Inhibitors. J. Med. Chem. 54 (24), 8440–8450. doi:10.1021/jm200911r

PubMed Abstract | CrossRef Full Text | Google Scholar

Schönberg, K., Rudolph, J., Vonnahme, M., Parampalli Yajnanarayana, S., Cornez, I., Hejazi, M., et al. (2015). JAK Inhibition Impairs NK Cell Function in Myeloproliferative Neoplasms. Cancer Res. 75 (11), 2187–2199. doi:10.1158/0008-5472.CAN-14-3198

PubMed Abstract | CrossRef Full Text | Google Scholar

Schwartz, D. M., Kanno, Y., Villarino, A., Ward, M., Gadina, M., and O'Shea, J. J. (2017). JAK Inhibition as a Therapeutic Strategy for Immune and Inflammatory Diseases. Nat. Rev. Drug Discov. 17 (12), 78–862. doi:10.1038/nrd.2017.20110.1038/nrd.2017.267

PubMed Abstract | CrossRef Full Text | Google Scholar

Sohn, S. J., Barrett, K., Van Abbema, A., Chang, C., Kohli, P. B., Kanda, H., et al. (2013). A Restricted Role for TYK2 Catalytic Activity in Human Cytokine Responses Revealed by Novel TYK2-Selective Inhibitors. J. Immunol. 191 (5), 2205–2216. doi:10.4049/jimmunol.1202859

PubMed Abstract | CrossRef Full Text | Google Scholar

Soth, M., Hermann, J. C., Yee, C., Alam, M., Barnett, J. W., Berry, P., et al. (2013). 3-Amido Pyrrolopyrazine JAK Kinase Inhibitors: Development of a JAK3 vs JAK1 Selective Inhibitor and Evaluation in Cellular and In Vivo Models. J. Med. Chem. 56 (1), 345–356. doi:10.1021/jm301646k

PubMed Abstract | CrossRef Full Text | Google Scholar

Stahl, M., Guba, W., and Kansy, M. (2006). Integrating Molecular Design Resources within Modern Drug Discovery Research: the Roche Experience. Drug Discov. Today 11 (7-8), 326–333. doi:10.1016/j.drudis.2006.02.008

PubMed Abstract | CrossRef Full Text | Google Scholar

Sterling, T., and Irwin, J. J. (2015). ZINC 15--Ligand Discovery for Everyone. J. Chem. Inf. Model. 55 (11), 2324–2337. doi:10.1021/acs.jcim.5b00559

PubMed Abstract | CrossRef Full Text | Google Scholar

Taha, M. O., Atallah, N., Al-Bakri, A. G., Paradis-Bleau, C., Zalloum, H., Younis, K. S., et al. (2008). Discovery of New MurF Inhibitors via Pharmacophore Modeling and QSAR Analysis Followed by In-Silico Screening. Bioorg. Med. Chem. 16 (3), 1218–1235. doi:10.1016/j.bmc.2007.10.076

PubMed Abstract | CrossRef Full Text | Google Scholar

Taldaev, A., Rudnev, V. R., Nikolsky, K. S., Kulikova, L. I., and Kaysheva, A. L. (2022). Molecular Modeling Insights into Upadacitinib Selectivity upon Binding to JAK Protein Family. Pharmaceuticals 15 (1), 30. doi:10.3390/ph15010030

CrossRef Full Text | Google Scholar

Thoma, G., Nuninger, F., Falchetto, R., Hermes, E., Tavares, G. A., VangrevelingheZerwes, E. H. G., et al. (2011). Identification of a Potent Janus Kinase 3 Inhibitor with High Selectivity within the Janus Kinase Family. J. Med. Chem. 54 (1), 284–288. doi:10.1021/jm101157q

PubMed Abstract | CrossRef Full Text | Google Scholar

Vazquez, M. L., Kaila, N., Strohbach, J. W., Trzupek, J. D., Brown, M. F., Flanagan, M. E., et al. (2018). Identification of N-{cis-3-[Methyl(7H-pyrrolo[2,3-d]pyrimidin-4-yl)amino]cyclobutyl}propane-1-sulfonamide (PF-04965842): A Selective JAK1 Clinical Candidate for the Treatment of Autoimmune Diseases. J. Med. Chem. 61 (3), 1130–1152. doi:10.1021/acs.jmedchem.7b01598

PubMed Abstract | CrossRef Full Text | Google Scholar

Vyas, V., Jain, A., Jain, A., and Gupta, A. (2008). Virtual Screening: a Fast Tool for Drug Design. Sci. Pharm. 76 (3), 333–360. doi:10.3797/scipharm.0803-03

CrossRef Full Text | Google Scholar

Wang, H., Aslanian, R., and Madison, V. S. (2008). Induced-fit Docking of Mometasone Furoate and Further Evidence for Glucocorticoid Receptor 17alpha Pocket Flexibility. J. Mol. Graph Model. 27 (4), 512–521. doi:10.1016/j.jmgm.2008.09.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, H. Y., Cao, Z. X., Li, L. L., Jiang, P. D., Zhao, Y. L., Luo, S. D., et al. (2008). Pharmacophore Modeling and Virtual Screening for Designing Potential PLK1 Inhibitors. Bioorg. Med. Chem. Lett. 18 (18), 4972–4977. doi:10.1016/j.bmcl.2008.08.033

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, T., Duffy, J. P., Wang, J., Halas, S., Salituro, F. G., Pierce, A. C., et al. (2009). Janus Kinase 2 Inhibitors. Synthesis and Characterization of a Novel Polycyclic Azaindole. J. Med. Chem. 52 (24), 7938–7941. doi:10.1021/jm901383u

PubMed Abstract | CrossRef Full Text | Google Scholar

Wen, W., Liang, W., Wu, J., Kowolik, C. M., Buettner, R., Scuto, A., et al. (2014). Targeting JAK1/STAT3 Signaling Suppresses Tumor Progression and Metastasis in a Peritoneal Model of Human Ovarian Cancer. Mol. Cancer Ther. 13 (12), 3037–3048. doi:10.1158/1535-7163.MCT-14-0077

PubMed Abstract | CrossRef Full Text | Google Scholar

Williams, N. K., Bamert, R. S., Patel, O., Wang, C., Walden, P. M., Wilks, A. F., et al. (2009). Dissecting Specificity in the Janus Kinases: the Structures of JAK-specific Inhibitors Complexed to the JAK1 and JAK2 Protein Tyrosine Kinase Domains. J. Mol. Biol. 387 (1), 219–232. doi:10.1016/j.jmb.2009.01.041

PubMed Abstract | CrossRef Full Text | Google Scholar

Xiang, Z., Zhao, Y., Mitaksov, V., Fremont, D. H., Kasai, Y., Molitoris, A., et al. (2008). Identification of Somatic JAK1 Mutations in Patients with Acute Myeloid Leukemia. Blood 111 (9), 4809–4812. doi:10.1182/blood-2007-05-090308

PubMed Abstract | CrossRef Full Text | Google Scholar

Xie, H. Z., Li, L. L., Ren, J. X., Zou, J., Yang, L., Wei, Y. Q., et al. (2009). Pharmacophore Modeling Study Based on Known Spleen Tyrosine Kinase Inhibitors Together with Virtual Screening for Identifying Novel Inhibitors. Bioorg. Med. Chem. Lett. 19 (7), 1944–1949. doi:10.1016/j.bmcl.2009.02.049

PubMed Abstract | CrossRef Full Text | Google Scholar

Xie, X., Wang, X., Shi, X., Zhang, Y., Laster, K. V., Liu, K., et al. (2021). Anwulignan Is a Novel JAK1 Inhibitor that Suppresses Non-small Cell Lung Cancer Growth. J. Cel Mol Med 25 (5), 2645–2654. doi:10.1111/jcmm.16289

CrossRef Full Text | Google Scholar

Yang, C. Y., Sun, H., Chen, J., Nikolovska-Coleska, Z., and Wang, S. (2009). Importance of Ligand Reorganization Free Energy in Protein-Ligand Binding-Affinity Prediction. J. Am. Chem. Soc. 131 (38), 13709–13721. doi:10.1021/ja9039373

PubMed Abstract | CrossRef Full Text | Google Scholar

Yang, S. M., Malaviya, R., Wilson, L. J., Argentieri, R., Chen, X., Yang, C., et al. (2007). Simplified Staurosporine Analogs as Potent JAK3 Inhibitors. Bioorg. Med. Chem. Lett. 17 (2), 326–331. doi:10.1016/j.bmcl.2006.10.062

PubMed Abstract | CrossRef Full Text | Google Scholar

Yeh, Y. T., Ou-Yang, F., Chen, I. F., Yang, S. F., Su, J. H., Hou, M. F., et al. (2007). Altered P-JAK1 Expression Is Associated with Estrogen Receptor Status in Breast Infiltrating Ductal Carcinoma. Oncol. Rep. 17 (1), 35–39. doi:10.3892/or.17.1.35

PubMed Abstract | CrossRef Full Text | Google Scholar

Zak, M., Hurley, C. A., Ward, S. I., Bergeron, P., Barrett, K., Balazs, M., et al. (2013). Identification of C-2 Hydroxyethyl Imidazopyrrolopyridines as Potent JAK1 Inhibitors with Favorable Physicochemical Properties and High Selectivity over JAK2. J. Med. Chem. 56 (11), 4764–4785. doi:10.1021/jm4004895

PubMed Abstract | CrossRef Full Text | Google Scholar

Zak, M., Mendonca, R., Balazs, M., Barrett, K., Bergeron, P., Blair, W. S., et al. (2012). Discovery and Optimization of C-2 Methyl Imidazopyrrolopyridines as Potent and Orally Bioavailable JAK1 Inhibitors with Selectivity over JAK2. J. Med. Chem. 55 (13), 6176–6193. doi:10.1021/jm300628c

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhan, C.-G., Nichols, J. A., and Dixon, D. A. (2003). Ionization Potential, Electron Affinity, Electronegativity, Hardness, and Electron Excitation Energy: Molecular Properties from Density Functional Theory Orbital Energies. J. Phys. Chem. A. 107 (20), 4184–4195. doi:10.1021/jp0225774

CrossRef Full Text | Google Scholar

Zhao, X., Chen, M., Huang, B., Ji, H., and Yuan, M. (2011). Comparative Molecular Field Analysis (CoMFA) and Comparative Molecular Similarity Indices Analysis (CoMSIA) Studies on α(1A)-adrenergic Receptor Antagonists Based on Pharmacophore Molecular Alignment. Int. J. Mol. Sci. 12 (10), 7022–7037. doi:10.3390/ijms12107022

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhong, H., Tran, L. M., and Stang, J. L. (2009). Induced-fit Docking Studies of the Active and Inactive States of Protein Tyrosine Kinases. J. Mol. Graph Model. 28 (4), 336–346. doi:10.1016/j.jmgm.2009.08.012

PubMed Abstract | CrossRef Full Text | Google Scholar

Keywords: JAK1, pharmacophore modeling, virtual screening, molecular docking, molecular dynamics simulation, density function theory

Citation: Babu S, Nagarajan SK, Sathish S, Negi VS, Sohn H and Madhavan T (2022) Identification of Potent and Selective JAK1 Lead Compounds Through Ligand-Based Drug Design Approaches. Front. Pharmacol. 13:837369. doi: 10.3389/fphar.2022.837369

Received: 16 December 2021; Accepted: 07 March 2022;
Published: 21 April 2022.

Edited by:

Leonardo L. G. Ferreira, University of São Paulo, Brazil

Reviewed by:

Prasanth Kumar, Gujarat University, India
Mariana Laureano De Souza, Universidade de São Paulo São Carlos, Brazil

Copyright © 2022 Babu, Nagarajan, Sathish, Negi, Sohn and Madhavan. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

*Correspondence: Thirumurthy Madhavan, thiru.murthyunom@gmail.com; Honglae Sohn, hsohn@chosun.ac.kr

Disclaimer: All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article or claim that may be made by its manufacturer is not guaranteed or endorsed by the publisher.