DFT study on the structural and chemical properties of Janus kinase inhibitor drug Baricitinib.

Baricitinib is a small molecule used to treat moderate to severe rheumatoid arthritis (RA) in adults. It is an inhibitor of Janus kinase 1 and 2 (JAK1 and JAK2). It has also been repurposed as a potential treatment for Covid 19. The current study has been carried out to understand the structural and chemical properties of this molecule. The molecule is optimized by using density functional theory (DFT) method. The DFT calculations are performed using Gaussian 09 W software package. The bond lengths and bond angles between atoms in the molecules are investigated. The intramolecular interaction within the molecule is identified using the natural bond orbital (NBO) study. The atom in molecule (AIM) study is performed using Multiwfn software. All the calculations are performed at B3LYP /6311G++ (d, p) level of theory. The molecular parameters, such as first-order hyperpolarizability, HOMO-LUMO energy gap, global electrophilicity index, dipole moment, chemical potential, hardness, ionization energy and electron affinity are determined from the calculation. The molecular docking analysis of Baricitinib is also carried out against different target proteins such as 6VSB, 6W9C and 6LU7.

Introduction

Janus kinase (JAK) plays an important role in the pathogenesis of several immune-mediated diseases by enabling the signal transduction of various kinds of cytokines [1]. Cytokines are soluble factors produced and secreted by both immune non-immune cells. The abnormal cytokine production causes loss of immune homeostasis, which can manifest in several immune-mediated inflammatory diseases [2]. Over the years, there has been a widespread effort to recognize and design a small molecule JAK inhibitor to address the unmet medical requirements [3]. Baricitinib is a small molecule that has been approved for the treatment of certain autoimmune and inflammatory disorders. It was first used to treat rheumatoid arthritis and is considered a potential treatment for several dermatologic diseases such as alopecia areata and atopic dermatitis [4]. The recent outbreak of Covid 19 was caused due to the involvement of a cell surface protein, identified as angiotensin-converting enzyme II (ACE2) in the receptor mediated endocytosis for SARS-CoV-2 entry to the cells. It was accompanied by increase in the levels of cytokine signalling through the Janus kinase-signal transducer and activator of transcriptor (JAK-STAT) pathway. At present, there is no drug explicitly identified for Covid 19 and the cytokine storm it causes. However, Baricitinib may interfere with the passage and intracellular assembly of SARS-CoV-2 into the target cells mediated by ACE 2 receptor and also suppress the cytokine storm caused due to covid 19 by suppressing the JAK1/JAK2 [5]. M. J. Hassan et al. [6] compared the clinical results of high dose of Baricitinib with its usual dose in patients with severe covid-19 pneumonia. They observed that the daily intake of increased Baricitinib in severe covid-19 patients results in early stabilization of the respiratory functions, which reduced the need for critical care supports, and decreased the number of rehospitalisations with a mortality rate compared to its usual dose. S. V. Gandhi et. al. [7] analyzed the Baricitinib drug in pure and dosage form by developing a simple and sensitive spectrophotometric method. A. S. Alshetaili [8] studied the solubility of Baricitinib in water, ethanol, polyethylene glycol-400 (PEG-400), ethyl acetate (EA), dichloromethane (DCM) and dimethyl sulfoxide (DMSO) at temperatures between 298.2 K and 323.2 K. They observed that Baricitinib was easily soluble in DMSO and PEG-400, soluble in ethanol and EA, slightly soluble in DCM and poorly soluble in water. Poly lactic-co-glycolic acid (PLGA) nanoparticles of Baricitinib were developed and evaluated by M. J. Ansari et. al. [9] to enhance in-vitro dissolution and performance. They observed that the PLGA polymer showed a more pronounced effect on sustaining the drug release than enhancing the encapsulated drug.

In the current work, we focus on the DFT based study of molecular structure, non-linear optical (NLO) behavior and natural bond orbital (NBO). We have also calculated the HOMO-LUMO energy gap, chemical potential, global electrophilicity index, hyperpolarizability and other thermodynamical properties of the molecule.

Method

The optimization of the molecular structure is performed by implementing DFT method in Gaussian 09 W software at B3LYP level with 6–311++G (d, p) basis set [10]. The topological parameters such as electron density (ρ), the Laplacian of electron density (∇2ρ), the energy density (H(r)) and Potential energy density V(r) at the bond critical points (BCPs) of Baricitinib is determined by using Multiwfn software [11].

Results and discussion

Structural parameters

The molecular structure of Baricitinib with its bond lengths is shown in Fig. 1 . The molecule’s point group is found as C1 with a dipole moment of 3.44 Debye. The self-consistent field (SCF) energy is calculated as −1552.05 Hartree. The optimized bond lengths of the molecule are shown in Supplementary Table 1. In the given molecule, the C—C bond length in the ring is observed to lie between 1.37 and 1.44 Å.The bond lengths of side chains C21-C22, C21-C28 and C21-C27 falls at 1.55 Å while for C5-C14 and C22-C23 it is observed at 1.46 Å. The bondlength for C37-C38 is calculated as 1.52 Å. The strongest C—N bond is observed as 1.15 Å for C23 N24 while the bonds such as N17-C21, C28-N29 and C27-N29 has bond lengths of 1.46 Å and 1.49 Å respectively. The C—N bonds lying within the rings have bond lengths in the range of 1.32- –1.38 Å. The bond N18-N17 shows a bond length of 1.36 Å. The bond lengths for S34 = O35 and S34 = O36 are calculated as 1.47 Å and 1.46 Å respectively. The C—H bond lengths are calculated at 1.08 Å and 1.09 Å respectively while the N9-H11 bond has a bond length of 1.01 Å. The calculated bond angles for the given molecule are shown in Supplementary Table 2.

Fig. 1

Molecular structure of Baricitinib.

Natural bond orbital (NBO) analysis

Natural Bond Orbital (NBO) is an important tool which provides information on intramolecular and intermolecular charge transfer. The interaction between non-Lewis and Lewis orbital and their corresponding stabilization energies can be attained through NBO calculation of the molecule [12].The stabilization energies of the interactions are estimated by using the given equation,

In the above equation E2 represents the stabilization energy, ∊i, ∊j represents the diagonal elements and qi represents the donor orbital occupancies while the diagonal Fock matrix elements are represented by Fij. The high value of stabilization energy direct towards stability and intensive interaction between donor and acceptor [12], [13]. In this work, we have carried out the NBO calculation of the molecule by using NBO 5.0 program under the framework of DFT. The calculated values of NBOs for Baricitinib are shown in Table 1 . The intramolecular interactions within a molecule occurs due to the overlapping between bonding and anti-bonding orbitals which causes intramolecular charge transfer (ICT) resulting in stabilization of the molecule [14]. The charge transfer between bonding orbital σ(C38-H43) and anti-bonding orbital σ*(C23-N24) (7898.57 kJ/mol) is found to be the strongest among all the interactions. The hyper-conjugative interaction of bonding orbitals σ(C37-H39), σ(C37-H40),σ(C37-C38) and σ(C38-H42) with anti-bonding orbital σ*(C23-N24) also shows a very high stabilization energy of 1140.76, 1102.82, 1864.20 and 3423.27 kJ/mol respectively. The high value of E2 (3269.07 kJ/mol) for intramolecular interaction between σ(C37-C38) and σ*(C21-C22) indicates towards its high stability. The orbital interactions such as *(C2-C4) → σ*(C21-C22), n1(O35) → *(C2-C4), n1(O35) → *(C28-H31), σ(C37-C38) → σ*(C21-C22) and σ(S34-C37) → σ*(C21-C22) shows higher stabilization energy values as compared to other interactions. Therefore, these intramolecular charge transfers within the molecule contributes largely towards to the stability of the molecule.

Table 1

Second Order perturbation theory analysis of Fock matrix in NBO basis for Baricitinib.

Sl. No.	Donor (i)	Acceptor (j)	E2 (kJ/mol)	E(j) − E(i) (kcal/mol)	Fij (a.u)
1.	(C2-C4)	*(C5-N6)	55.56	0.16	0.084
2.	(C5-N6)	*(N1-C3)	34.14	0.30	0.0925
3.	(C14-C15)	*(C5-N6)	29.10	0.16	0.064
4.	(C14-C15)	*(C16-N18)	29.30	0.26	0.080
5.	(C16-N18)	*(C14-C15)	60.86	0.26	0.119
6.	(C16-N18)	σ*(C22-H26)	122.74	0.15	0.123
7.	σ(C21-C22)	σ*(C23-N24)	99.85	0.22	0.132
8.	σ(C22-H25)	σ*(C23-N24)	63.34	0.10	0.072
9.	σ(C22-H26)	σ*(C23-N24)	63.47	0.10	0.071
10.	σ(C27-N29)	σ*(C21-C22)	48.05	0.14	0.075
11.	σ(C28-N29)	σ*(C21-C22)	66.57	0.14	0.086
12.	σ(N29-S34)	σ*(C21-C22)	64.37	0.19	0.098
13.	σ(N29-S34)	σ*(C23-N24)	59.60	0.32	0.124
14.	σ(S34-O35)	σ*(C28-H31)	151.09	0.10	0.110
15.	σ(S34-O36)	σ*(C21-C22)	45.06	0.41	0.123
16.	σ(S34-O36)	σ*(C28-H31)	285.33	0.10	0.150
17.	σ(S34-C37)	σ*(C21-C22)	485.55	0.10	0.195
18.	σ(S34-C37)	σ*(C23-H24)	311.06	0.23	0.242
19.	σ(C37-C38)	σ*(C14 = C16)	133.83	0.82	0.297
20.	σ(C37-C38)	σ*(C21-C22)	3269.07	0.08	0.462
21.	σ(C37-C38)	σ*(C23-N24)	1864.20	0.22	0.567
22.	σ(C37-H39)	σ*(C23-N24)	1140.76	0.10	0.302
23.	σ(C37-H40)	σ*(C23-N24)	1102.82	0.10	0.299
24.	σ(C38-H42)	σ*(C23-N24)	3423.27	0.06	0.410
25.	σ(C38-H42)	σ*(C22-H26)	285.54	0.31	0.265
26.	σ(C38-H43)	σ*(C23-N24)	7898.57	0.05	0.545
27.	n1(O35)	*(C2-C4)	1628.26	0.07	0.337
28.	n1(O35)	*(C28-H31)	3302.45	0.06	0.404
29.	*(C2-C4)	*(C23-N24)	2293.76	0.11	0.887
30.	*(C2-C4)	σ*(C21-C22)	860.87	0.31	0.905

Frontier molecular orbital (FMO) analysis

The HOMO and LUMO of a molecule together are called frontier molecular orbital (FMO). HOMO is the highest occupied molecular orbital while LUMO is termed as the lowest unoccupied molecular orbital. The energy difference between HOMO and LUMO gives us the HOMO-LUMO energy gap [15]. The HOMO-LUMO energy gap plays a significant role in determining the bioactivity due to charge transfer within the molecule [16]. A molecule with high HOMO-LUMO energy gap is more stable and less reactive while a low HOMO-LUMO energy gap indicates towards its reactivity. The reactivity of a molecule decreases with increase in the HOMO-LUMO gap [17]. Quantum chemical parameters such as hardness (η), global electrophilicity index (ω), chemical potential (μ), electron affinity and ionization energy can be obtained and elucidated from the calculated HOMO-LUMO energy. The pictorial representation of HOMO-LUMO energy gap is shown in Fig. 2 . The high HOMO-LUMO energy gap (4.57 eV) of Baricitinib in the gas phase direct towards its stability and low reactivity. The chemical potential is calculated as −4.05 eV. The negative value of chemical potential implies that the molecule does not decompose instantaneously into constituent atoms. Thus, the calculated values for HOMO-LUMO gap and chemical potential predicts that the molecule has low reactivity and high stability. The other chemical parameters such as hardness, electron affinity, ionization energy, dipole moment along with thermo-dynamical parameters are shown in Table 2 .

Fig. 2

Frontier Molecular Orbital of Baricitinib.

Table 2

Thermodynamical and chemical parameters of optimized Baricitinib.

Parameters	Values
SCF energy (Hartree)	−1552.05
Total energy (thermal) (kcal mol_1)	221.60
Zero-point vibrational energy(kcal mol_1)	206.99
Rotational constants (GHz)
A	0.40
B	0.08
C	0.08
EHOMO	−6.33
ELUMO	−1.76
ELUMO – EHOMO	4.57
Hardness (η) =1/2 (ELUMO – EHOMO)	2.29
Chemical potential (μ) = 1/2 (ELUMO + EHOMO)	−4.05
IE = −EHOMO	6.33
EA = −ELUMO	1.76
Global electrophilicity index (ω) = μ2/2η	2.74

Analysis of non-linear optical (NLO) behaviour

NLO materials due to its wide application in telecommunication system, switching, optical modulation and other applications in technology has gained the interest of researchers [18]. The non-linear optical behaviour of a material can be predicted by the first order hyperpolarizability. The first order hyperpolarizability is a tensor of rank 3 and has 27 components. These 27 components get reduced to 10 components by Kleinman symmetry [19].

The equation given below is used for evaluating the dipole moment and the first order hyperpolarizability of the molecule.andwhere,

βtot represents the total hyperpolarizability and μ represents the total dipole moment.

The calculated hyperpolarizability expressed in atomic units (a. u) is converted into electrostatic units (e. s. u) (1 a. u = 8.639 × 10−33) as shown in Table 3 .The hyperpolarizability for Baricitinib is found to be 1.51 × 10−30esu which is 7.76 times greater than the standard prototype material urea (0.1947 × 10−30) [19], [20]. It indicates the potential NLO application of the molecule.

Table 3

First order hyperpolarizability and Dipole moment of Baricitinib.

Parameters	Baricitinib
βxxx	−92.62
βyyy	−67.12
βzzz	−15.08
βxyy	26.32
βxxy	7.91
βxxz	−128.43
βxzz	10.94
βyzz	10.32
βyyz	−15.32
βtotal	1.51 × 10−30esu
Dipole moment	3.44 Debye

Atom in molecule (AIM) analysis

The AIM analysis is a powerful means for investigating the nature of intramolecular hydrogen bond. This method requires the value of electron density (ρ) at bond critical point (BCP) and its Laplacian of electron density (∇2ρ) to set the criteria for the existence of hydrogen bond [21]. The energy density at BCP (HBCP) helps in describing the nature of hydrogen bonds more precisely [22].The hydrogen bonds can be divided into three groups: a) for strong H bonds and covalent in nature; ∇2 ρBCP < 0, HBCP < 0b) for medium H bonds and partially covalent in nature; ∇2 ρBCP > 0, HBCP < 0c) for weak H bonds and electrostatic in nature; ∇2 ρBCP > 0, HBCP > 0 [23], [24]. From Table 4 it can be seen that the electron density (ρBCP) of the molecule is 0.0065 which lies within the hydrogen bond range of 0.002–0.04 au [11]. The Laplacian of electron density (∇2 ρBCP) for Baricitinib is calculated to be 0.024 au while the energy density (HBCP) is calculated as 0.13 au. We observe that the values for both ∇2 ρBCP and HBCP are positive indicating that the given molecule is electrostatic in nature.

Table 4

The electron density (ρ), the Laplacian of electron density (∇2 ρ) and the energy density (H) at the bond critical point (BCP) of the hydrogen bond of Baricitinib.

Hydrogen bond	ρBCP	∇2 ρBCP	HBCP
C17-H13-H19	0.0065	0.024	0.13

Molecular docking

The interaction between the ligand and protein receptor can be understood properly by means of molecular docking. The drug activity of the ligand depends on the binding of the drug with receptor site [25], [26]. This study has been carried out for assessing the efficiency of the drug against Covid-19 by targeting different proteins. The docking between ligand and protein receptor was performed using Autodock 4.0 program and visualized through Biovia Discovery Studio 2021. The selected proteins were obtained from RSCB protein data bank. The downloaded proteins are prepared for docking by removal of water. The Kollman and Gasteigner charges and hydrogen atoms are then added as required. The Lamarckian Genetic Algorithm (LGA) feature available in the Auto Dock software is used for the molecular docking process [27]. There are 10 docked conformations obtained with different binding energies for each docking process. The higher negative value of binding energy signifies a better docking process. Therefore, the conformation with the highest negative binding energy for ligand protein interaction is chosen.The pre-fusion spike glycoprotein (PDB ID: 6VSB) is essential for the entry of virus into the host cell while papain like protease PLpro (PDB ID: 6W9C) is significant for the cleavage and maturation of viral polyproteins, assembly of the replicase-transciptase complex and distruption of host responses [28]. The coronavirus main protease Mpro (PDB ID: 6LU7) also plays an important role in mediating the viral replication and transcription of the virus [29]. Therefore, in the present study, molecular docking of Baricitinib is carried out against 6VSB, 6W9C and 6LU7 receptors. The molecular docking results of drug with different target proteins are shown in Fig. 3 (a-c). The binding energy values for 6LU7, 6VSB and 6W9C are calculated as −7.81, −5.71 and −7.33 kcal/mol respectively (Table 5 ). The results suggest that Baricitinib has high binding affinity with the selected receptors. From the 2D representation of the ligand-receptor interaction, it can be observed that the atoms: H11, N24, O35 and O36 forms conventional hydrogen bond interactions with THR190, HIS163, CYS145, SER144 and GLY143 of 6LU7 protein. The binding residues ARG188 and HIS164 forms Pi-Donor Hydrogen bond interactions with H12 and C27 atoms respectively (Fig. 4 a). The amino acid MET 165 interacts with the pyrazole and pyrrole rings through the Pi-Alky interaction. The binding residue GLU166 was also observed to form a Carbon-Hydrogen bond interaction with the pyrazole ring. The amino acids TYR305 and TYR213 in 6W9C receptor interacts with the S34 atom of the ligand via Pi-Sulfur interaction (Fig. 4b). The atoms N18 and H11 shows two conventional hydrogen bond interaction with GLU214 and GLU252. The residues TYR305 and TYR251 forms three Pi—Pi stacked interactions with the pyrrole, pyrimidine and pyrazole rings of the ligand. A Pi-Alkyl interaction also takes place between the binding residue LEU253 and pyrrole ring as shown in Table 5. For the 6VSB receptor, the atom O35 forms a conventional hydrogen bond interaction with the binding residue ASN331 (Fig. 4c). An alkyl interaction is detected between the C23 atom and amino acid PRO521. The pyrazole and pyrimidine rings of Baricitinib forms four Pi-Alkyl interactions with ARG577, PRO579 and VAL576 residues respectively. In addition, ASN544 binding residue forms two Carbon-Hydrogen interactions with the pyrimidine and pyrrole rings.

Fig. 3

a. Interaction of Baricitinib with 6LU7. b. Interaction of Baricitinib with 6W9C.c. Interaction of Baricitinib with 6VSB.

Table 5

Amino acid residue and binding energy analysis of Baricitinib.

Ligand	Receptor (PDB IB)	Binding energy (kcal/mol)	Binding residue	Atoms	Interactions
Baricitinib	6LU7	−7.81	A: THR190	H11	Conventional Hydrogen bond
			A: HIS163	N24	Conventional Hydrogen bond
			A: CYS145	O35	Conventional Hydrogen bond
			A: SER144	O35	Conventional Hydrogen bond
			A: GLY143	O36	Conventional Hydrogen bond
			A: GLU166	Pyrazole	Carbon-Hydrogen bond
			A: MET165	Pyrazole	Pi-Alkyl
			A: MET165	Pyrrole	Pi-Alkyl
			A: ARG188	H12	Pi-Donor Hydrogen bond
			A: HIS164	C27	Pi-Donor Hydrogen bond
	6VSB	−5.71	C: ASN331	O35	Conventional Hydrogen bond
			C: PRO521	C23	Alkyl
			C: ARG577	C16	Pi-Donor Hydrogen bond
			C: ARG577	Pyrazole	Pi-Alkyl
			C: PRO579	Pyrazole	Pi-Alkyl
			C: PRO579	Pyrimidine	Pi-Alkyl
			C: PRO579	N6	Pi-Donor Hydrogen bond
			C: VAL576	Pyrimidine	Pi-Alkyl
			C: ASN544	Pyrimidine	Carbon-Hydrogen bond
			C: ASN544	Pyrrole	Carbon-Hydrogen bond
	6W9C	−7.33	C: GLU214	N18	Conventional Hydrogen bond
			C: GLU252	H11	Conventional Hydrogen bond
			C: TYR305	S34	Pi-Sulfur
			C: TYR305	Pyrazole	Pi-Pi Stacked
			C: TYR213	S34	Pi-Sulfur
			C: TYR251	Pyrrole	Pi-Pi Stacked
			C: TYR251	Pyrimidine	Pi-Pi Stacked
			C: LEU253	Pyrrole	Pi-Alkyl
			C: THR257	C15	Carbon-Hydrogen Bond
			C: SER212	C16	Carbon-Hydrogen Bond

Fig. 4

a. 2D representation of 6LU7-Baricitinib interaction. b. 2D representation of 6W9C-Baricitinib interaction.c. 2D representation of 6VSB-Baricitinib interaction.

Conclusion

In the present study, the molecule is optimized by DFT method at B3LYP/6-311G++(d, p) level of theory. The structural analysis of the molecule is done by analysing the bond lengths and bond angles of the molecule. The intramolecular charge transfer within the molecule is investigated by using NBO 5.0 program. The maximum stabilization energy is obtained for charge transfer between σ(C38-H43) and σ*(C23-N24) with a value of 7898.57 kJ/mol. The negative value of chemical potential (−4.05 eV) and high HOMO-LUMO energy gap (4.57 eV) of the molecules indicates that the molecule is stable. The hydrogen bond identified by AIM analysis shows that the molecule is electrostatic in nature. The molecular docking analysis shows that the drug posses highest binding energy towards the target protein 6LU7 with a value of −7.81 kcal/mol.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.