In-Silico Development of 1, 7-Dihydrodipyrrolo [2,3-b:3’,2’-e] Pyridine -3-carboxamide Derivatives as Candidate Janus Kinase Inhibitors to Control Rheumatoid Arthritis

Janus kinase inhibitors are in the front line among the therapeutics used for rheumatoid arthritis. In the 2016 update of EULAR recommendations, Janus kinase (JAK)

inhibitors had been endorsed for patients with inadequate response to methotrexate (MTX) or biological disease-modifying antirheumatic drugs (Smolen et al., 2017). However, in the 2019 update, the JAK inhibitors have been been given equal weightage to that of biologic drugs based on the efficacy reports of JAK inhibitors in recent years. On failure of the conventional drugs, patients can be switched either to biologic drugs or JAK inhibitors (Harigai, 2019; Taylor, 2019; Smolen et al., 2020). Biologic drugs such as TNFα inhibitors, IL-1 inhibitors, IL-6 inhibitors, and co-stimulatory signal inhibitors entered the rheumatoid arthritis therapeutics with the advent of the twenty-first century and have been very successful in halting the disease progression (Smolen et al., 2017; Camean-Castillo et al., 2019). Although biologic drugs have been a breakthrough, there are some drawbacks associated with the biologics. Biologic drugs produce an appreciable response in a small proportion of patients. A larger proportion of the patients don’t respond to biologics and require a treatment change to second-line biologics (Kotyla, 2018). Immunogenicity, the production of antibodies against the drugs, causes loss of response on prolonged use leading to secondary lack of efficacy while long term medication is needed for biologics to achieve disease remission. Almost 50 - 80% of the patients redevelop the symptoms on discontinuation of biologic drugs (Smolen et al., 2013, 2018; Nishimoto et al., 2014; Emery et al., 2015; Huizinga et al., 2015; Tanaka et al., 2010, 2013, 2015). There is an increased risk of infection and adverse effects related to the intravenous or subcutaneous mode of administration of biologics (Kotyla, 2018). All these factors along with the higher cost of biologic drugs paved the path for small molecule kinase inhibitors.

The JAK inhibitors have been proved to be superior to methotrexate (MTX) (Onuora, 2014; Taylor, 2019). JAK kinases have four members JAK1, JAK2, JAK3, and TYK2. All these Janus kinases along with the STAT (signal transducer and activators of transcription) proteins regulate different immune responses by controlling the JAK-STAT pathways leading to the gene expression or gene silencing (O’Shea and Plenge, 2012). JAK-STAT signaling has been found to regulate the complex relationship between cellular metabolism and inflammation in rheumatoid arthritis (Atzeni et al., 2018; Schwartz et al., 2017). JAK3 is of central importance to control immune responses. It is expressed only in myeloid and lymphoid cells which form the most important part of the immune system (Roskoski, 2016). The importance of JAK3 in the immune response is emphasized by the development of severe combined immunodeficiency (SCID) due to the loss of function of JAK3 (Russell et al., 1994, 1995). The deficiency of JAK3 also has been found to block innate lymphoid cells (ILC) differentiation in the bone marrow.

Tofacitinib, the first FDA approved JAK inhibitor, blocks JAK1/2/3/ and TYK2 (to a lesser extent). It reduces the production of pro-inflammatory cytokines in animal models. In clinical studies, it reduces the production of natural killer cells. In-vivo Tofacitinib has been found to inhibit mitochondrial membrane potential, production of reactive oxygen species by fibroblast-like synoviocytes and important metabolic genes such as Glut-1, PFK3B, PDK1, HK2, and GSK3A. It also inhibits IL-6, IL-8, IL-1, ICAM-1, VEGF, and MMP-1 (McGarry et al., 2018). Tofacitinib has been found to impair the production of IFN-γ from human intraepithelial ILC1 without affecting the production of IL-22. It impairs the ILC1 and ILC3 proliferation and differentiation in vitro (Robinette et al., 2018). In clinical investigations, it improves the disease parameters in adult patients with moderate to severe disease activity (Bergrath et al., 2017). The second FDA approved JAK inhibitor, Baricitinib, is JAK1/2 specific with a lesser effect on JAK3 and TYK2 (Atzeni et al., 2018; Markham, 2017). It has been found to inhibit the expression of nuclear factor-κB ligand and IL-8 regulated neutrophil chemotaxis, two important pathways in RA pathogenesis (Mitchell et al., 2017; Murakami et al., 2017). However, it is being evaluated for safety in long term use. Both these inhibitors have adverse effects associated with them and the FDA has issued a black boxed warning for these (Atzeni et al., 2018; Markham, 2017). So, there is an unmet need for the more efficacious JAK inhibitors with better safety and efficacy profile.

The search for new therapeutics has been accelerated by the advent of computer-aided drug designing techniques. Computer-aided drug designing (CADD) techniques have revolutionized the pharmaceutical industry in the last two decades. Structure-based drug designing and ligand-based drug designing are the two main approaches of CADD. In a structure-based drug designing approach, the structure of the binding pocket of the receptor is analyzed and the drug molecule is designed to develop interactions with the key residues of the binding cavity (Ferreira et al., 2015). The ligand-based techniques are used when the 3D structure of the protein target is not known. It uses the structural information of the already known inhibitors to build new inhibitor molecules. Many computational drug discovery tools have been introduced to assist with the process (Tabeshpour et al., 2018). Here based on the structure-based CADD techniques, we report the identification of 1,7,dihydrodipyrrolopyridine as the hinge binding fragment for the development of new JAK inhibitors with higher predicted efficacy and some selectivity against JAK2. JAK2 is involved in the signaling of many hematopoietic cytokines along with the inflammatory cytokines. Its inhibition has the risk of developing leukopenia, thrombocytopenia, and anemia (Neubauer et al., 1998; James et al., 2005; Scott et al., 2007). So, preference is given to some selectivity against JAK2. However, the FDA approved inhibitor Baricitinib is a JAK1/2 inhibitor (Markham, 2017). Virtual screening has been used to identify the molecular framework showing binding properties and shape complementarity in the JAK3 binding cavity. The modification and docking of selected virtual screening hit resulted in the identification of 1, 7, dihydrodipyrrolopyridine -3-carboxamide based inhibitor molecules. Molecular dynamic simulation studies and binding energy calculations (MM-PBSA) prove the designed inhibitors able to form conformationally stable complexes with JAK3.

The structure of JAK3 in complex with the approved JAK inhibitor ‘Tofacitinib’ was downloaded from the protein data bank (PBD ID: 3LXK). The binding pocket was analyzed using PyMol molecular viewer 2.1.0 and binding interactions were elicited with Discovery studio client (BIOVIA, 2016). PyMol was used to remove the ligand from the enzyme-inhibitor complex to get the structure of the free enzyme in its active form. The protein structure thus obtained was saved as a PDB file. The protein structure was prepared for docking by removing water molecules and adding polar hydrogens using Autodock tools (ADT) (Morris et al., 2009). ChemSketch was used to draw the ligand structures(ACD/ChemSketch (freeware) 2012). These were energy minimized using Avogadro software with the GAFF force field. (Hanwell et al., 2012). Ligands were prepared for docking with ADT. Docking was run using AutoDock Vina with default parameters (Trott and Olson, 2010). The search space was specified at the binding pocket of the JAK3 using grid points. The grid points and docking protocol was validated by running test experiments of docking ‘Tofacitinib’ and other reported inhibitors against JAK3 (Boggon et al., 2005; Chrencik et al., 2010; de Vicente et al., 2014; Farmer et al., 2015; Goedken et al., 2015; Hennessy et al., 2016; Jaime-Figueroa et al., 2013; Lynch et al., 2013; Soth et al., 2013; Thoma et al., 2011). The success of a docking protocol is assessed by its ability to reproduce the crystal structure pose of ligands in the binding pocket. To evaluate the success of the test experiments, the interactions and poses of reported ligands in the binding pocket of JAK3 obtained by docking results were compared with the reported crystal structure poses and interactions. The crystal structure pose as the top-ranked pose was obtained in 76% of the test docking experiments. The literature reported success rate of autodock vina to generate the crystal structure pose is 78% which is very close to our results validating our docking protocols. (Trott and Olson, 2010).

After analyzing all the interacting residues and the docking studies of reported inhibitors, we concluded that the ligands having more than one structural feature to develop hydrogen bonding interactions in the binding pocket (dual binding inhibitors) are better inhibitors. Such inhibitors better fit in the binding pocket and restrict the high energy ligand orientations. Mcule online virtual screening program was used to screen a library of 1000 compounds (Mcule library of purchasable compounds) against JAK3 based on diversity selection and using filters of Lipinski rules which include maximum 10 hydrogen bond acceptors, 5 hydrogen bond donors, the maximum mass of 500 daltons and distribution coefficient less than 5 (Lipinski, 2004). The number of rotatable bonds was set to a maximum of 10. Diversity selection is a workflow of the Mcule database that screens a diverse set of compounds having requested features. The docking workflow was set to Autodock Vina. The virtual screening output was set to be displayed in the descending order of binding affinity values of the hit compounds against JAK3. So, the top 100 ligands were selected for further investigation and were manually docked against JAK3 using Autodock Vina. Our interest was to design a ligand with dual or multiple binding sites. Some ligands were showing a very high value of binding affinity (ΔG) but their pose in the binding pocket did not represent a favorable fit while others had a good pose but lower ΔG value. A favourable fit means the ligand is able to develop appreciable interactions to the ATP binding residues and can adapt the shape complementarity of the binding cavity accomodating its apparent pockets. A compound showing a good compromise between the favourable fit and value of binding affinity was selected. It was then optimized by replacing its different structural components with promising analogue fragments to get a compound with a better fit in the binding cavity, favorable interactions, and appreciable binding affinity. Thus almost 61 different inhibitors were designed and docked against JAK3. The docking results were analyzed using PyMol molecular viewer 2.1.0 and discovery studio client 2016.

The molecular dynamic simulation studies of the designed inhibitors in complex with JAK3 were performed for 10 ns using Gromacs 5.1.5 (Abraham et al., 2015). The topology files of the protein (Janus kinase 3) were prepared with ‘gmxpdb2gmx’, a gromacs builtin tool and AMBER99SB force field was used for it while topology files of the ligands were prepared using Acpype (Sousa da Silva and Vranken, 2012). The topologies and coordinate files of the protein and ligand were merged to create the input files for the protein-ligand complex. TIP3P water model and the cubic box was used to solvate the complex and sodium ions were added to neutralize it. The steepest descent algorithm was used to minimize the system for 100 ps. The rest of the parameters were kept default as given in gromacs tutorial documentation for protein-ligand complex (Lemkul, 2019).

The minimized protein-ligand system was then equilibrated at 300K using a Berendsen thermostat for 100 ps of NVT (constant Number of particles, Volume, and Temperature) ensemble and 100 ps of NPT (constant Number of particles, Pressure, and Temperature) ensemble (using Parrinello-Rahman barostat). Verlet cutoff scheme, XYZ periodic boundary conditions, and Particle Mesh Ewald was used. Van der waals cutoff was set to 1.4 nm. The production run was simulated for 10 ns time-period with each timestep of 0.002 ps and 2 ps interval was used to record trajectories. Results were analyzed using Gromacs built-in tools. Xmgrace was used to draw the plots (Paul, 2015).

The trajectories obtained from molecular dynamic simulations were used to calculate the binding energies of the enzyme-inhibitor complex. Molecular mechanics Poisson−Boltzmann surface area (MM-PBSA) approach was applied to determine the binding energies.‘g_mmpbsa (v 1.6)’ tools were used for the purpose (Kumari et al., 2014). The frames from the trajectories were read every 100 ps to calculate the binding free energy of the enzyme-inhibitor complex. Solvent accessible surface area model (SASA) was used for the calculation of non-polar solvation energy. All other parameters for the determination of polar solvation energy and vacuum potential energy terms were kept default (Kumari, 2019).

A schematic representation of the overall workflow and tools used for drug designing is given in Figure 1.

The vowing features of the binding pocket include a lower pocket defined by LEU 956 and ALA 966, an upper hydrophobic pocket surrounded by main chain β sheets and gatekeeper residue MET 902, a hinge region containing GLU 903, TYR 904 and LEU 905 and a glycine loop at the opening of the binding grove with LYS 830, GLY 831 and ASN 832 (Figure 2). The hydrophobic interactions with CYS 909 from the lower lobe at the base of the binding pocket are important for selectivity against JAK1 and JAK2 as these have serine residue at the equivalent position.

Compound I (Figure 3) was selected from the virtual screening results. It represents a good compromise between the binding pose and the value of binding affinity (ΔG). It was able to accommodate the binding cavity in different orientations. It had good shape complementarity and interactions with the residues from the upper and lower lobe around the binding pocket as well as the residues at the opening of the binding groove. However, it did not form appreciable interactions with the hinge residues, the ATP binding site of the enzyme. This problem was resolved by replacing the benzene ring of the selected compound I, that was facing hinge, by the core of ‘Tofacitinib’ (compound II), the FDA approved JAK inhibitor. It resulted in compound III (having Tofacitinib core) and IV (with slightly modified Tofacitinib core). Both these were docked against JAK3 using Autodock Vina and analyzed for favorable interactions.

Appreciable results obtained for ‘IV’ that formed favorable interactions with the ATP binding hinge region and the residues from the upper and lower lobe (Figure 4), urged us to explore more hinge binding motifs.

The compound ‘IV’ was selected as the lead molecule for the structure-activity relationship (SAR) studies to search for more promising inhibitors. The modification of the compound ‘IV’ was done in four steps. In the first step, the hinge binding motif was replaced by different fragments having hydrogen bond donor and acceptor atoms at appropriate places to engage the ATP binding residues of the binding cavity. In the next step, the piperidine ring of compound ‘IV’ was replaced by different analogues to get the best fit in the front pocket of the binding cavity. It was followed by the modifications of the terminal furan ring and the linker between the furanamide moiety and the piperidine ring. The basic purpose of all these modifications was the search of an inhibitor that can develop interactions with residues, surrounding the binding cavity, from the upper and lower lobe while keeping the interactions with ATP binding residues intact. Better interactions and fit in the binding cavity is associated with better efficacy of the inhibitors.

In the first step of SAR studies, the basic hinge binding motif was replaced by a number of other molecules containing hydrogen bond donor and acceptor atoms at positions able to interact with the ATP binding hinge residues. Once again the Mcule library of molecules was used to get a data set of compounds based on similarity search against pyrrolopyridine molecule, the core fragment of compound ‘IV’. From the results obtained, the compounds having hydrogen bond donors and acceptors at the appropriate position to interact with the ATP binding residues were selected. In this way, 30 different compounds were designed and docked against JAK3 to get a promising binding affinity value, better fit in the binding pocket, and favorable interactions.

The compounds 1, 4, 7, and 8 (Figure 5) reduced the binding affinity value. ‘2’ although restored the value of binding affinity but it did not show any hinge binding pose. ‘3’ had shown a good fit in the binding pocket but the binding affinity did not show any improvement. ‘5’ and ‘6’ show improved score but other than the first pose all other poses place the molecule in the outer region of the binding pocket suggesting it to be flexible enough not to stay in the binding cavity. Compound 9 and 11 have an inverted position in the binding cavity with good hinge binding but as the orientation of the molecule changes, the shape complementarity is lost and it fails to fill the hydrophobic pockets of the binding site that results in reduced binding affinity. Compound 10 produced the highest value for binding affinity but molecular orientation is changed and it does not sit well in the binding pocket instead popes out and makes interactions with residues outside the binding pocket (cause of higher binding affinity) rendering it a poor inhibitor. The cyclopropyl group was introduced in 12 with the view of filling the back hydrophobic pocket but it did not work. In ‘13’ and ‘14’, N, H, and O atoms of the amide group were expected to act as hydrogen bond donor and acceptor to the hinge residues but it did not give any appreciable interactions and fit to the binding pocket. To rule out the possibility that orientation change moves the ligand molecule bit away from the hinge loop and it cannot interact with it, relatively larger groups were introduced as core fragments in compounds 15 and 20. An unexpected orientation change occurred in ‘15’ with the furanamide carbonyl group now facing the hinge and interacting with LEU 905. The furan ring faces the upper (back) hydrophobic pocket where some steric clashes are seen. Compound 16 makes the molecule too large to fit in the binding pocket so a poor inhibitor. In ‘18’, the molecule turns back on itself and fits in the inner part of the binding pocket with no interaction with the outer lip of the binding groove. Compounds 16 and 19 to 24 show lower values of binding affinity. Compound 24 has a truncated five-membered ring to mimic the case as if the orientation of the molecule inverts but no fruitful results are obtained.

Some tricyclic hinge binding fragments were also tried to get the desired fit and shape complementarity. Compounds 25, 26, 29, and 30 orients in such a way in the binding pocket that their hydrogen bond donors and acceptors are not in a position to bind with the hinge residues, and most of the interactions are in the region outside the binding groove. ‘29’ orients such that half of it pops out of the binding pocket. Compound 28 gives a low value of binding affinity. Compound 27, having a linear arrangement of fused tricycles, show good binding affinity value as well as best hinge binding and shape complementarity. The tricyclic ring fills up the back hydrophobic pocket and interacts with GLU 903 and LEU 905 of the ATP binding site through its hydrogen bond donor and acceptor atoms while the piperidine ring interacts to the residues in the front hydrophobic pocket of the binding site. The furanamide ring oxygen makes a hydrogen bond with the ARG 953 from the lower lobe (Figure 6). The hydrogen of the amide group attached to the dipyrrolopyridine moiety makes a hydrogen bond with ALA 966. The dipyrrolopyridine ring makes pi-sigma interactions with the LEU 956 and LEU 828. CYS 909 develops interactions with the piperidine ring while ALA 853 and VLA 836 forms pi-alkyl interactions with the dipyrrolopyridine ring. After analyzing all these results we selected compound 27, with best results so far, for further work. The binding affinity values for all the designed compounds are given in Table 1.

Compound 27 was further investigated to find the effect of the piperidine ring. So it was replaced with different piperidine and pyrimidine isomers and other analogues (Figure 7). In this way, 20 more compounds were designed and docked against JAK3. The results obtained were analyzed for better values of binding affinity, appropriate fit into the binding cavity, and favorable interactions. The binding affinity values for all the compounds obtained by the replacement of the piperidine ring are given in Table 1.

Compounds 31- 38 containing the benzene, pyridine, or pyrimidine ring with different positions of ‘N’ atom did not show any favorable effect. Instead, the presence of these rings reoriented the molecule in the binding pocket with the dipyrrolopyridine moiety facing the outer end of the binding groove and it lost its hinge binding. This may be attributed to the planarity of the ring and electronic conjugation that changed the interaction preferences. Compound 39 has shown the best value and sits well in the binding groove. But surprisingly in ‘40’ and ‘41’, pyrimidine ring is pushed out of the binding pocket. Similarly, compounds 42, 43, and 44 also failed to produce a better binding affinity. Compound 45 with cyclohexyl ring has shown good binding affinity and favorable interactions much similar to that of ‘39’. Overall, the binding mode was better when piperidine and piperazine rings were substituted. This may in part be attributed to the non-planner nature of these rings which bend to better fit in the front hydrophobic pocket. The attempt to substitute five-membered rings was also not successful in producing the required results. This all work suggests the compound 39 and 45 (Figure 8) as the best inhibitors.

In ‘39’ nitrogen of the dipyrrolopyridine ring makes a hydrogen bond with GLU 903 of the hinge loop while the oxygen of the amide group attached to the furanamide ring makes a hydrogen bond with the ASN 832 of the glycine loop. The furan ring oxygen develops hydrogen bonding interactions to the ARG 953 of the lower lobe. The piperidine ring makes hydrophobic interactions with the VAL 836. The core makes pi-sigma hydrophobic interactions with the LEU 828, LEU 956, and alkyl and pi-alkyl hydrophobic interactions with ALA 853, VAL 884, MET 902 and ALA 966. Due to the presence of water-mediated hydrogen bonding between the ‘NH’ of the piperidine ring and nitrogen of the amide group attached to furanamide moiety, the piperidine ring orients to better fit in the front pocket and interacts with the VAL 836 from the upper floor. This orientation change brings the oxygen of the amide group attached to dipyrrolopyridine moiety in close proximity of CYS 909 and a hydrogen bond is established there. This interaction is not seen in compound 27. Although the interaction of dipyrrolopyridine ring with LEU 905 is lost, the other favorable interactions developed compensates for this and better value of binding affinity is achieved.

Compound 45 has a similar binding pose and interactions as that of ‘39’ (Figure 8b). Although the absence of ‘NH’ of piperidine ring changes the orientation and carbonyl attached to the piperidine ring now interacts with LYS 830 from the glycine loop. This carbonyl was involved in intramolecular hydrogen bonding with ‘NH’ of piperidine ring in compound 39. Replacement of the piperidine ring with the cyclohexyl ring in ‘45’ removed such intramolecular hydrogen bonding and the inhibitor molecule is able to adjust according to the shape of the binding site. Three hydrogen bonds have been seen between the dipyrrolopyridine moiety and the hinge residues GLU 903 and LEU 905. Although the strong hydrogen bond with the CYS 909 is lost, the cyclohexyl ring sits very well to fill the front pocket at the floor of the binding cavity and develops pi-alkyl hydrophobic interactions with the CYS 909. The furanamide moiety is oriented in such a way that its amide oxygen and the ring oxygen both make hydrogen bonds with LYS 855 from the upper lobe, the type of interaction not seen in any other inhibitor yet. Dipyrrolopyridine ring is in such a position to develop pi-sigma interactions with LEU 956 and LEU 828 while it forms pi-sulphur interactions with the MET 902. All this makes compound 45 a good inhibitor.

Next, the terminal furan ring was replaced with other nitrogen and oxygen analogues. The fragments were designed to have oxygen and nitrogen (comparable electronegativity) atoms at different positions of the ring to see its effect on binding affinity (Figure 9). The resulting compounds were docked against JAK3 and binding modes were analyzed. In compounds 51, 53, 54, and 55, the piperidine ring is pulled a little away from the front hydrophobic pocket and it’s not able to interact with the pocket. Also, the interactions with LEU 956 and CYS 909 are lost. Compound 52 shows the comparable binding affinity value and interacts with the binding pocket in a similar manner as that of 39 and 45. In compound 56 the orientation change alters the position of the piperidine ring. The hydrogen of the piperidine ring now interacts with the ASP 912 from the lower lobe while the carbonyl attached to piperidine ring hydrogen bonds to the CYS 909. The oxygen of the amide group attached to the furan ring hydrogen bonds to the ARG 911 from the lower lobe. The interactions with LYS 830 and ASN 832 seen in compound 39 and 45 are lost. As a result, compound 56 does not fit well in the binding groove but pops out of the pocket. Some hexacyclic rings were also substituted for terminal furan ring but favorable results were not obtained. Perhaps due to the steric effects, the larger rings are not able to fit well at the lip of binding groove while maintaining the interactions of dipyrrolopyridine moiety with the hinge residues.

The linker group between the piperidine and the furan ring was modified to see the role of its different components in enzyme binding (Figure 10). As the best results were obtained for furan ring so the linkers were substituted for compound 39. The binding affinity values for all the compounds thus obtained are given in Table 1. Smaller lengths of linkers were not favorable for the required interactions and better fit in the binding groove. The removal of the amide group attached to the furan ring reduced the score significantly as it is involved in the interactions of inhibitor with LYS 830 and ASN 832. However, the removal of a carbonyl adjacent to the piperidine ring in ‘57’ did not produce much effect and the interactions observed were similar to those of ‘39’.

All this work resulted in two more potent inhibitors, compound 52 (with pyrrole ring substituted for furan ring) and 57 (with a modified linker), both the modifications of compound 39. Similar modifications were also made for compound 45 but these did not show appreciable results.

In compound 52, the interactions of the core fragment are similar to that of ‘39’. The nitrogen of the dipyrrolopyridine moiety makes a hydrogen bond with GLU 903 of the hinge loop. The piperidine ring makes hydrophobic interactions with the VAL 836 (Figure 11a). The dipyrrolopyridine core makes pi-sigma hydrophobic interactions with the LEU 828, LEU 956, and alkyl and pi-alkyl hydrophobic interactions with ALA 853, VAL 884, MET 902 and ALA 966. The methyl group of the linker, between the piperidine ring and the furanamide moiety, forms carbon-hydrogen non-conventional interactions with the ASN 954. The orientation of the piperidine ring is also similar to what is observed in ‘39’. However, the pyrrole ring in place of furan has changed the orientation at the outer lip of the binding grove. The oxygen of the amide group attached to the terminal ring that was making hydrogen bonds with the ASN 832 of the glycine loop now develops hydrogen bonding interactions with the ARG 953 and ASN 954 of the lower lobe while its NH interacts with LYS 830 from the upper lobe. The orientation of the pyrrole ring is slightly different than that of furan but it does not produce any clashes.

The orientation of the furanamide moiety of ‘57’ is the same as that of ‘39’. The furan ring oxygen makes a hydrogen bond to the ARG 953 of the lower lobe while the oxygen of this amide group makes a hydrogen bond with the ASN 832 of the glycine loop. Due to the absence of carbonyl group attached to the piperidine ring, water-mediated intramolecular hydrogen bonding between it and ‘NH’ of the amide of dipyrrolopyridine moiety is lost and the molecule is in relatively extended form. The dipyrrolopyridine ring is close enough to the hinge binding region to make three hydrogen bonds to the GLU 903 and LEU 905 instead of only one hydrogen bond with the hinge region as was seen in the case of ‘39’. This ring is also close enough to MET 902 to develop pi-sulphur interaction with it. The dipyrrolopyridine moiety makes pi-sigma hydrophobic interactions with the LEU 828, LEU 956, and GLY 908 along with some weak alkyl and pi-alkyl hydrophobic interactions with ALA 853, VAL 884, VAL 836 and ALA 966. As the molecule is in slightly extended form the piperidine ring has lost its interactions with VAL 836 (Figure 11b).

Table 1: Values of binding affinity (obtained from Autodock Vina) for the compounds obtained by the modification of the designed hybrid compound IV.

In surface view, compound 39, 45, 52, and 57 sit well in the apparent pockets of the binding cavity. The dipyrrolopyridine ring fits well in the upper hydrophobic pocket of the binding groove while the piperidinyl/cyclohexyl rings sit well in the front pocket. In 39, 52, and 57, the terminal ring fills the third pocket lying at the glycine loop end while in 45, it goes under the glycine loop. An overlay of the selected inhibitors in the binding cavity of JAK3 has been shown in Figure 12a. The Tofacitinib in the JAK3 binding cavity has been shown for comparison (Figure 12b). Table 2 summarises the interactions developed by the selected inhibitors in the binding cavity of JAK3. The interacting residues for ‘Tofacitinib’ were taken from crystal structure PDB ID: 3LXK

All the selected inhibitors were checked for drug likeliness using Lipinski rules to rule out any violations gained during modification of the lead molecule. All these were found suitable as drug candidates with no violation of Lipinski rules. Some of their parameters are given in Table 3.

All the designed inhibitors were also docked against JAK1, JAK2, and TYK2 to check their selectivity within the JAK family. These were found to have some selectivity against JAK2. So, the work that was started with compound 27 resulted in compounds 39, 45, 52, and 57 as the potential inhibitors of JAK3. Although compound 27 itself was well fit as a JAK3 inhibitor, it was not selective against JAK2. The promising inhibitors ‘52’ and ‘57’ that emerged from this study are (Figure 13) predicted to be more selective for JAK3 and JAK1 over JAK2 as compared to Tofacitinib, the FDA approved JAK inhibitor (Table 4). Tofacitinib was successfully docked into the binding pocket of the JAK1, JAK2, JAK3 and TYK2 reproducing the crystal structure pose with a free energy of binding of -8.9, -8.1, -9.1 and -8.0 respectively which is in accordance with the reports about the selectivity of tofacitinib given in the literature. In solution-phase assays, the IC50 values of tofacitinib have been reported to be 3.2 nM, 4.1 nM, 1.6 nM, and 34.0 nM for JAK1, JAK2, JAK3, and Tyk2, respectively (Meyer et al., 2010). Its functional selectivity for JAK1 and JAK3 against JAK2 in cellular assays has also been reported (Flanagan et al., 2010). This emphasizes the impact of the slight difference in the predicted binding affinity value on the experimental selectivity and thus supports the predicted selectivity and efficacy of our designed inhibitors.

Table 2: Residues interacting with the designed inhibitors within the binding cavity of JAK3.

Table 3: The parameters of the designed inhibitors showing their drug likeliness according to the Lipinski rules.

Table 4: Values of binding affinity in kcal/mol for selected inhibitors.

Molecular dynamic simulations of designed enzyme-inhibitor complexes over 10 ns time period have established the dynamic stability of our designed inhibitors within the binding cavity of the JAK3. It has proved the ability of the inhibitors to retain within the binding cavity of the enzyme during conformational changes. This is supported by the RMSD (root mean square deviation), the radius of gyration (Rg), and the number of hydrogen bonds between the ligand and protein per time frame obtained from trajectories of 10 ns MD simulations.

The steady convergence of RMSD from the initial structure for all the four enzyme-inhibitor complexes and free enzyme (taken as reference) show their conformational stability (Figure 14).

The radius of gyration (Rg) predicts the compactness of the protein structure. An irregular value of Rg for a protein represents its unfolding or denaturation of secondary structure. The radius of gyration for our proposed inhibitor complexes remains steady throughout the simulation period ruling out unfavorable changes in protein structure (Figure 15).

The number of intermolecular hydrogen bonds (within 0.35Å) between the inhibitors and JAK3 in the enzyme-inhibitor complex was obtained through the trajectories of 10 ns simulations. The plots of hydrogen bonds against the time show the formation of three to four hydrogen bonds per time frame (Figure 16). This not only validates the docking results but also confirms the stability of the enzyme-inhibitor complex formed.

Table 5: Values of binding energy for proposed inhibitors as calculated by MM-PBSA.

The binding free energies of the designed inhibitors in complex with JAK3 are calculated to estimate the relative binding affinities. Molecular mechanics Poisson–Boltzmann surface area (MM-PBSA) approach has been used. The binding energies thus calculated by these computational techniques do not necessarily match with the experimental results but gives good agreement with the ranking of a number of ligands (Hou et al., 2011).

The binding energies for the designed inhibitors and the reference drug ‘Tofacitinib’ are given in Table 5. The values indicate the proposed inhibitors to form stable complexes with JAK3 showing better efficacy as compared to the reference drug. These results also validate our docking results.

A number of therapeutics are available for the treatment of rheumatoid arthritis including conventional drugs, biologic drugs, and Janus kinase inhibitors. A few JAK inhibitors have been approved by the FDA for rheumatoid arthritis to deal with moderate to severe disease. There are many adverse effects associated with approved therapies available and a large number of patients remain unresponsive to the available therapies. This shows the search of more efficacious drugs to deal with the disease. Using structure-based CADD techniques, we have developed JAK inhibitors showing promising efficacy. Detailed investigation of the docking results and structural analysis of 1, 7-dihydrodipyrrolo[2,3-b:3’,2’-e]pyridine -3-carboxamide derivatives within the binding groove of the JAK3 led us to conclude the dihydrodipyrrolopyridine -3-carboxamide derivatives as potent candidates for selective JAK3 inhibition. The designed inhibitors show higher potency and better selectivity than Tofacitinib (FDA approved inhibitor) as predicted by the docking results. Molecular dynamic studies confirm the conformational stability of the designed inhibitors within the binding cavity of the Janus kinase 3. The proposed inhibitors 52 and 57 show better selectivity against JAK2 providing a safety window against the risk factors associated with JAK2 inhibition. These can be pursued for synthesis and further cytotoxic studies to evaluate their in-vivo and in-vitro efficacy. The work suggests the non-aromatic rings (piperidine and cyclohexane) at a position appropriate to interact with the front pocket represent a good fit in terms of the shape complementarity. The introduction of polar groups at the outer end of the inhibitor molecule develops appreciable interactions with the glycine loop and lower lobe residues enhancing binding affinity. This work identified 1,7-dihydrodipyrrolopyridine molecule as a novel hinge binding motif to develop ATP competitor inhibitors. The use of 1, 7-dihydrodipyrrolo[2,3-b:3’,2’-e]pyridine ring has shown a boost for efficacy. Although the dipyrrolopyridines have been used in some enzyme inhibitors but in this arrangement their properties have not been explored for drug designing. 1, 7-dihydrodipyrrolo[2,3-b:3’,2’-e] pyridine moiety can adapt with different conformations of the inhibitor in the binding cavity while maintaining the interactions with the hinge ATP binding residues. These structural features can be investigated in combination with other reported molecular features to develop inhibitors against other protein kinases of medicinal importance.

The authors declare that there is no conflict of interest regarding the publication of this article.