By taking, the observed IC50 values of all compounds in M were converted into pIC50 = -Log (1/IC50), which was then used as the dependent variable. For the categorization of the compounds into training and test sets in this research, QSARINS v2.2.4 (Gramatica et al., 2013) was used in the most random manner possible. The data set has a homogenous distribution and includes two sets: a training set (20 compounds, 70%) and a test set (11 compounds, 30%). The PubChem CID number, IUPAC name, and inhibitory activities were listed in Table 1. The open3DALIGN tools’ docking-based alignment was used to superimpose the compounds (Tosco et al., 2011). Compound 4 was chosen as the format to line up other compounds because of its strong autodocking rating and the alignment was finished by the open3DALIGN software as depicted in Figure 1. The calculation of molecular field descriptors was then applied to the compound alignment.

The electrostatic field was determined utilizing a volume-less +1 charge test (Tosco and Balle, 2011), though the steric field was determined utilizing a carbon molecule test for the Open3DQSAR programming (Kumar et al., 2015). Some erroneous variances were eliminated to shrink the partial least square (PLS) framework. As part of the data pre-treatment operation, additional N-level variable removal was carried out (Kumar et al., 2015). These N-level variables posit a distribution of only N values across a constrained number of training set objects. The Smart Region Definition (SRD) algorithm, which is rooted in factor connectedness in 3D space, was then used to group the variables (Pastor et al., 1997). Then, a variety of variable selection techniques, such as Fractional Factorial Design (FFDSEL) and Uninformative Variable Elimination-Partial Least Square (UVE-PLS) variable selection (Baroni et al., 1992; Centner et al., 1996), were used to develop the best PLS models. Cross-validated LOO (Leave-One-Out), LTO (Leave-Two-Out), or LMO (Leave-Many-Out) paradigms were used to compute PLS models (Tosco and Balle, 2011; Abdel Samee et al., 2012; Abdel Samee, 2020; Nawaz et al., 2022). Last but not least, the Maestro graphics package was used to visualize PLS coefficient grid maps or the activity-correlating molecular regions in the form of iso-contour maps.

For the present study, two different crystallized structures of type 1 diabetes (PDB ID: 1JK8 and 1XW7), rheumatoid arthritis (PDB ID: 2AXJ and 2FSE), one for protein tyrosine phosphatase, nonreceptor type 22 (lymphoid) isoform 1 (PTPN22) (PDB ID: 4J51), and one for SARS-CoV-2 main protease (PDB ID: 6LU7) were used to test the validity of Autodocking vina docking program with Assisted Molecular Docking with AutoDock4 and Autodocking Vina (AMDock) v1.5.2 (Valdés-Tresanco et al., 2020).

Molecular docking simulations, as a right approach, were employed using the AMDock v1.5.2 (Valdés-Tresanco et al., 2020) to discover the structural interaction mechanism between ligands and the diabetes type 1, rheumatoid arthritis, and SARS-CoV-2 main protease receptors. The 3D crystallized structure of diabetes type 1, rheumatoid arthritis, PTPN22, and SARS-CoV-2 main protease receptors was made available by Protein Data Bank. To create docking simulations, the docking rule’s setup procedures were used. The co-crystallized ligand and all water molecules were separated from the protein using the Discovery Studio client software, and the binding modes were calculated using Optimal box size 1.1 (Feinstein and Brylinski, 2015) and AutoLigand (Harris et al., 2008) embedded in AMDock software. The AMDock experimental tool with the autodocking Vina method (Trott and Olson, 2010) was employed to find the best-docked ligands with the receptors. The Discovery Studio 2020 client software was used to analyze the results of the studies of ligand-protein interactions. The results of the docking were represented by the binding affinity scores as affinity/ ∆ G (Kcal/mol), and they were further converted to the estimated inhibition constants (Ki). The following formula was used to calculate the Ki parameters from the binding affinity values for each docked pose: K I = exp ∆ G R T (1) where ΔG is the binding affinity or the calculated docking score value, R was Boltzmann gas constant (= 1.987 cal/mol/K), and T was the temperature (= 298 K), respectively.

In comparison to the affinity, the estimated Ki is a much more useful quality because it is more closely related to frequently measured experimental parameters. On the other hand, ligand efficiency (LE) is a crucial informative factor when choosing a lead compound (Kenny, 2019). The following equation is used to determine LE in this situation: L E = − ∆ G / H A (2) where HA is the ligand’s total number of heavy (non-hydrogen) atoms. Potential lead compounds are indicated by compounds with LE 0.3 (Schultes et al., 2010). For additional molecular dynamics simulations, the receptor-ligand binding conformation with the largest negative docking score (binding affinity) was examined.

To assess the stability, Molecular Dynamics (MD) simulations were run and probed the dynamic conformational changes of the selected complexes using Nano Scale Molecular Dynamics (NAMD v2.14) software (Phillips et al., 2005). The receptor-ligand complex obtained from molecular docking simulations was the initial structure used for MD simulations. In CHARMM-GUI, the topologies and parameter files for the complexes (protein-ligand) were generated using Solution Builder (Jo et al., 2008; Jo et al., 2014; Lee et al., 2016; Kim et al., 2017). The system is solvated by adding the TIP3P model to the solvation box, and the counter ions (NaCl 0.15M) were used to neutralize the simulation box. The complex is then subjected to the CHARMM36m force field (Huang et al., 2016) at a constant number of molecules, volume, and temperature (NVT). The steepest descent algorithm was used to complete 20,000 steps of energy minimization. The energy-minimized frameworks were then put to use for simulations at a constant number of molecules, pressure, and temperature (NPT) using Langevin dynamics parameters, the temperature at 310K, and for 10 ns under constant Periodic Boundary conditions to compare trajectories. The MD trajectories were analyzed using the script in the VMD v1.9.3 software (Humphrey et al., 1996), including root mean square deviation (RMSD), root mean square variation (RMSF), solvent accessible surface area (SASA), and radius of gyration (Rg).

The Molecular Mechanics/General-Boltzmann Surface Area (MM-GBSA) method was performed to estimate the binding free energy ( ∆ G b i n d i n g ) of the protein-ligand complex, which is a popular endpoint technique for calculating free energy (Wang et al., 2019). Using MolAlCal (Bai et al., 2021), the binding free energies of the complex in liquid/solution were calculated by Eq. 3: δ G b i n d = δ G c o m p l e x − δ G p r o t e i n + δ G l i g a n d (3) where, δ G b i n d , δ G c o m p l e x , δ G p r o t e i n , and δ G l i g a n d are calculated free binding energy, calculated free binding of the complex, calculated free binding of the protein, and calculated free binding energy of the ligand, respectively.

Drug research time has been significantly shortened in the past few years by the creation of computerized (in silico) modeling techniques to evaluate absorption, distribution, metabolism, excretion, and toxicity (ADMET) attributes. It is easier to exclude compounds with prospective ADMET issues when these attributes can be predicted quickly and accurately. This aids researchers in making decisions about which compounds to generate and test first (Mary et al., 2021)

CoMFA figure plots were created to depict the regions in 3D space around the compound where changes in the steric and electrostatic fields were predicted to increase or decrease activity to speculate about the information content of the derived 3D QSAR model. A full assessment of the obtained contours identifies the crucial physicochemical factors governing the activity and explores the crucial role played by the complex formation in their 3D orientation. The CoMFA contour map was created using Compound 4 as the reference structure. Figure 3A depicts the steric contour map for the CoMFA FFDSEL model. Areas, where the large groups do not encourage activity, are indicated by red contour maps (negative region), while these areas that encourage large groups are indicated by blue contour maps (positive region). Electrostatic interactions are depicted in Figure 3B by contours in maroon and green colors. Green regions only increase activity with negative charges, while maroon regions indicate that positive charges are preferred. To generate higher an-ti-autoimmune disease activity, a bulky or large group is not preferred, according to the region of the steric contour map with the large percentage of red color. For instance, due to fewer bulky groups in that area, the other compounds have higher anti-autoimmune disease activity than compound 4. The reference compound (compound 4) is surrounded by a maroon region on the electrostatic contour map, which suggests that the activity of the electron-withdrawing substituent may be improved. The anti-autoimmune disease activity of certain compounds is found to be increased by an electron-withdrawing group (such as 9, 13, 19, 22, 25, 27, and 29). The involvement of an electron-rich substituent is advantageous at this position, as shown by the green electrostatic contour surrounding the reference compound 4. It suggests that alkoxy groups, such as -OCH3, are preferred in that position over the -CO group because they are more electron-rich and easily accept donated electrons while being incapable of ionization.

The 3D contour maps created using the CoMFA UVEPLS model to suggest a 3D-QSAR model on the target are shown in Figure 4. Red and blue contours stand in for steric interactions, while maroon and green contours stand for electrostatic interactions. The blue contours in this diagram indicate the areas where the involvement of bulkier groups (53% contribution) would contribute to improving biological activity, while the red contours (47% contribution) show the areas where such bulkier groups have the opposite effect and cause biological activity to decrease. The influences of the steric field are shown in Figure 4A. The maroon and green contours on the CoMFA UVE-PLS electrostatic contour maps are depicted in Figure 4B. The blue contours (44% contribution) show the regions where the involvement of electron-rich (electronegative) groups would enhance biochemical activity, while the maroon contours (contribution of 56%) show the areas where the positively charged (or reduced negatively charged) group leads to an increase in biochemical activity.

The thirty-one selected compounds and some FDA-approved drugs (as standard) were studied in silico to select the least binding affinity with the active site of some autoimmune diseases and SARS-coronavirus-2 Mpro, using the AMDock software with AD-Vina algorithm. The docking scores, estimated inhibition constants (K_i), and ligand efficiency (LE) of the binding site of all the receptors used in these studies are summarized in Supplementary Table S1. At first, to choose the best hit-lead compound, the binding affinity of the PTPN22 receptor was used as a criterion to select the least docking score. Based on the findings of the docking screening, eight (International Federation Federation, 2019) compounds with the strongest binding affinity than the standard drugs were decided to portray the binding mode of some autoimmune diseases and SARS-CoV-2 Mpro. The top eight (International Federation Federation, 2019) compounds are presented in Table 4. Whereas compounds (Assiri et al., 2013; Elkhalifa et al., 2018; Zhou et al., 2020; Zhou et al., 2020; Mrid et al., 2022) were found to be consistent as the first eight least binding affinities in the other receptors as presented in Table 4.

The pharmacological properties and toxicity of inhibitors in living beings were calculated using SwissADME (http://www.swissadme.ch/) and ADMETlab 2.0 (https://admetmesh.scbdd.com/) web servers on the top 8 compounds. According to Lipinski’s “rule of five,” (Lipinski et al., 2001) good absorption or permeation is more likely when the molecular weight (MW) is 500 Da or less, there are five or fewer hydrogen bond donors present, LogP is five, and there are ten or more hydrogen bond acceptors present. Two additional pertinent descriptors were found by Veber et al. (2002) to be the number of rotatable bonds (NBR) 10 and the polar surface area (PSA) 140 Å2. The results showed that all inhibitors satisfied Lipinski and Veber Rule. All the compounds also show moderately soluble lipid and water solubility except compounds (18 and 23) (Table 5). According to the literature, drugs typically have seven rotatable bonds, whereas toxins only have three (Khanna and Ranganathan, 2009). We discovered that the number of rotatable bonds in the compounds (Rieck et al., 2007; Elkhalifa et al., 2018; Zhou et al., 2020; Zhou et al., 2020) ranges from 0 to 7 (Table 5). The BOILED-Egg is an easy-to-use model for foretelling the biodistribution of organic compounds (Daina and Zoete, 2016). Two compounds were located inside the egg, which represents an appropriate physicochemical space for biodistribution, when we mapped WLOGP and TPSA of the virtual screening hits to the BOILED egg (Figure 5). Lead compounds found in egg white, which indicates human intestinal absorption (HIA) without blood-brain barrier (BBB) permeation, would be favored for faster drug development related to the treatment of autoimmune diseases. While 4 compounds are still close to egg white and would have better bioavailability profiles during a drug development phase, 6 compounds are in the gray area and fall into this category. Compounds inside the egg white are molecules 1 and 7, that is compounds 4 and 27, respectively.

CYP1A2, CYP2C19, CYP2C9, CYP2D6, and CYP3A4 are five major cytochrome P450 isoforms that have a significant role in drug metabolic activities and abolition. As a result, these isozymes play an important role in controlling drug reactions, which in turn can determine effectiveness and undesirable effects. Data on the ability of top virtual screening hits to inhibit significant CYP isozymes can be found in Supplementary Table S2. We discovered that compound 4, which would inhibit four isozymes, is likely to show more response to therapy. The remaining compounds may display two or three isozymes. We then calculated the toxicity statuses of the chosen compounds using the ADMETlab webserver and OSIRIS Property Explorer (Dong et al., 2018). Eight of the chosen compounds could have high toxicities, as shown in Supplementary Table S3. It should be noted that compound 4 is expected to have high reproductive effectiveness. Compound 10 is predicted to have high irritants. Compounds with high hepatotoxicity are compounds (4, 10, 16. 18, 23, and 27). Compounds with high ames mutagenicity are compounds 14 and 16. Compound 4 appears to be a better lead without significant toxicity, on the other hand.

One compound stood out when we looked at the ADMET profiles of the top eight hits from the target identification (Figure 3). These molecules are predicted not to cross the BBB, have rheological properties suitable for absorption and bioavailability, and appear to carry no or lower toxicity risks. The radar charts show that the properties of these biomolecules that favor absorption and bioavailability are lipophilicity, size, polarity, insolubility, instauration, and flexibility (Figure 6). It’s noteworthy that com-pound 4 occupies the entire physicochemical region for absorption and bioavailability. It is also tempting to mention that, out of the top eight hits, compound 4 poses only reproductive effective as toxicity risks (Supplementary Table S3). Compound 4’s analyses, as seen in Table 3; Supplementary Table S2, revealed no offenses of these guidelines, indicating that they would exhibit well-behaved absorption or permeation.

Under the results of toxicity and ADME analyses of the top eight compounds, compound 4 with favored oral bioavailability attributes was selected for docking visualization/interactions with the receptors. Figure 7 displays numerous interactions of compound 4 with the least docking score with PTPN22 (4J51.pdb), diabetes type 1 (1JK8.pdb), and rheumatoid arthritis (2FSE.pdb) also interact with the critical residues of the SARS-CoV-2 main protease (6LU7.pdb). Molecular docking calculations confirmed that compound 4 can occupy the catalytic sites of all the receptors and produce a net of hydrogen (Figures 7A–D). From the results, it has been observed that compound 4 formed two conventional hydrogen bonds and one carbon-hydrogen bonds interactions with the PTPN22 receptor was shown in Figure 7A. Docking analysis of the PTPN22 receptor with compound 4 allowed us to recognize certain residues, viz. ARG266, TYR44, PRO45, LYS42, LYS39, and SER35, within the PTPN22 receptor binding pocket, which is crucial to ligand binding affinity. The modes of interaction between compound 4 and the active site diabetes type 1 (PDB ID: 1JK8) are shown in Figure 7B, which shows four conventional hydrogen bonds and three carbon-hydrogen bond interactions between the ligand and amino acids GLY20, VAL34, TYR22, LYS147, SER19, GLU134, and SER136 respectively. Thus, the one pi-pi T-shaped bonds with the amino acids PHE137. In addition, two pi-alkyl-type interactions are observed between the compound and the amino acids VAL116 and TYR33, respectively. However, compound 4 does not show any carbon-hydrogen-like bond with rheumatoid arthritis (PDB ID: 2FSE). The two-dimensional visualization indicates that compound 4 inter-acted with the amino acid TYR13 via conventional hydrogen bonds. Compound 4 was also bound to the ILE82 site by an alkyl bond. Additional van der Waals bonds also occurred at the ASP66, PHE12, GLU11, HIS143, PHE154, SER113, THR74, ASN15, LEU70, and LEU14 sites (Figure 7C). These results provide evidence of the critical role of the amino acids tyrosine (TYR), lysine (LYS), and serine (SER) for the stability of the compound 4 in the active site of autoimmune receptors. Therefore, these amino acid residues should be taken into account to improve the biological inhibitory activity of compound 4 analogs against autoimmune disorders.

The two crucial peptides in the catalytic site of the SARS-CoV-2 target site are CYS145 and HIS41. Studies have shown that the catalytic dyad formed by CYS145 and HIS41 increases the reactivity of the nucleophile by acting as a base and a nucleophile, respectively. The SARS-CoV-2 Mpro is then released after the inhibitor attacks the dyad to create an alternate complex. Compound 4 interacts with CYS145 and HIS41 in our study, indicating that these substances may have inhibitory activity against SARS-coronavirus-2 Mpro. Asymmetric aromatic disulfides may fight with the substrate in the SARS-CoV-2 Mpro cysteine protease as an inhibitor (competitive inhibitor). We think compound 4 has the potential to be an enzyme inhibitor and, as a result, could be a drug candidate for both autoimmune disease and the SARS-CoV-2 virus. An in vitro study against autoimmune disorder and SARS-CoV-2 Mpro is advised to verify this hypothesis.

In the case of the PTPN22 enzyme, the control compound (Sulfasalazine) has the highest docking affinity than the rest of the standard drugs. The sulfasalazine-4J51 complex was stabilized by the formation of conventional hydrogen bonds with SER35 (4.22 Å), SER271 (4.17 Å), and THR46 (3.59 Å) of protein with PTPN22. The complex was also stabilized by hydrophobic interactions, including pi-pi T-shaped interaction with TYR44, Amide-pi stacked bond with TYR66, and pi-alkyl bonds (with PRO270, PRO45, LYS39, and LYS42), as shown in Supplementary Figure S1A. The amino acid residue for compound 4 and sulfasalazine complexes revealed from docking simulations results shows that (GLY20, LYS147, VAL34, TYR22, SER35, SER271, and THR46) are essential for the binding of ligands to the enzyme of PTPN22.

Furthermore, for the type 1 diabetes receptor, among the two control drugs, sotagliflozin has the highest docked score of −7.3 kcal/mol. Sotagliflozin has no single hydrogen bond but was stabilized with the formation of hydrophobic interactions, including Alkyl interactions with TYR33, VAL116, TYR16, pi-Alkyl interaction with LEU36, and pi-pi- T-shaped interaction with PHE137. The most common amino acid residue involved in the interactions of compound 4 and sotagliflozin drugs is TYR33, LEU36, PHE137, VAL116, and TYR16. The docking interactions of both compounds revealed that GLY20, VAL34, TYR55, SER19, TYR33, LEU36, PHE137, VAL116, and TYR16 are essential for the binding of ligand to the enzyme of type 1 diabetes (PDB id: 1JK8).

Here, for rheumatoid arthritis protein, sulfasalazine (standard drug) has the highest docked score of −7.7 kcal/mol. In the complex of PDB id 2FSE with sulfasalazine, three conventional Hydrogen bonds were formed by sulfasalazine with ARG44 (5.08 Å), TYR150 (6.92 Å), and ASP29 (3.55 Å) of 2FSE (Supplementary Figure S1C). Other than conventional H-bonds, the hydrophobic interactions such as pi-Alkyl bonds (with LEU45, ILE31, ALA52), pi-pi stacked bonds (with PHE28), pi-pi T-shaped bonds (with PHE28), and electrostatic interaction (pi-cation) with ARG44 were involved in the stabilization of the complex. The amino acid residue for compound 4 and sulfasalazine complexes against rheumatoid arthritis revealed from docking simulations results shows that TYR13, ARG44, TYR150, ASP29, and ARG44 are essential for the binding of ligands to the enzyme (PDB id 2FSE).

Lastly, for the SARS-CoV-2 protein (PDB id: 6LU7), the control drug hydroxychloroquine has a binding affinity of −6.3 kcal/mol. Hydroxychloroquine formed two conventional hydrogen bonds with THR190 (4.22 and 4.38 Å) of SARS-CoV-2 protein in the hydroxychloroquine-6LU7 complex (Supplementary Figure S1D). Further, this complex was also stabilized by one carbon-hydrogen bond with GLN189 (4.43 Å), one pi-cation ionic interaction with HIS41, and hydrophobic (pi-Alkyl and Alkyl) interactions with CYS145 and MET49. The common amino acids involved in the interactions with compound 4 and the standard drug hydroxychloroquine are HIS41, CYS145, and MET49. Therefore, the interactions revealed that amino acids ASN142, MET165, THR190, GLN189, HIS41, CYS145, and MET49 are essential for the binding of ligands to the enzyme (PDB id 6LU7). Dhankhar et al. (2020), Kumar et al. (2021), and Dhankhar et al. (2021), examined research on ligand access channels in SARS-CoV-2 receptor and found that amino acid residues Met49, Cys145, and Gln189 in particular may play a significant function in ligand binding interactions. Our data corroborate that hypothesis and identify Met49, Cys145, and Gln189 as an amino acid residues implicated in the binding of SARS-CoV-2 Mpro.

The energetic motion and level of stability of the complexes were examined and understood using the MD simulation. MD simulation for all the complexes of compound 4 is performed for 10 ns (Figure 8). RMSD of alpha carbon atoms, the radius of gyration (RG), RMSF, and solvent accessible surface area (SASA) of the complexes are investigated. Figure 8A shows how compound 4 complexes are stable before 10 ns and how the RMSD of the protein backbone and compound 4 are related. Interestingly, none of the complexes showed RMSD values greater than 2, which supports the strict specificity of the most active complexes. Throughout the remainder of the simulation, these complexes kept their RMSD profile stable. When compound 4 is tested against the PTPN22 protein, it exhibits the lowest RMSD value compared with the remaining complexes, indicating that it is more stable and stays in the protein’s pouch (Table 6). To determine the stability of receptors through MD simulation, the RG of complexes with compound 4 was also estimated (Figure 8B). In comparison to other receptors, the rheumatoid arthritis receptor generally exhibits less variation in RG values. The average RG values of receptors with compound 4 are 20.19, 18.82, 18.70, and 22.52 Å, respectively. The complexes’ stability is indicated by this small value (Table 6). The RMSF study provided insightful data on the structural fluctuations of various protein regions. Increased variations in the residues indicate that the protein is less stable. The consistency of compound 4 with the sample protein is indicated by the fact that the RMSF of the complexes, which is depicted in Figure 7C, remained below 0.5 for most of the receptors’ amino acids. Nevertheless, a few variations are seen at the terminal, which might be a result of these residues’ high plank position in the rheumatoid arthritis receptor. The SARS-CoV-2 Mpro receptor exhibits some variations as well, which could be brought on by compound 4’s dynamic characteristics in the bonding zone. To further examine the trajectories of compound 4, the average RMSF of com-pound 4 was also calculated. This calculation revealed some variation, indicating a kinetic shift from their start point. Compound 4 must therefore be categorized as a drug. Additionally, following MD simulation, the lowest RMSF value of compound 4 against type 1 diabetes (PDB ID: 1JK8) was also discovered. Calculating the changes in SASA allowed for further confirmation of the stability. The SASA of the receptors with compound 4 are shown in Figure 8D. The SASA of diabetes type 1 with compound 4 complex is similar to rheumatoid arthritis with compound 4. Likewise, PTPN22 and SARS-CoV-2 with compound 4 are similar, further confirming the stability of compound 4 with all the protein crystal structures. The average SASA values for all the receptors-compound 4 are presented in Table 6.

The MM/GBSA approach is used to evaluate the calculated free binding energy of compound 4 from each receptor. The free binding energy for each complex was calculated using the molecular dynamics trajectory from the previous 2 ns (100) frames. The effects of other non-bonded interaction energies were predicted along with the calculated free binding energy of each compound 4-receptor complex. With diabetes type 1 protein, compound 4 has a calculated free binding energy of −45.1538 ± 0.3823 kcal/mol, whereas, with rheumatoid arthritis protein, compound 4 has a calculated binding energy of −42.0993 ± 0.4699 kcal/mol. Interactions like internal energy, van der Waals energy, electrostatic, and solvation energy is governed. Across all types of interactions, the van der Waals and the sum of electrostatic and the calculated free binding energy were mainly influenced by solvation energies. The internal energy, on the other hand, made the smallest contribution to the final calculated free binding energies. Moreover, the internal energy interaction values of compound 4 and SARS-CoV-2 Mpro protein complexes demonstrated high contribution after van der Waals energy (Table 6). Figure 9A shows that compound 4 in the catalytic pocket of the PTPN22 protein has undergone a significant geometric change in the pose after simulation (10 ns), and Figures 9B–D shows that compound 4 in the catalytic pocket of DT1 and RA has also undergone a significant elongated change after the simulation (10 ns) (curved to straight). Effectively receptor—ligand acquisition and interactions with residues result from these conformational changes, which increase stability and binding energy (shown in Table 6). Therefore, MM-GBSA calculations from MD simulation trajectories were well supported by the binding affinity found in the docking results. In addition, the last frame (10 ns) of MD simulation showed the positional change of the compound 4 and protein targets as compared to docking (shown in Figure 7), indicating the best binding position for effectively fitting in the catalytic pocket of the protein (see Figure 9).

To further analyze the binding between compound 4 and the receptor complexes, the contactFreq.tcl module in VMD (cut-off of 4 Å) was used. The proportion of compound 4 contact frequency with the receptors’ binding amino acid residues was reported, and the results indicate that certain residues are engaged in compound 4’s ongoing interactions with those residues (Supplementary Table S4).