Unraveling the Molecular Mechanism of Recognition of Selected Next-Generation Antirheumatoid Arthritis Inhibitors by Janus Kinase 1.

Rheumatoid arthritis (RA) is a chronic immune-related condition, primarily of joints, and is highly disabling and painful. The inhibition of Janus kinase (JAK)-related cytokine signaling pathways using small molecules is prevalent nowadays. The JAK family belongs to nonreceptor cytoplasmic protein tyrosine kinases (PTKs), including JAK1, JAK2, JAK3, and TYK2 (tyrosine kinase 2). JAK1 has received significant attention after being identified as a promising target for developing anti-RA therapeutics. Currently, no crystal structure is available for JAK1 in complex with the next-generation anti-RA drugs. In the current study, we investigated the mechanism of binding of baricitinib, filgotinib, itacitinib, and upadacitinib to JAK1 using a combined method of molecular docking, molecular dynamics simulation, and binding free energy calculation via the molecular mechanics Poisson-Boltzmann surface area (MM-PBSA) scheme. We found that the calculated binding affinity decreases in the order upadacitinib > itacitinib > filgotinib > baricitinib. Due to the increased favorable intermolecular electrostatic contribution, upadacitinib is more selective to JAK1 compared to the other three inhibitors. The cross-correlation and principal component analyses showed that different inhibitor bindings significantly affect the binding site dynamics of JAK1. Furthermore, our studies indicated that the hydrophobic residues and hydrogen bonds from the hinge region (Glu957 and Leu959) of JAK1 played an essential role in stabilizing the inhibitors. Protein structural network analysis reveals that the total number of links and hubs in JAK1/baricitinib (354, 48) is more significant than those in apo (328, 40) and the other three complexes. The JAK1/baricitinib complex shows the highest probability of the highest-ranked community, indicating a compact network of the JAK1/baricitinib complex, consistent with its higher stability than the rest of the four systems. Overall, our study may be crucial for the rational design of JAK1-selective inhibitors with better affinity.

Introduction

Cytokines are essential for host defense and immunoregulation.[1,2] They also play a significant role in the immunopathogenesis of autoimmune diseases, such as rheumatoid arthritis (RA).[3] Like other inflammatory diseases, RA is a chronic immune-related condition (primarily of joints), highly disabling and painful, and reduces the quality of life. In the Western population, the estimated prevalence of RA is 1–2%,[4] and the worldwide prevalence is believed to be 1%.[5] Over the last 30 years, many therapies including anti-IL-6,[6] anti-CD20, anti-TNF,[7] and CTLA4-Ig[8] agents have been developed to treat RA. Although these significantly improve care, there are some limitations to maximizing these therapies.[9−12] Despite their success, anti-IL-6, anti-CD20, anti-TNF, and CTLA4-Ig agents that target key cytokines are not entirely effective in all patients. Also, they can cause several potentially life-threatening side effects, including liver and lung damage and decreased ability to fight off infections. Recently, the inhibition of Janus kinase (JAK)-related cytokine signaling pathways using small molecules (inhibitors) has received significant attention recently.[13] JAKs are essential growth factors and intercellular parts of cytokine signaling pathways. The JAK kinase phosphorylates and activates the signal transducers and activators of transcription (STATs), which translocate to the nucleus and induce the transcription of effector genes.[14,15]

The JAK family belongs to nonreceptor cytoplasmic protein tyrosine kinases (PTKs), including four basic mammalian types, JAK1, JAK2, JAK3, and TYK2 (tyrosine kinase 2).[15] JAK1 has received significant attention after being identified as a promising target for anti-RA drugs. In addition, other inflammations are closely related to the JAK/STAT signaling pathway.[16] The human JAK1 protein JH1 (JAK homology 1) domain comprises 290 (842–1132) residues. Similar to all protein kinases, JAK1 is divided into a large C-terminal lobe (residues 957–1132) and a small N-terminal lobe (residues 842–956), which were first observed for PKA (protein kinase A).[17] The small N-terminal lobe is comprised of a five-stranded antiparallel β-sheet (β1−β5),[18] 3/10A-helix, and conserved α-helix (residues 919–931). The conserved ATP-phosphate-binding glycine-rich loop (P-loop) (879RDLGEGHFGKV889) is located between its β1 and β2 sheets (Figure A,B). The C-terminal lobe starts from a conserved hinge region (residues 956–968). Moreover, the large C-lobe of JAK1 contains two significant catalytic areas (1001HRDLAARN1008) and an activation segment (residues 1021–1051). JAK1 can be activated through the phosphorylation of Tyr1034 and Tyr1035, located in the activation segment. The structural adjustment induced by the phosphorylation of the activation segment occurs in the Asp–Phe–Gly (DFG) motif. Kornev et al.[19,20] first uncovered eight hydrophobic residues as a catalytic or C-spine and another four hydrophobic residues that constitute a regulatory or R-spine occurring in the small N-terminal lobe and the large C-terminal lobe (see in Figure S1). These hydrophobic clusters also define the active and inactive states of JAKs.

Figure 1

(A) Typical active Janus kinase 1 (JAK1) structure with different vital segments color-coded; (B) schematic illustration of the JAK1 structure; and (C–F) three-dimensional structures of inhibitors with each heavy atom: baricitinib, filgotinib, itacitinib, and upadacitinib, respectively.

The discovery of JAK kinase inhibitors opens a plethora of treatment options for various autoimmune diseases. Tofacitinib (CP-690550) is a type-I pan-JAK inhibitor that the FDA first approved to treat RA.[2,9,18] In the face of extensive kinome selectivity of tofacitinib, it has limited selectivity within the JAK family of kinases. It led to the discovery of JAK1-specific type-I inhibitors, such as baricitinib (INCB028050),[2] upadacitinib (ABT-494),[2] filgotinib (GLPG0634),[21,22] and itacitinib (INCB039110).[23,24] Baricitinib is used to treat moderately to severely active RA in adults.[13] Upadacitinib is recommended to treat RA and other autoimmune diseases.[25] It is highly specific to JAK1 and shows 58-fold and 74-fold (in vitro) selectivity for JAK1 over JAK3 and JAK2, respectively.[24] Currently, filgotinib and itacitinib are in the phase-III clinical trial. Itacitinib shows >20-fold selectivity for JAK1 over JAK2 and >100-fold selectivity over JAK3 and TYK2.[24]

Unfortunately, the mechanism of recognizing these drug molecules by JAK1 is not known entirely due to the absence of crystal structures of the complex. Therefore, it is necessary to explore the binding mechanisms of these antirheumatoid arthritis inhibitors by JAK1. Herein, we sought to investigate the mechanism of binding of type-I inhibitors, such as baricitinib, filgotinib, itacitinib, and upadacitinib to JAK1 by employing comprehensive computational protocols involving molecular docking, molecular dynamics (MD) simulations, and free energy calculations.

Results and Discussion

We combined molecular dynamics simulations and energetic analysis to elucidate the mechanism underlying the binding of baricitinib, filgotinib, itacitinib, and upadacitinib to the JAK1 kinase. The production simulations of these systems were stable in terms of the total and potential energies (data not shown) and the root-mean-square deviation (RMSD) from the corresponding initial structure. Here, we divided the last 100 ns trajectory into four equal segments of 25 ns trajectory and calculated the average RMSD for each segment (see Table S1). We found that the average RMSD for successive segments lies within the standard deviation (SD). This suggests the convergence of each simulation during the last 100 ns. The same trend was observed for replica simulations also. We also evaluated the stability of the hallmark salt bridge between K908@NZ and E925@CD (see Figure S2).

Structural Stability and Flexibility Analysis of JAK1

The temporal root-mean-square deviations (RMSDs) of the backbone atoms of the JAK1 kinase for the corresponding initial equilibrated structure for all systems were calculated and are shown in Figure A. A similar trend was also observed for the replica simulation (see Figure S3). The average RMSD values are reported in Table . The last 100 ns of both runs were considered for the calculation of the mean. It is evident from Figure A that all complex systems, including apo, reached equilibrium within the initial 100 ns of the simulation. The average RMSD value of the apo system was estimated to be 2.10 Å, and a similar value was obtained for all of the complex systems. The superimposed initial and final conformations of each complex are also shown in Figure S4. Next, we investigated the structural stability of the binding pocket by evaluating the time-series RMSDs of the binding pocket (5 Å radius around the inhibitor) residues (Figure B). The corresponding average RMSD values are reported in Table . In the cases of JAK1/baricitinib and JAK1/filgotinib, the pocket was the most stable (∼0.64 Å), while a slightly larger deviation was obtained for JAK1/itacitinib (1.25 Å).

Figure 2

Evolution of the JAK1 (A) protein backbone; (B) root-mean-square deviations (RMSDs) of the backbone atoms of the ligand-binding pocket over simulation time along with the apo system. The residual root-mean-square fluctuations (RMSFs) are shown in (C). Different binding areas are highlighted with a dotted line rectangle box.

Table 1

Average RMSD during the Last 100 ns Production Simulations (Run1 and Run2) of the JAK1/Inhibitor Complexes and Standard Deviations (SD) in the Parenthesesa

RMSD/systems	protein residues	pocket residues	ligand heavy atoms
JAK1/apo	2.10 (0.15)	0.97 (0.11)
JAK1/baricitinib	2.12 (0.27)	0.70 (0.12)	1.46 (0.18)
JAK1/filgotinib	2.26 (0.17)	0.64 (0.14)	0.44 (0.30)
JAK1/itacitinib	1.96 (0.14)	1.25 (0.14)	2.77 (0.16)
JAK1/upadacitinib	2.10 (0.18)	0.91 (0.15)	1.64 (0.13)

All values are in the Å unit.

We also calculated the temporal RMSD of different important regions of JAK1 complexed with these four inhibitors, as shown in Figure S5. The corresponding average values are recorded in Table S2. It is evident from Figure S5 and Table S2 that the P-loop region drifted similarly for all cases except for JAK1/itacitinib. In the case of JAK1/itacitinib, the lowest deviation was obtained (1.37 Å) while varying between 1.97 and 2.17 Å for all other cases, including apo. The catalytic (∼0.40 Å) and hinge (0.29–0.57 Å) regions were stable throughout the simulation for all systems. The same trend was observed for both the C-spine (0.45–0.79 Å) and R-spine (0.50–0.77 Å) regions. In the cases of JAK1/itacitinib and JAK1/upadacitinib, a marginally lower deviation (1.5 Å) was obtained for the activation region compared to the apo system (2.14 Å). Overall, no significant structural rearrangements compared to the apo system were obtained due to the ligand binding.

Subsequently, we measured the residual flexibility of each complex system by estimating the root-mean-square fluctuations (RMSFs) of Cα atoms and compared them with those of the apo system, as shown in Figure C. It is evident from Figure C that all four complexes shared similar trends in RMSF patterns. The nonbinding site regions, including C and N terminals and the activation segment, were more dynamic than the other areas. The P-loop, catalytic loop, and hinge areas were very stable and displayed less flexibility for all complexes. The JAK1/filgotinib system showed higher fluctuations around residues 910–922, 1038, and 1096 compared to other systems.

We also studied the structural compactness of each complex by calculating the radius of gyration (Rg) from the entire production simulation, and the temporal evolution of Rg is shown in Figure S6. The average value of Rg for all systems was found to vary between 19.63 and 19.98 Å (see Table S3). This suggests that no significant change in Rg occurred due to the ligand binding. In addition, we also calculated the solvent-accessible surface area (SASA) for each system, including ligands (see Figure S7A,B), and the corresponding mean SASA values are displayed in Table S3. The average SASA values for all systems range from 13613.88 to 14241.07 Å2. The JAK1/baricitinib complex system showed a relatively higher SASA value (14241.07 Å2) than the other three complexes, including the apo JAK1 kinase. Further, it is evident from Table S3 that the JAK1/itacitinib complex displayed the lowest SASA value. On the other hand, the lowest SASA was obtained for upadacitinib (55.63 Å2), while the highest SASA value was observed for filgotinib (174.50 Å2).

Ligand Dynamics

We investigated the conformational dynamics of each ligand in the JAK1-bound state by calculating the corresponding RMSD of heavy atoms, as shown in Figure A. We aligned the protein with Cα atoms and then calculated the RMSD of the ligand. The average value of the RMSD is listed in Table . Baricitinib and upadacitinib displayed a similar trend in the RMSD, and the mean RMSD is almost the same for these two ligands. Itacitinib showed the maximum positional deviation inside the binding pocket during the simulation, and the average deviation was found to be 2.77 Å. On the other hand, filgotinib showed the least deviation (0.44 Å) during the simulation. Moreover, we compared the JAK1/baricitinib complex structures with the crystal structure of JAK2/baricitinib (PDB ID: 6VN8, resolution: 1.9 Å) and an RMSD value of 1.18 Å for baricitinib (see Figure S8).

Figure 3

Time evolution of the RMSD of ligand heavy atoms (A) and the potential of mean force (PMF) for all four inhibitors (B), with RMSD as a reaction coordinate. The PMF is obtained by calculating the free energy profile along with their respective RMSD values of inhibitors.

Subsequently, we generated the potential of mean force (PMF) for each ligand with respect to the corresponding ligand RMSD, as shown in Figure B. The PMF is generated to determine the relatively low energy states sampled during the production simulation. It is evident from Figure B that the PMF profile is characterized by a single narrow peak for baricitinib, filgotinib, and upadacitinib. However, the PMF profile of itacitinib was wide and characterized by a single global minimum and two additional secondary minima. As shown in Figure B, the free energy minimum values were observed at ∼1.5, ∼0.4, ∼2.9, and ∼1.5 Å for JAK1/baricitinib, JAK1/filgotinib, JAK1/itacitinib, and JAK1/upadacitinib, respectively. In the case of the JAK1/itacitinib complex, two additional secondary local minima were found at ∼1.0 and ∼1.7 Å.

Next, to understand the ligand movement inside the binding cavity, we generated the free energy landscape (FEL) for each ligand, as shown in Figure S9. We considered the ligand RMSD and radius of gyration (Rg) as reaction coordinates for the FEL construction. We identified each minimum energy conformation with the K-mean[26] clustering algorithm, as shown in Figure S9. As evident from Figure S9, baricitinib and filgotinib exhibited two probable conformations in the cavity, and the probabilities of these conformations were 68 and 32% for baricitinib, while the probabilities were 83 and 17% for filgotinib. Similarly, itacitinib and upadacitinib exhibited the three most probable structures, and the sampling spaces were relatively more expansive than the other two. The probabilities of these three conformations occurring were estimated to be 56, 22, and 22% for itacitinib and 78, 11, and 11% for upadacitinib. Overall, our study suggests that itacitinib and upadacitinib are more flexible inside the binding cavity than the other two ligands.

Correlation and Anticorrelation Motions

To study the internal dynamics of JAK1, the dynamic cross-correlation maps (DCCM) for all systems were generated from the corresponding MD trajectory and are displayed in Figure A–E. The DCCM analysis uncovers the relationships between residues and distinct areas by quantifying their relative movement. Highly correlated movement regions or highly positive movement regions are described by yellow to green colors shown in the color scale given in the figure. In contrast, red to sky blue colors indicate the anticorrelation movement regions or highly negative movement regions and zero strength region on the color bar corresponding to no correlation movement regions. The positive and negative values are associated with the same and opposite displacements of residues, respectively. It is evident from Figure A–D that R1, R2, and R3 regions displayed a notable strong negative correlation in apo, JAK1/baricitinib, JAK1/filgotinib, and JAK1/itacitinib complexes, and a moderate negative and positive correlation was observed for the JAK1/upadacitinib complex (Figure E). The catalytic and activation segment displayed a negative correlation concerning the P-loop region for JAK1 complexed with baricitinib, filgotinib, and itacitinib (Figure B–D).

Figure 4

Cross-correlation matrices of the coordinates’ fluctuations for Cα atoms around their mean positions after the equilibrium of production simulation and shows the relative motions of protein residues. The extent of correlated motions and anticorrelated motions are color-coded: (A) apo system, (B) JAK1/baricitinib, (C) JAK1/filgotinib, (D) JAK1/itacitinib, and (E) JAK1/upadacitinib.

In contrast, the JAK1/upadacitinib complex showed a moderate negative or anticorrelation movement (Figure E). The notable strong anticorrelation in apo, JAK1/baricitinib, JAK1/filgotinib, and JAK1/itacitinib complexes demonstrate the direct impact of JAK1 on the overall structure, further enhancing the flexibility of the overall conformation. In the region R5 and its surroundings, the positively correlated motion increases for upadacitinib, and the hinge region or the active site area displays no correlation.

The diagonal movements are highly positive for each system, and off-diagonals show different directional movements after the inhibitor binding. Jonniya et al.[27] also reported that the anticorrelated movement increases for the WNK–ligand complexes after inhibitor binding. However, different inhibitors produced different motion modes, which is clear from this study. In the R4 region, the anticorrelation motions are present in JAK1/itacitinib, while the anticorrelated motions diminish in JAK1/baricitinib, JAK1/itacitinib, and JAK1/upadacitinib complexes. Overall, compared to the apo system, all complex systems display a decrease in the anticorrelated movements.

Principal Component Analysis (PCA) and Dynamics of Free Energy Landscapes (FELs)

To further characterize the structural variations of the JAK1/inhibitor complexes, principal component analysis (PCA) was performed for all systems to extract the principal modes contributing to the complexes’ functional dynamics. The MD trajectories of JAK1 from its free simulation and complexes were employed in PCA by taking the Cα atom into account. The PCA results in terms of eigenvalues (3N, 3 × 290 = 870) and eigenvectors were obtained by diagonalizing the covariance matrix for all five systems and displayed in decreasing order (see Figure S10). The first 10 PCs contained 76–86%, while the first two PCs depicted 35–57% of the overall motion for all systems. The PC1 porcupine plot of all systems is shown in Figure S11A–E, where the arrowhead indicates the direction of motion, and the length indicates the amplitude of the motion modes.

Nowadays, the free energy landscape (FEL) is extensively used to explore the structure–function relationship of proteins through the protein lower energy minima conformations during the molecular interaction due to different inhibitor bindings. Qualitative observations of this FEL can indicate the energy barriers across the conformational basins. The K-means[26] cluster analysis was performed in each basin to investigate the relationship between the structure and FEL further. The center of the most populated cluster was chosen as the representative structure of each FEL’s minimum basin. Herein, the FEL is drawn to understand the conformational dynamics of JAK1 after different types of inhibitor bindings. The FEL was generated utilizing the first two principal components (PC1, PC2) as collective variables and shown in Figure A–E. It depicts changes in the free energy minima in apo JAK1 and JAK1 bound to baricitinib, filgotinib, itacitinib, and upadacitinib. The noteworthy difference significantly indicates the transition pathway during the entire simulation trajectories and the time variations of PC1 and PC2, as shown in Figure S12A–E, which gives an idea of the sampled space of each system.

Figure 5

Free energy landscape (FEL) analysis of five systems: (A) the apo system, (B) JAK1/baricitinib, (C) JAK1/filgotinib, (D) JAK1/itacitinib, and (E) JAK1/upadacitinib. The FELs were obtained using the projections of JAK1 Cα-atom position vectors onto the first two principal components as reaction coordinates. As indicated by the color bar, free energy values calculated from the density of distribution are given in kcal/mol.

As shown in Figure B, a single global minimum is observed at the bottom of the FEL of JAK1/baricitinib and shows three distinct local minima, separated by ∼1–2 kBT. The three distinct local minima correspond to the thermodynamic ground state and encompass different configurations. The representative structure (63%) of the most populated basin was extracted and is shown in Figure B. Similarly, the other three systems show a significantly different transition pathway (see Figure S12C–E) compared with JAK1/baricitinib. There is a noticeable difference among them based on free energy minima. The JAK1 bound to filgotinib, itacitinib, and upadacitinib showed two distinct deep energy minima separated by a shallow energy barrier, as shown in Figure C–E. The corresponding structures from each system’s most populated basin were extracted and are displayed in Figure A–E. Therefore, the unequal landscape reflects the multiple native states of protein and different possible binding geometries with the inhibitors.

Protein Structure Network (PSN) Analysis

A PSN depicts a network of nodes and their edges. Amino acids represent nodes, and the edges are represented as interactions among the nodes; it could be short-range or long-range interactions. Hubs are mainly defined as highest degree nodes, and the network community represents the more connected region of nodes.

Network properties such as the number of nodes, links, hub, and community were analyzed for both apo and ligand-bound JAK1. The average network properties are shown in Supporting Information Table S4, and hubs and communities are shown in Figure A–J for all systems. Interestingly, we found a significant difference in the total number of hubs and communities among all five systems. The apo structure possesses 40 hubs, whereas the inhibitor-bound structures, JAK1/baricitinib, JAK1/filgotinib, JAK1/itacitinib, and JAK1/upadacitinib, consist of 48, 33, 29, and 29 hubs, respectively. As expected, the total number of links and hubs in JAK1/baricitinib (354, 48) is more significant than apo (328, 40) and the other three bound structures (see Figure B and Table S4). The JAK1/baricitinib complex shows the highest probability of the highest-ranked community (C1, red), indicating a compact network of the JAK1/baricitinib complex, consistent with its higher stability than the rest of the four systems. The highest-ranked C1 community of apo, baricitinib, filgotinib, itacitinib, and upadacitinib complexes (see Figure F–J) stabilizes the binding pocket and C-terminal regions, which reduces the flexibility in those areas and correlated well with the RMSF results (see Figure D). Hubs near the binding site area are assumed to be essential for the inhibitor bindings. Our previous findings (Figure C,D) suggest that the active site (mainly hinge and catalytic region) of itacitinib and upadacitinib shows higher stability and lower fluctuations, further verified by hubs and community PSN analysis. Moreover, the rearrangement in the network properties suggests that each system has global conformations, which define their sensitivity and selectivity toward inhibitors.

Figure 6

Hubs and community structures of residue contact networks in (A, F) unbound apo, (B, G) baricitinib, (C, H) filgotinib, (D, I) itacitinib, and (E, J) upadacitinib bound to JAK1, respectively. The hubs are represented as spheres centered on the Cα atoms on the JAK1 protein. The community structure describes a fundamental feature of residue contact networks. Colored according to the community nodes; the largest community ranked 1 is shown in red.

The community analysis suggests that the activation segments of apo and JAK1/filgotinib and JAK1/upadacitinib complexes were involved in the first largest community and those of JAK1/baricitinib and JAK1/itacitinib were involved in the second-largest community. The network community of the inhibitor-bound structure suggests that inhibitors induced perturbations in the network connectivity. These comparative PSN results provide insight into the conformational dynamics that could be essential guidance toward JAK1-specific novel drug developments to treat rheumatoid arthritis-like autoimmune diseases.

Energetics of Inhibitor Binding to JAK1

The binding free energy (ΔGbind) of four antirheumatoid arthritis inhibitors against the JAK1 kinase was estimated using the molecular mechanics Poisson–Boltzmann surface area (MM-PBSA) method. The details of the different components of binding affinity of JAK1-inhibitors are summarized in Table and graphically shown in Figure S13. It is evident from Table that van der Waals (ΔEvdW), electrostatic interactions (ΔEelec), and nonpolar desolvation energy (ΔGnp) favor the complex formation, while polar desolvation energy (ΔGpol) and the configurational entropy (−TΔSconf) penalty disfavor the protein–ligand binding. The favorable intermolecular electrostatic interaction is overcompensated by the desolvation of polar solvation free energy. Consequently, the total polar contribution (ΔGpol+elec) of each system is unfavorable to the protein–ligand complexation. It means that the primary binding interaction is the van der Waals interaction.

Table 2

Energetic Components of the Binding Free Energy for JAK1/Inhibitor Complexes in kcal/mol

components	ΔEvdW	ΔEelec	ΔGpol	ΔGnp	ΔGsolva	ΔGpol+elecb	ΔEinternalc	–TΔSd	ΔG′e	ΔGbindsim
Baricitinib
run 1	–45.72	–26.55	46.86	–4.10	42.76	20.31	–72.27	17.49	–29.51	–12.02
run 2	–49.91	–30.90	53.88	–4.24	49.64	22.98	–80.81		–31.17	–13.68
avg	–47.82	–28.73	50.37	–4.17	46.20	21.65	–76.54		–30.34	–12.85
Filgotinib
run 1	–50.53	–29.93	51.47	–4.79	46.68	21.54	–80.46	17.53	–33.78	–16.25
run 2	–50.32	–29.92	51.51	–4.76	46.75	21.59	–80.24		–33.49	–15.96
avg	–50.43	–29.93	51.49	–4.78	46.72	21.56	–80.35		–33.64	–16.11
Itacitinib
run 1	–63.30	–12.95	45.53	–5.74	39.79	32.58	–76.25	17.88	–36.46	–18.58
run 2	–52.62	–23.02	43.98	–5.24	38.74	20.96	–75.64		–36.90	–19.02
avg	–57.96	–17.99	44.76	–5.49	39.26	26.77	–75.94		–36.68	–18.80
Upadacitinib
run 1	–50.24	–40.37	56.21	–4.16	52.05	15.84	–90.61	19.70	–38.56	–18.86
run 2	–49.05	–41.21	56.12	–4.10	52.02	14.91	–90.26		–38.24	–18.54
avg	–49.65	–40.79	56.16	–4.13	52.03	15.38	–90.43		–38.40	–18.70

ΔGsolv = ΔGnp + ΔGpol.

ΔGpol+elec = ΔEelec + ΔGpol.

ΔEinternal = ΔEvdW + ΔEelec.

ΔS = configuration entropy.

ΔG′ = ΔEinternal + ΔGsolv.

Further, the predicted binding free energy (ΔGbindsim) (see Table ) reveals that upadacitinib (−18.70 kcal/mol) and itacitinib (−18.80 kcal/mol) have the highest binding affinity toward JAK1, followed by filgotinib (−16.11 kcal/mol) and baricitinib (−12.85 kcal/mol). It is encouraging to note that the rank of the experimental affinities[9,15,21,24] of all of the complexes is well-matched with that of our calculated binding free energies. The total polar contributions (ΔGpol+elec) for upadacitinib are less unfavorable than those for the other three inhibitors. The nonpolar desolvation energy ranges from −4.13 to −5.49 kcal/mol for each complex. Table also suggests that ΔEMM (molecular mechanics) and total nonpolar energy (vdW + np) favor the complexation.

In general, biological complexations are opposed by reducing the binding partners’ configurational entropy.[28] Similarly, the inhibitor binding leads to entropy loss, reducing the binding affinity between the protein and inhibitors. For JAK1/inhibitor systems, the configurational entropy varies from 17.5 to 19.7 kcal/mol. However, upadacitinib and itacitinib’s highest affinity compared to the other two is due to the increased favorable electrostatic interactions and van der Waals interactions, respectively. Upadacitinib is more potent against JAK1 compared to baricitinib. This is because, in the case of JAK1/upadacitinib, ΔEvdW (−49.65 kcal/mol) and ΔEelec (−40.79 kcal/mol) are more favorable than those for JAK1/baricitinib (ΔEvdW = −47.82 kcal/mol, ΔEelec = −28.73 kcal/mol). Further, the net polar energy (ΔGpol+elec = 15.38 kcal/mol) is less unfavorable in the case of JAK1/upadacitinib compared to JAK1/baricitinib (ΔGpol+elec = 21.65 kcal/mol). On the other hand, itacitinib displays a better binding affinity than baricitinib and filgotinib due to an increased contribution of ΔEvdW in itacitinib (ΔEvdW = −57.96 kcal/mol) relative to baricitinib (ΔEvdW = −47.82 kcal/mol) or filgotinib (ΔEvdW = −50.43 kcal/mol). Overall, our study suggests that among these four inhibitors, upadacitinib and itacitinib are the most potent inhibitor against JAK1. Mainly, the complexation is driven by the van der Waals interactions, implying that the hydrophobic amino acids in the binding pocket played a vital role in the inhibitor binding. Gao and his co-worker[16] also theoretically investigated the binding mechanisms of a pyrrolotriazine derivative with JAK2 and obtained a similar result as we observed for JAK1 in the current study.

Essential Residues for Inhibitor Binding

To gain a more in-depth insight into the JAK1/inhibitor interaction patterns, the total binding energy was decomposed at the residual level using the MM-GBSA scheme for all complexes and is shown in Figure A–D. It reflects the average energetic contribution of each residue obtained from both MD runs of each complex to the total binding free energy. Each contributing residue (larger than |ΔG| ... 1.0 kcal/mol) along with its energetic components (elec, vdW, pol, np) are reported in Table S5. As shown in Figure A–D, predominantly hydrophobic residues play a significant role in binding and recognizing inhibitor by JAK1. Further, it is revealed from Figure A–D that the overall significant attractive binding contributions come from the P-loop, hinge loop, catalytic segment, and the starting part of the activation segment (A-loop). It can further be noted that the inhibitor/JAK1 interaction spectra of JAK1/baricitinib, JAK1/filgotinib, and JAK1/upadacitinib are similar. However, the interaction spectrum is slightly different for JAK1/itacitinib. This may be due to the asymmetric nature of the inhibitor. Sanachai et al. observed similar trends for the binding of tofacitinib to JAK1/2/3.[29] In the case of JAK1/upadacitinib, one of the HRD motif residues, Arg1007, contributes unfavorably to the binding, which is in contrast to JAK1/baricitinib, JAK1/itacitinib, and JAK1/filgotinib. The same trend was observed for JAK1/tofacitinib.[29]

Figure 7

Decomposition of the binding free energy for the (A) JAK1/baricitinib, (B) JAK1/filgotinib, (C) JAK1/itacitinib, and (D) JAK1/upadacitinib, respectively. Residues with a free energy of more than 1 kcal/mol are labeled here, and the color diagram also shows their actual position in the JAK1 structure.

It is evident from Table S5 and Figure A that in the case of the binding of baricitinib to JAK1, the hotspot residues are Leu1010, Val889, Leu959, Phe958, Leu881, Glu957, and Gly882. The catalytic side residue Leu1010 and hinge residue Val889 contribute to the binding by −2.59 and −2.44 kcal/mol, respectively. For JAK1/filgotinib, the most prominent residues are Leu959, Leu881, Val889, Phe958, Leu1010, Arg1007, and Ala906 (see Table S5 and Figure B). The energetic contributions of the first three hinge residues are −3.60, −2.27, and −2.26 kcal/mol, respectively. Similarly, for JAK1/itacitinib and JAK1/upadacitinib, contributions mainly come from Leu959, Val889, Leu1010, Arg1007, Phe958, Leu881, and Phe886 for itacitinib (Figure C) and from Leu1010, Val889, Leu881, Phe958, Leu959, Gly882, Glu957, and Ala906 for upadacitinib (Figure D). The contributions toward binding of residues Leu959, Val889, and Leu1010 are −2.78, −2.65, and −2.50 kcal/mol, respectively (see Table S6). Similarly, residues Leu1010, Val889, and Leu881 contributed −2.65, −2.61, and −2.28 kcal/mol toward the binding of upadacitinib, respectively (see Table S5).

Furthermore, the detailed comparison of each complex’s per-residual energy components for each complex’s seven most highly contributing residues was analyzed and is shown in Figure S14. The total energy, total polar energy, and the total nonpolar energy of the predicted highly contributing residues for each complex are also shown in Figure S14A–C, respectively. As suggested by Figure S14A, the highest total energy contributing residues are Leu881, Val889, Phe958, Leu959, and Leu1010. The total polar contribution (see Figure S14B) of residues Leu881, Phe886, and Leu1010 disfavor the ligand binding for all cases. The residue Glu957 contributes favorably to baricitinib and upadacitinib and disfavors filgotinib and itacitinib binding. The total nonpolar contributions of these residues contribute favorably to the ligand binding (see Figure S14C).

We also calculated the energetic contribution of the four main loops (P-loop, hinge loop, catalytic loop, and activation loop) for all four complexes (see Figure S15A–D and Table S6). It is evident from Figure S15A that in the case of the P-loop, the total nonpolar (vdW + np), side-chain, and backbone contributed favorably to the binding. However, the total polar (ele + pol) interactions disfavor the ligand binding. Overall, the P-loop highly co-operates with the itacitinib binding in the pocket. Figure S15B suggests that in the case of the hinge loop, all of the energetic components, total polar (ele + pol), total nonpolar (vdW + np), and both side-chain and backbone contributed favorably to the ligand binding. However, in the case of filgotinib, the total polar contributions disfavor the binding. Finally, the energetic contribution of the hinge loop is relatively high for filgotinib. In the case of the catalytic loop (see in Figure S15C), we found that the total polar (ele + pol) part opposes the ligand binding. Overall, P-loop and catalytic loop display a similar type of energetic contribution. The activation loop (Figure S15D) yielded different results from the P-loop, hinge, and catalytic loop. Only the total nonpolar interactions (vdW + np) contributed favorably, while other energetic components disfavored the binding. Overall, this loop’s total energetic contribution is unfavorable to the ligand binding.

Hydrogen Bonds and Hydrophobic Interactions

Hydrogen bonds (H-bonds) play a critical role in biochemical reactions, such as protein folding, protein–ligand, protein–lectin, or protein–protein interactions.[30,31] To complement the energetic analysis, we also investigated the H-bond interactions between the protein and inhibitor for all four complexes using the corresponding MD trajectories (see Figure S16A–E). The percentage occupancy of H-bond obtained from the simulations determined the strength and stability of the H-bond interactions. In Table , the average H-bond occupancy obtained from run1 and run2 is recorded. The time evolution of H-bonds is shown in Figure S17A–D. In our simulations, we observed that both Glu957 and Leu959 residues formed strong H-bond interactions with inhibitors. Williams et al.[18] also experimentally found these two direct H-bonds in the case of the pan-JAK inhibitor tofacitinib.

Table 3

Hydrogen Bonds Formed between JAK1 and Inhibitors and the Corresponding Average Distance and Percent Determined Using the Simulation Trajectories of Run1 and Run2

		molecular dynamics
binding couples			occupancy (%)b
acceptor	donor	avg. distance (Å)a	run 1	run 2	avg
JAK1/Baricitinib
Glu957@O	Lig@NAT	2.82	88.92	91.08	90.0
Lig@N1	Leu959@N	2.94	22.03	21.46	21.75
JAK1/Filgotinib
Leu959@O	Lig@N23	2.80	87.63	86.42	87.02
Lig@N22	Leu959@H	2.94	81.70	5.03	43.36
JAK1/Itacitinib
Leu959@O	Lig@N7	2.84	64.78	64.59	64.68
Lig@N8	Leu959@N	2.93	50.54	25.09	37.82
JAK1/Upadacitinib
Glu957@O	Lig@N6	2.82	88.80	87.91	88.36
Asp1021@OD1	Lig@N4	2.85	20.27	22.71	21.49

The hydrogen bonds are determined by an acceptor···donor distance of ≤3.5 Å and an acceptor···H–donor angle of ≥120°.

Occupancy of the hydrogen bond is defined as the percentage of simulation time that a specific hydrogen bond exits.

In JAK1/baricitinib, we observed two H-bond interactions (Glu957@O–Lig@NAT and Lig@N1–Leu959@N) with occupancies of 90.0 and 21.75%, respectively (see Figure S17A and Table ). In an earlier experimental study, it was found that baricitinib formed two H-bonds with JAK2[32] (see Figure S16E). Further, it can be seen from Figure S17A that initially, Arg1007 had a strong H-bond interaction with baricitinib, but after 10 ns, this H-bond broke due to positional deviations. Figure S18A shows that H-bonds and other interactions are also involved in the baricitinib and JAK1 complex. The gatekeeper residue M956 interacts with baricitinib with a Pi–sulfur bond (see Figure S18 and Table S7). Similarly, residues Leu881, Val889, Ala906, and Leu1010 were involved in Pi-alkyl or Pi-sigma hydrophobic interactions with baricitinib. The hydrophobic surface of the binding pocket for all complexes is shown in Figure S19A–D. The residues that participate in the van der Waals interactions are Val938, Lys908, Gly887, Lys888, Gly884, Glu883, Asp1021, Gly882, Arg1007, Asn1008, Ser963, Gly962, and Phe958.

In the case of the filgotinib-bound system, we noticed two strong H-bond interactions (Leu959@O–Lig@N23 and Leu959@H–Lig@N22) with an occupancy higher than 43.36% (see Figure S17B). Residues such as Val889, Leu1010, Ala906, and Leu881 are involved in Pi-alkyl or Pi-sigma hydrophobic interactions (see Table S7), and the hydrophobicity surface is shown in Figure S19B. Filgotinib also interacts with JAK1 via van der Waals interactions and Ser963, Gly882, Met956, Val938, Glu957, Phe958, Gly962, Pro960, Ser961, Lys970, Arg1007, Gly883, Gly884, Glu966, and Lys965 play an essential role in this interaction (see Figure S18B).

As is evident from Figure S16C and Table , the Leu959 residue has two stable H-bonds with itacitinib (Leu959@O–Lig@N7 and Leu959@N–Lig@N8) with occupancies of 64.68 and 37.82%, respectively. As seen in Figure S17C, the temporal evolution of these H-bonds is stable and within the 3.5 Å distance. On the other hand, His885 forms a new H-bond with itacitinib after 175 ns simulation. In this case, Leu881, Ala906, Leu1010, and Val889 are involved in hydrophobic (Pi-alkyl or Pi-sigma) interactions, while Phe886, Lys965, Gly962, Phe958, Met956, Lys908, Gly1020, Asp1021, Val1009, Gly882, and Gly884, are involved in the van der Waals interactions.

The important H-bonds between JAK1 and upadacitinib are reported in Table and are shown in Figures S16D and S17D. The gatekeeper residues Met956 and Glu883 participate in the Pi–sulfur and halogen interactions, and residues Val938, Val889, Ala906, and Leu1010 are involved in Pi-alkyl or Pi-sigma hydrophobic interactions (see Figure S18D). Similarly, residues Phe958, Gly1020, Ser963, Asn1008, Val1009, Leu964, Arg1007, Gly884, Gly882, and Gly881 are very crucial for the van der Waals and simple carbon–hydrogen bonding interactions, respectively. Overall, the van der Waals interactions stabilize JAK1/inhibitor complexes.

Atomic Contact at the JAK1 Active Site and Solvent Accessibility

We further explored the number of atomic contacts within a sphere of 3.5 Å for each inhibitor during the entire 200 ns simulations, and the same is shown in Figure (left panel). The average number of atom contacts, calculated from the last 60 ns MD simulations, is 17.7 ± 2.7, 23.9 ± 3.4, 14.1 ± 2.3, and 24.5 ± 3.7 for JAK1/baricitinib, JAK1/filgotinib, JAK1/itacitinib, and JAK1/upadacitinib, respectively. This means that upadacitinib showed higher number of atom contacts, followed by filgotinib and baricitinib. It is known that the binding pocket’s water molecules are replaced by inhibitors and could affect protein–ligand interactions.[33] Consequently, to provide insights into the water accessibility of the JAK1 binding pocket, the solvent-accessible surface area (SASA) was calculated by considering a sphere of 3.5 Å around each inhibitor and is shown in Figure (right panel). The average SASA values were found to be 476.2 ± 55.7, 874.8 ± 66.8, 617.4 ± 117.8, and 407.1 ± 54.8 Å2 for JAK1/baricitinib, JAK1/filgotinib, JAK1/itacitinib, and JAK1/upadacitinib complexes, respectively. The SASA calculation indicates that in the case of JAK1/filgotinib, more solvent molecules can access the drug-binding pocket. On the other hand, upadacitinib has less water molecule accessibility, indicating that the binding efficiency of upadacitinib is more significant than the other three inhibitors. Overall, the number of atomic contacts and active site SASA values consistently support the calculated total binding affinity (Table ).

Figure 8

Time evolution of the number of atomic contacts within 3.5 Å of each ligand during the last 200 ns simulations (left panel). The solvent-accessible surface area (SASA) of the residues within 5 Å of each competitive inhibitor (right panel).

Conclusions

Herein, we studied the structural and molecular mechanisms of binding of type-I second-generation competitive inhibitors against the active JAK1 kinase at the atomic level by employing comprehensive protocols, including molecular docking, molecular dynamics simulations, and free energy measurements using the MM/PB(GB)SA method. First of all, the RMSD results suggest that all systems are in good equilibrium during the 200 ns production simulations. The stability of the binding pocket is increased after the inhibitor binding with JAK1. The overall flexibility of residual motions was affected by inhibitor binding, and correlated movements were increased in upadacitinib compared to the other three inhibitor complexes. The PCA analysis suggests that the structural rearrangement has occurred after the binding of inhibitors to JAK1, and these structural changes decrease the hydrogen bond occupancy and affect the distance-dependent interactions. In addition, the PSN analysis on the quaternary structure of JAK1 suggests structural and network changes after inhibitor binding. The MM-PBSA and hydrogen bonds results revealed that the inhibitor’s affinity decreases in the order upadacitinib > itacitinib > filgotinib > baricitinib, which is in good agreement with experimental finding ranks. The selectivity of the inhibitor upadacitinib toward the JAK1 arises mainly due to an increase in total internal energy and a decrease in the size of total polar energy to other inhibitors. The main interacting profile was observed for hydrogen-bonding interactions with hinge region residues Glu957 and Leu959 and some strong hydrophobic interactions with residues Leu881, Val889, Ala906, Val938, Leu1010, Met956, Phe958, and Leu959. Most of the critical binding residues provided the inhibitor stabilizations through the van der Waals interactions, such as Phe886, Lys965, Gly962, Phe958, Met956, Lys908, Gly1020, Asp1021, Val1009, Gly882, and Gly884. Therefore, the structural and energetic study from MD simulations provides significant insights and knowledge into experiments for designing new inhibitors and analogous compounds against JAK1 to treat rheumatoid arthritis.

Materials and Methods

Molecular Docking and System Preparation

When this work was started, there were no complex X-ray crystallographic structures of JAK1 with baricitinib, filgotinib, itacitinib, and upadacitinib in the Protein Data Bank (PDB) (https://www.rcsb.org). Therefore, we used the rapid Lamarckian Genetic algorithm (GA 4.2) search method in the AutoDock 4.2 software suite to prepare the initial complex structures. For the target receptor, JAK1, three-dimensional coordinates were taken from the PDB database (PDB entry: 6AAH; resolution, 1.8 Å).[34] We removed the crystal waters, and the missing residues (913–916) were built using the UCSF Chimera[35] plugin MODELLER.[36] Polar hydrogen atoms and Kollman charges were added by AutoDock Tools (ADT), while nonpolar hydrogen atoms were merged. The ligand atomic coordinates were collected from the PubChem[37] database. The Gasteiger charges and rotatable bonds were assigned to the ligands, and the grid spacing was set at 0.375 Å. AutoDock generated results in the PDBQT format, and optimal geometric conformation with the best binding energy was selected from all possible conformations. Chimera[35] was used to prepare the protein and ligand complex for further molecular dynamics simulations study. The Cartesian coordinates of the docked structures are provided in the Supporting Information. Each ionizable amino acid protonation state was characterized by PROPKA 3.0.[38]

Molecular Dynamics (MD) Simulations

All of the MD simulations were carried out using the AMBER18 MD suite.[39] The force fields used for the ligands and proteins were the updated Generalized Amber Force Field (GAFF2)[40] and ff14SB,[41] respectively. The AM1-BCC[42] atomic charges were calculated for ligands with the antechamber module.[43] For the phosphorylated tyrosine residues, the Phosaa10 force field[44] was used. Each system was solvated using an explicit TIP3P[45] water model, and an octahedron periodic box with a distance cutoff of 10 Å was used. An appropriate number of counterions was added to neutralize the solvated system and maintained a physiological pH of 0.1 M.

Initial minimization of each solvated complex was achieved by performing 500 cycles of steepest descent, followed by another 500 cycles of the conjugate gradient with a harmonic restraint of 2.0 kcal mol–1 Å–2. The second minimization stage was conducted without restraining the solute for 100 cycles of the steepest descent, followed by another 900 cycles of the conjugate gradient. Subsequently, all complexes were gradually heated from 0 to 300 K for 50 ps at the NVT ensemble with a restraint force constant of 2.0 kcal mol–1 Å–2 on the solute.

We performed the density equilibration of each complex with a restraint force of 2.0 kcal mol–1 Å–2 under the NVT ensemble (300 K, 50 ps). Each system was then equilibrated at the NPT ensemble for 1.0 ns without any restrain. Finally, the production simulation was run for 200 ns with a time-step of 2.00 fs under the NPT ensemble. We performed 2 × 200 ns simulations for each system, starting with a different initial velocity. The coordinates were saved for every 10 ps, leading to 20,000 snapshots. During simulations, all bond lengths were constrained, including that of the hydrogen atom, using the SHAKE algorithm.[46] Under the periodic boundary conditions (PBCs), the particle mesh Ewald (PME)[47] method was used to treat the long-range electrostatic interactions with a cutoff of 10.0 Å. The temperature was maintained at 300 K using the Langevin thermostat[48] with a collision frequency of 2.0 ps–1, and the pressure was controlled with the Berendsen Barostat.[49] The postprocessing of trajectories was done using the Cpptraj module[50] of AmberTools19.

Dynamic Cross-Correlation Map (DCCM)

The pairwise cross-correlation coefficients, or dynamic cross-correlation maps,[51] depict the correlation between pairs of atoms. We considered only the Cα atomic coordinates of each residue. The positive and negative values in the map represented the correlated and anticorrelated behavior, respectively (see the Supporting Information). The last 140 ns trajectories of both runs were used for this analysis using Cpptraj.

Principal Component Analysis (PCA) and Free Energy Landscape

To gain insight into the complex motions of all systems of JAK1, the principal component analysis PCA[52] was performed (see the Supporting Information for details). It is a dimensional reduction method and could extract the dominant modes from the whole trajectory. Similar to the DCCM analysis, only Cα atomic coordinates were used for this analysis. The MD trajectories of both runs were employed for the PCA analysis. The cosine content of the first few PCs was obtained via GROMACS (g_covar, g_anaeig, and g_analyze modules).[53] It is used to ascertain the convergence of each trajectory.

We used the first two principal components (PC1 and PC2) as reaction coordinates to draw the free energy landscape (FEL). In three-dimensional FEL, the low-energy “valleys” denote metastable conformational states, while “hills” correspond to the energetic barriers connecting those valleys.

Protein Structure Networks (PSNs)

We used a mixed protein structure network (PSN)[54] and an elastic network model-normal mode analysis (ENM-NMA)-based method to investigate structural communications in the form of network in proteins, which are otherwise not easily detectable. The PSN analysis enumerates network features, like nodes, edges, hubs, links, community, etc., and the shortest path communication from stable states of molecular dynamics (MD) trajectories. In PSN, nodes correspond to the residue, and a hub represents a highly linked residue that generally refers to the degree of a node. Community or modularity is the region in a network where nodes are more connected to each other. Here, network differences were explored for the ligand unbound, apo, and the four ligand-bound JAK1 complexes using WebPSN v2.0.[55] All of the analyses of the hubs and network communities were visualized with VMD.[56]

MM-PBSA Calculations and Ligand-Residue Energy Decomposition Analysis

To evaluate the binding free energy of all the four inhibitors against JAK1, the widely used molecular mechanics Poisson–Boltzmann surface area (MM-PBSA) method[57−61] was employed. In this study, a total of 7000 structural frames in an interval of 2 ps from the last 140 ns trajectory were used for the MM-PBSA calculation. Details of the MM-PBSA methodology are discussed in our previous works,[27,62−65] and here, we followed the same protocol. A brief description of the method is provided in the Supporting Information. The configurational entropy (−TΔS) was calculated using the normal mode analysis (NMA)[66] method (see the Supporting Information). We considered only 130 configurations from the last 140 ns trajectory for the entropy calculation due to the high computational cost. Further, we evaluated the significance of each residue in drug binding by performing the per-residue decomposition of the total binding free energy using the molecular mechanics generalized Born surface area (MM-GBSA) scheme.