In-silico discovery of common molecular signatures for which SARS-CoV-2 infections and lung diseases stimulate each other, and drug repurposing.

The SARS-CoV-2 virus is responsible for COVID-19 which was the first outbreak in Wuhan city, Hubei province, China, in December 2019 [ 1 ]. Its outbreak became a terrible form rapidly whole over the world and the WHO officially announced it to be a global epidemic on March 11, 2020 [ 1 ]. Although COVID-19 has hurt almost all countries, notably the United States, India, France, Germany, Brazil, South Korea, and Japan are the top seven countries affected by SARS-CoV-2 ( https://www.worldometers.info/coronavirus/ ). According to the WHO report, until December 31, 2023, over 6.9 million people out of 700 million SARS-CoV-2 infected people have died. Though infection rates are gradually decreasing worldwide due to the impact of vaccination, however, some people are yet infecting [ 2 ]. It may have happened due to the unstable RNA pattern of SARS-CoV-2 and the weak immunity of the patients. Not every patient with SARS-CoV-2 infection suffers in the same way. Some patients become more vulnerable, who are already suffering from one or more comorbidities like cardiovascular diseases [ 3 ], diabetes [ 4 ], hypertension [ 3 ], and different lung diseases including chronic obstructive pulmonary disease (COPD) [ 5 ], idiopathic pulmonary fibrosis (IPF) [ 6 ], interstitial lung disease (ILD) [ 7 ], asthma [ 8 ], tuberculosis [ 9 ], cystic fibrosis [ 4 ], pneumonia [ 1 , 10 ], emphysema [ 11 ], and bronchitis [ 12 ]. Chronic obstructive pulmonary disease (COPD) is a chronic inflammatory lung disorder that encompasses chronic bronchitis and emphysema, and it is characterized by restricted airflow with symptoms of breathing problems and cough with mucus due to abnormalities in the airways or air sacs of the lungs [ 13 ]. It is the third leading cause of death globally [ 14 ]. Idiopathic pulmonary fibrosis (IPF) is a long-term, progressive lung disorder characterized by lung scarring or fibrosis that leads to respiratory failure [ 15 ]. Patients with IPF survive only about 3–5 years after diagnosis with symptoms of dry cough and shortness of breath [ 16 ]. Respiratory failure is responsible for death related to IPF [ 16 ]. Interstitial lung disease (ILD) is a lung disorder that can stimulate both the vulnerability and severity of COVID-19 [ 7 ]. Tuberculosis (TB) constitutes the predominant cause of mortality associated with respiratory infections. Furthermore, TB significantly augments the susceptibility to COVID-19 while simultaneously exacerbating the severity of the disease [ 9 ]. People with cystic fibrosis (CF), a chronic lung disease involving mucus blockage and persistent airway inflammation, are highly vulnerable to COVID-19 due to the increased risk of severe viral respiratory infections [ 4 ]. Moreover, pneumonia is a type of lung disease characterized by inflammation in the tiny air sacs within the human lungs, leading to the accumulation of fluid and resulting in breathing difficulties [ 17 ]. The S-protein of SARS-CoV-2 has a higher interaction with ACE2 (angiotensin-converting enzyme 2); however, a significant amount of ACE2 is found in lung disease patients [ 13 ]. Therefore, in this study, we have considered chronic lung disease patients as a high-risk group for COVID-19 complications [ 7 , 13 ]. Several studies explored SARS-CoV-2 infections causing key genes (KGs). Some of these studies detected shared key-genes (sKGs) to disclose common pathogenetic processes of SARS-CoV-2 infections with one or two lung diseases including COPD [ 13 ], IPF [ 6 ], COPD and IPF [ 13 ], ILD [ 7 ], asthma [ 18 ], tuberculosis [ 19 ], cystic fibrosis [ 20 ], pneumonia [ 21 ], emphysema [ 13 ], and bronchitis [ 13 ]. Few of these studies recommended sKGs-guided common drug molecules in which molecules (curcumin, triclosan, tamoxifen, deguelin) were recommended for the treatment of SARS-CoV-2 infections with COPD [ 13 ], molecules (tegobuvir, nilotinib, digoxin, proscillaridin, simeprevir, sorafenib, torin 2, rapamycin, vancomycin and hesperidin) with IPF [ 6 ], molecules (suloctidil, estradiol, prenylamine, clioquinol) with asthma [ 18 ], molecules (rituximab, bevacizumab, bosentan, sitaxentan, and macitentana) with tuberculosis [ 19 ], molecules (imiquimod and raloxifene) with cystic fibrosis [ 20 ]. However, so far, there is no study that explored sKGs/sDEGs to disclose common pathogenetic mechanisms and associated drug molecules for SARS-CoV-2 infections and different lung diseases. Therefore, the main objective of this study is to explore potential sKGs/sDEGs to reveal the common pathogenetic mechanisms of SARS-CoV-2 infection and different lung diseases to adopt a common treatment plan. The workflow of this research is shown in Fig 1 .

Three RNA-Seq count datasets ( GSE147507 , GSE150392 , and GSE152075 ) were analyzed by using three methods DESeq2, edgeR, and LIMMA (voom), and detected a total of 888 shared DEGs between COVID-19 and control samples as [ ( ⋂ i = 1 3 C i D E S e q 2 ) ∪ ( ⋂ i = 1 3 C i e d g e R ) ∪ ( ⋂ i = 1 3 C i L I M M A ( v o o m ) ) ] = 888 Then, we analyzed GSE64913 , GSE34608 , GSE107846 , GSE42830 , GSE1122 , GSE22148 , GSE53845 , GSE40839 , and GSE100281 datasets for Asthma, Tuberculosis, Cystic Fibrosis, Pneumonia, Emphysema, Bronchitis, IPF, ILD, and COPD, respectively using three methods KW, SAM and LIMMA and combinedly identified total 3698 DEGs for these lung diseases as [ ⋂ j = 1 9 ( L j K W ∪ L j S A M ∪ L j L I M M A ) ] = 3698 Then we found 267 shared DEGs (sDEGs) between COVID-19 and different lung diseases using Eq 1 ( Fig 2A ) . (A) Venn diagram of DEGs for COVID-19 and lung diseases to display shared DEGs (sDEGs). The blue and green circles represent DEGs from three datasets of COVID-19 and nine datasets of lung diseases, respectively. (B) Protein-protein interaction network of 267 sDEGs, where the large size maroon color octagonal node indicates hub-sDEGs. (C) Venn diagram among hub-sDEGs, reviewed DEGs of COVID-19, and reviewed DEGs of lung diseases. The PPI network of 267 sDEGs was constructed which contained 267 nodes and 291 edges displayed in Fig 2B to detect the biomarker genes named hub-sDEGs. Based on the degree of importance, we selected the top-ranked 11 genes: STAT1 , TLR4 , CXCL10 , CCL2 , JUN , DDX58 , IRF7 , ICAM1 , MX2 , IRF9 and ISG15 as the hub-sDEGs and used for further analysis. Among them, 5 genes ( CCL2 , CXCL10 , ICAM1 , JUN , TLR4 ) are common with text mining of COVID-19 and lung disease presented through the Venn diagram in Fig 2C . To investigate the different disease risk factors of SARS-CoV-2 infections from the genetic viewpoint, the interaction network analyses of hub-sDEGs with different diseases were performed which revealed that 6 hub-sDEGs ( ICAM1 , MX2 , CCL2 , IRF7 , JUN , STAT1 ) out of 11 are significantly associated with different diseases including asthma, atherosclerosis, bronchiectasis, cardiovascular diseases, brain ischemia, diabetes mellitus, obesity, schizophrenia, urticaria, melanoma, pneumonia, pulmonary fibrosis, tuberculosis, lung injury, lung neoplasms, rheumatoid arthritis, liver cirrhosis, hypertension, etc. ( Fig 3 and S3 Table ) . Here the green circle node indicates significant disease risk and the maroon-colored octagonal node indicates hub-sDEGs. GO functions are classified into three subsections: biological processes (BPs), molecular functions (MFs), and cellular components (CCs). The few top common enriched GO functions for each of the three subsections (BPs, MFs, and CCs) based on two online databases (DAVID and Enrichr) are presented in Table 1 . According to the BP-terms, almost all the hub-sDEGs were mainly enriched with the defense response to virus, positive regulation of interferon-alpha/beta production, positive regulation of transcription by RNA polymerase II, positive regulation of DNA-templated transcription and cellular response to type II interferon, etc. Among the enriched CC-terms, endosome membrane, cytoplasmic vesicle membrane, nucleus, euchromatin, and intracellular membrane-bounded organelle are the top enriched GO functions. In the case of MFs, hub-sDEGs are significantly enriched in chemokine activity, transcription regulatory region sequence-specific DNA binding, RNA polymerase II core promoter proximal region sequence-specific DNA binding, double-stranded DNA binding, etc. Pathway enrichment analysis of the hub-sDEGs based on four databases (KEGG, WikiPathways, Reactome, and BioCarta) identified some crucial pathways ( S4 Table ) . Top-ranked four pathways (Coronavirus disease, Interferon signaling pathway, Immune responses, and Measles virus infection) supported by at least two databases are shown in Table 2 . To identify key transcriptional and post-transcriptional regulatory factors of hub-sDEGs, we constructed the interaction network among miRNAs, TFs, and hub-sDEGs as depicted in Fig 4 . From this network, we selected top-ranked six significant TFs proteins (MYC, SOX2, CEBPA, NANOG, RELA, and MSX1) and six significant miRNAs (hsa-miR-16-5p, hsa-miR-129-2-3p, hsa-miR-21-3p, hsa-miR-27a-5p, hsa-miR-1-3p, hsa-miR-155-5p) as the key transcriptional and post-transcriptional regulatory factors of hub-sDEGs, based on the degree score. The key miRNAs are indicated by larger round square nodes. To explore candidate repurposable drugs as a therapy for SARS-CoV-2 infections with one or more lung diseases, we choose our suggested 11 hub-sDEGs based proteins (STAT1, TLR4, CXCL10, CCL2, JUN, DDX58, IRF7, ICAM1, MX2, IRF9, and ISG15) and its regulatory 6 key TFs proteins (MYC, SOX2, CEBPA, NANOG, RELA, and MSX1) as the m = 17 target proteins and 184 meta-drug agents (ligands). The 3D structures of 11 receptor proteins (CCL2, DDX58, ICAM1, ISG15, JUN, MYC, NANOG, RELA, SOX2, STAT1, TLR4) were downloaded from the PDB with source code (1dok, 2lwd, 2oz4, 3phx, 1jun, 1nkp, 2kt0, 4kv1, 2le4, 3wwt, 3fxi) and remaining 6 receptor proteins (CEBPA, CXCL10, IRF7, IRF9, MSX1, MX2) were obtained from AlphaFold source using Uniprot IDs P49715 , P02778 , Q92985 , Q00978 , P28360 , P20592 . Then molecular docking analyses were performed between m = 17 drug target proteins and n = 184 drug agents to obtain the binding affinity score (kcal/mol) for each target protein with each drug agent ( S5 Table ) . Fig 5A represents the ordered top-ranked 30 drugs binding affinity score matrix out of 184 drugs. We observed that the first four top-ordered compounds (Nilotinib, Entrectinib, Imatinib, and SCH-772984) produce highly significant binding affinity scores (≤ -7) kcal/mol with at least 14 receptor proteins out of 17 and their average binding affinity score is less than or equal to −8.0 kcal/mol. The remaining six top-ordered compounds (Icotinib, Regorafenib, Zanubrutinib, Sorafenib, Ibrutinib, and Dabrafenib) also produced significant binding with at least 12 receptors out of 17, and their average binding affinity score (≤ −7.4) kcal/mol. Hence, we considered these top-ranked 10 compounds (Nilotinib, Entrectinib, Imatinib, SCH-772984, Icotinib, Regorafenib, Zanubrutinib, Sorafenib, Ibrutinib, and Dabrafenib) as the most viable candidate drugs to treat COVID-19 as well as patients with lung diseases and highlighted them in Fig 5A . To investigate the resistance performance of the considered drug molecules, and compare to the already published molecules against the state-of-the-arts alternatives top-ranked independent receptors published by other studies, we reviewed 27 published articles related to COVID-19 ( S6 Table ) and 27 published articles associated with lung diseases ( S7 Table ). Then, we selected eight top-ranked target proteins that are commonly reported in at least two articles in both SARS-CoV-2 and lung disease-related literature and highlighted them in the 5th column in S6 and S7 Tables . These 8 target proteins were considered as the top-ranked independent meta-receptors to examine the resistivity of the considered drugs compared to the selected top-ranked published drugs by molecular docking analysis. To examine the resistivity of the considered 10 candidate drugs against the top-ranked eight independent meta receptors, we downloaded the 3D structures of the seven (ICAM1, JUN, MMP1, CXCL8, CXCL1, VEGFA, and IL6) independent meta receptor proteins from the PDB with source codes 2oz4, 1jun, 1ayk, 1o7b, 1msh, 1kat, and 1il6, respectively and the protein SOSC3 downloaded from the AlphaFold database using Uniprot ID O14543 . Then, we performed molecular docking analysis of published drugs with the top-ranked 8 independent meta-receptor proteins, and their docking results are given in Fig 5B and S8 Table . We observed that among the considered drugs, 9 drugs (Nilotinib, Entrectinib, Imatinib, SCH-772984, Icotinib, Regorafenib, Sorafenib, Ibrutinib, and Dabrafenib) belong to the top-ranked ten drugs against the independent receptors also. Hence, we selected these nine drugs as the potential drugs for the therapy of SARS-CoV-2 infections and lung diseases. To verify the significance of the protein-ligand interaction for our proposed target proteins and the selected ligands (drugs), we also took into account an alternative docking tool known as the GLIDE module of Schrödinger and calculated binding scores (Glide g-score, Docking score, and Glide e-model) ( Fig 5C and S9 – S11 Tables ). (A) Image of BAS of the top-ranked 30 candidate drugs out of 184 meta-drug agents in the X-axis with the proposed top-ranked 17 drug target proteins in the Y-axis. (B) Image of BAS of the top-ranked 30 candidates published drugs in the X-axis against the published top-ranked 8 independent meta receptors in the Y-axis. (C) Image of Glide g-score of the top-ranked selected nine candidate drugs in the X-axis with the proposed top-ranked 17 drug target proteins in the Y-axis. Selected nine compounds from molecular docking analysis were further optimized through density functional theory (DFT). For that, Frontier molecular orbitals (FMO) properties (such as HOMO and HUMO) of nine compounds were computed to determine the significance of charge-transfer interactions at the protein binding site ( Table 3 ). Among the selected candidate drugs—Nilotinib, Entrectinib, and Imatinib were the top-ranked three candidate drugs, and we choose RELA for MDS, since RELA vs . Nilotinib, RELA vs . Entrectinib, and RELA vs . Imatinib complexes produce highest binding affinity scores compare to the other complexes ( Fig 5 ). Therefore, these top three drug-target complexes were selected for their stability analysis with the molecular dynamics (MD) simulation for a duration of 100 nanoseconds. To perform MD simulation using YASARA, the initial energy minimization of each simulation system included 33769, 33564, and 44604 atoms for the RELA vs . Nilotinib, RELA vs . Entrectinib, and RELA vs . Imatinib complexes, and the average RMSDs for these systems were 1.42 Å, 1.37 Å, and 1.89 Å, respectively. These three complexes were almost stable throughout the whole simulation and all of the systems projected RMSD between 0.36 Å and 2.95 Å, though in the case of RELA vs . Nilotinib complex, a steady movement was found from 38 ns to 43 ns ( Fig 6A ). To examine the stability of these three complexes by another software, we considered the Desmond module of Schrödinger. The conformational stability of biological molecules were investigated using the Root Mean Square Deviation (RMSD) [ 43 ]. The protein backbone RMSDs were used to assess the deviation of the RELA protein bound with proposed inhibitors. The RMSD values of RELA for each frame computed through MD simulation production bound with all three drug molecules were plotted against the time scale, as shown in Fig 6B , which shows that the three complexes RELA vs . Entrectinib, RELA vs . Imatinib, and RELA vs . Nilotinib average RMSDs were 1.85 Å, 1.89 Å, and 2.47 Å, respectively. These results indicated that the complexes of Entrectinib and Imatinib with RELA are almost stable throughout the whole simulation and all of the systems projected an RMSD between 1.0 Å and 2.75 Å. However, the RELA vs . Nilotinib complex shows a slight fluctuation than the other two complexes between 0 to 35 ns, after reaching at 40ns, it seems quite stable and no complex has the RMSD greater than 3.0 Å. Thus, in the case of RMSD, it is clearly observed from both software that none of those three complexes has more fluctuation and they formed stable interactions. (A) The RMSD analysis results with YASARA. (B) The RMSD analysis results with the Desmond module of Schrödinger. To investigate the stability of the top-ranked three complexes, we also determined the MM-PBSA and MM-GBSA binding energies for these complexes ( Fig 7 ) . From Fig 7A , we observed that the average binding energies obtained from the RELA vs . Nilotinib, RELA vs . Entrectinib, and RELA vs . Imatinib complexes were 132.44 kJ/mol, 180.87 kJ/mol, and 79.50 kJ/mol, respectively. These positively average binding energies with YASARA indicate the significant stability of these three complexes [ 44 ]. To examine the stability of those complexes with another software, we calculated MM-GBSA binding free energy (ΔG) for the MD simulation with the Desmond module of Schrödinger ( Fig 7B ) . Fig 7B indicates that the average ΔG of RELA vs . Entrectinib, RELA vs. Imatinib, and RELA vs . Nilotinib are -146.24 kJ/mol, -148.59 kJ/mol and -94.84 kJ/mol, respectively. Also, the results ( S1 Fig ) show that the overall binding energy for RELA vs . Imatinib is -35.51 kcal/mol after adding the solvation term of 27.65 kcal/mol. Moreover, RELA vs . Entrectinib and RELA vs . Nilotinib showed the average VDW energy as -39.79 kcal/mol and -36.47 kcal/mol respectively, whereas the electrostatic energy showed -65.83 and -14.43 kcal/mol, respectively. The ΔG Total and solvation energy for RELA vs . Entrectinib is -35 and 70.67 kcal/mol. These negative average binding energies with the Desmond module indicate the significant stability of those three complexes as before [ 45 ]. (A) The MM-PBSA binding energies for the MD simulation with YASARA (B) The MM-GBSA binding energies for the MD simulation with the Desmond module of Schrödinger. We also investigated the drug-target binding positions in the proposed top-ranked three complexes docked by AutoDock Vina in Fig 8 . The 2D schematic diagram of the receptor-ligand interactions is given in the fourth column. The third column displays the 3D view of the protein-ligand complex, and the interactive key amino acids in the docked complex are mentioned in the fifth column.

In this study, we identified the top-ranked 11 genes ( STAT1 , TLR4 , CXCL10 , CCL2 , JUN , DDX58 , IRF7 , ICAM1 , MX2 , IRF9 , and ISG15 ) as the hub of the shared differentially expressed genes (hub-sDEGs) highlighting their molecular mechanisms for which SARS-CoV-2 infections and different lung diseases stimulate each other. The GO enrichment analysis with hub-sDEGs revealed that these hub genes are involved in cellular response to interferon-gamma, defense response to virus, positive regulation of interferon-alpha or beta production, etc., and pathways were enriched in coronavirus disease, TNF signaling pathway, interferon signaling pathway, etc., that are associated with both COVID-19 and lung diseases. Some TFs proteins (MYC, SOX2, CEBPA, NANOG, RELA, and MSX1) and miRNAs (hsa-miR-16-5p, hsa-miR-129-2-3p, hsa-miR-21-3p, hsa-miR-27a-5p, hsa-miR-1-3p, and hsa-miR-155-5p) were also detected as the key transcriptional and post-transcriptional regulatory factors of key/hub-sDEGs. Then, hub-sDEGs guided top-ranked three repurposable drug molecules (Entrectinib, Imatinib, and Nilotinib) were detected by molecular docking analysis with AutoDock Vina and GLIDE module of Schrödinger. Additionally, molecular dynamic simulations confirmed the stability of the top-ranked three drugs (Entrectinib, Imatinib, and Nilotinib) in complex with the RELA protein over a 100 ns period. Therefore, the findings of this study would be useful resources for the diagnosis and therapies of COVID-19 patients who are also suffering from chronic lung diseases.

In this study, we analyzed three RNA-Seq profile datasets for SARS-CoV-2 infections and nine microarray gene expression datasets for nine types of lung diseases (COPD, IPF, ILD, Asthma, Tuberculosis, Cystic Fibrosis, Pneumonia, Emphysema, and Bronchitis) that were collected from the Gene Expression Omnibus (GEO) platform of National Center for Biotechnology Information (NCBI) database [ 22 ]. The detail information of these datasets is given in S1 Table . We considered three methods DESeq2, edgeR and LIMMA (voom) for the identification of DEGs between SARS-CoV-2 infections and control groups separately from RNA-Seq profile datasets [ 23 , 24 ]. To identify DEGs between lung disease and control samples based on microarray gene expression profile datasets, separately we utilized KW, SAM, and LIMMA approaches [ 25 ]. R-software (Version 4.0.5) has been used to implement these approaches for data analysis. Let C i D E S e q 2 , C i e d g e R and C i L I M M A ( v o o m ) are three DEGs sets computed from i th RNA-Seq count dataset by DESeq2, edgeR and LIMMA (voom), respectively, for COVID-19 ( i = 1,2,3). Again, let L j K W , L j S A M and L j L I M M A are three DEGs sets computed from microarray gene expression profiles by KW, SAM and LIMMA, respectively, for j th lung disease ( j = 1,2,…,9). Then we computed shared DEGs (sDEGs) between COVID-19 and different lung diseases as follows, s D E G s = [ ( ⋂ i = 1 3 C i D E S e q 2 ) ∪ ( ⋂ i = 1 3 C i e d g e R ) ∪ ( ⋂ i = 1 3 C i L I M M A ( v o o m ) ) ] ∩ [ ⋂ j = 1 9 ( L j K W ∪ L j S A M ∪ L j L I M M A ) ] , (1) Since different methods (utilized in this study) identify DEGs based on different assumptions on expression data; however, a method fails to detect some potential DEGs when its assumption is not fully satisfied by the dataset [ 26 ]. To construct protein-protein interaction (PPI) network, we inserted sDEGs in the ‘Search Tool for the Retrieval of Interacting Genes (STRING)’ database (version 11.5). Then, we used Cytoscape (version 3.10.0) software with a confidence score ≥ 0.90 to analyze and visualize the PPI network. We implemented a Cytoscape plugin Network Analyzer (version 4.4.8), to investigate the interactions among sDEGs. Then, we used the cytoHubba (version 0.1) plugin in Cytoscape to identify significant nodes or hub-sDEGs based on the degree scores [ 25 , 27 ]. To investigate the disease risk factors of SARS-CoV-2 infections through the hub-sDEGs, we performed ‘gene-disease’ interaction network analysis by using the web-tool ‘NetworkAnalyst (version 3.0)’ with the database ‘DisGeNET’ [ 28 ]. A disease has been considered as significantly associated with hub-sDEGs if adjusted p -value < 0.05. The GO functional and pathway enrichment analysis of hub-sDEGs were performed by using ‘Enrichr’ and ‘DAVID’ databases [ 6 , 24 ]. Then, we extracted common GO terms and pathways from both databases. The adjusted p -value < 0.05 was considered as the threshold value for this analysis. To obtain key transcriptional and post-transcriptional regulators of hub-sDEGs, we performed network analysis of hub-sDEGs with TFs and miRNAs from TF2DNA and TarBase (v8.0) databases, respectively, and their interaction networks were constructed by using STRING (version 11.5) database and visualized in Cytoscape (version 3.10.0) software [ 24 ]. The significant key TFs and miRNAs were identified via the ‘CytoHubba’ plugin in Cytoscape based on the highest degree scores. To explore hub-sDEGs mediated receptor-proteins guided few potential repurposable drug molecules for the treatment against SARS-CoV-2 infections with one or more lung diseases by molecular docking analysis, a total number of n = 184 candidate drug molecules ( S2 Table ) were accumulated from different published articles associated with SARS-CoV-2 infections and/or different lung diseases. Before going to the molecular docking analysis, both the receptor-proteins and drug/ligand structures were prepared for molecular docking. The crystallographic structure of the target protein was obtained from the Protein Data Bank (PDB) and AlphaFold databases [ 6 ]. The ‘PubChem’ database was used to retrieve the 3D structures of those meta-drug agents. To identify potential binding sites within the target protein’s catalytic site (target pocket), we employed the ‘PrankWeb’ tools which employs a template-free machine learning method called ‘P2Rank’ for predicting ligand binding sites on solvent-accessible protein surfaces and identifying receptor protein target pockets [ 29 ]. Then, ‘PyMOL’ was used to visualize and select amino acid residues for the active sites of target proteins [ 24 ]. AutoDock tools 1.5.7 were utilized to process the receptor-proteins by removing water molecules, adding charges, and centering the grid box on the active site with specific coordinates [ 24 ]. The drug agents/ligands were pre-processed by minimized energy through the ‘Avogadro’ software and setting torsion tree using ‘AutoDock tools 1.5.7’. Both prepared receptor and ligands were converted into the PDBQT format. Subsequently, using the ‘AutoDock Vina’ [ 24 ], the binding affinities between the drug agents and target proteins were determined. Let A i j be the binding affinity score (BAS) between i th target protein ( i = 1,2,…, m ) and j th drug agent ( j = 1,2,…, n ). The target proteins and drug agents were arranged according to the descending order of row sums ∑ j = 1 n A i j ; i = 1 , … , m , and column sums ∑ i = 1 m A i j ; j = 1 , … , n , respectively, to pick the best-performing drugs as the candidate drugs. Then, we constructed the image of BAS corresponding to the ordered drugs and targets using TBtools [ 30 ]. Finally, the discovery studio visualizer (v21.1.0.0) and ‘PyMOL’ were used to display receptor-ligand interactions involving amino acids and their interactive position in the docked complex molecule. At first, the ligands were prepared by utilizing LigPrep module of Schrödinger [ 31 ] which involves the minimization of the appropriate bond angles and using force fields OPLS3 to minimize the structure’s energy. Next, several steps were performed to prepare proteins using protein preparation wizard tools. These steps involve adding hydrogen, utilizing prime to complete side chain gaps, replacing any absent loops with prime, removing water molecules located farther than 5.00 angstroms from the specified HET group, and producing the protonation state of heteroatoms using Epik, with a pH range of 7.0 ± 2.0. Finally, docking analysis was performed between the active site of the proteins and the drugs. Drug compounds were optimized using Gaussian 09 [ 32 ], a powerful computational chemistry software suite renowned for performing geometry optimizations and electronic structure calculations. The resulting data were then visualized and analyzed using GaussView 05 [ 33 ], a user-friendly interface that seamlessly integrates with ‘Gaussian 09’. Key parameters such as the highest occupied molecular orbital (HOMO) energies, lowest unoccupied molecular orbital (LUMO) energies, electron affinity, and electrophilicity index were calculated. These parameters play a crucial role in explaining the magnitude of ligand interaction in the binding pocket of the receptor protein. To evaluate the stability and flexibility of the top-ranked protein-ligand complexes, we considered two software’s (YASARA and Desmond) for molecular dynamics (MD) simulation for a duration of 100 nanoseconds. To investigate the dynamic behavior of the top-ranked protein-ligand complexes, MD simulations were performed using the AMBER14 force field by the YASARA Dynamics software [ 34 , 35 ]. Before starting the simulation, the hydrogen bonding network of the target-drug combination was solvated and tuned by a TIP3P water model [ 36 ]. Using the steepest gradient approach (5000 cycles), a simulated annealing method was used for the initial energy minimization of each simulation system. A 100 ns MD simulation was carried out under a Berendsen thermostat and constant pressure [ 37 ]. The YASARA macro’s default script and SciDAVis ( http://scidavis.sourceforge.net/ ) were used for this analysis. Then, using the following formula. B i n d i n g f r e e E n e r g y = E p o t Re c e p t o r + E s o l v Re c e p t o r + E p o t L i g a n d + E s o l v L i g a n d − E p o t C o m p l e x − E s o l v C o m p l e x the binding free energy of each snapshot was calculated by the molecular mechanics Poisson–Boltzmann surface area (MM-PBSA) function of the YASARA software. Initially, the Schrödinger software’s system constructor wizard was used to solvate the protein-ligand complex (PLC) with water molecules [ 38 ]. The study employed a transferable intermolecular potential 3P (TIP3P) solvent system [ 39 ]. The computational analyses were conducted within an orthorhombic box of 10 × 10 × 10 Å, utilizing periodic boundary conditions. A constant salt concentration of 0.15 M was maintained. The entire model system was subjected to a 100 ps minimization process, and the resulting trajectory was utilized for conducting a molecular dynamics investigation. The dynamics investigations were conducted for a duration of 100 nanoseconds, with a recording interval of 50 picoseconds. This resulted in a total of 5000 frames for the whole study. The simulations utilized a time step of 2 femtoseconds. The binding free energy of each snapshot in the protein-ligand interactions during high-throughput molecular dynamics (MD) simulations was calculated with MM-GBSA using gmx_MM-PBSA tools. The binding free energy (Δ G bind ) was calculated by using the following equation: Δ G b i n d = E c o m p l e x − ( E p r o t e i n + E l i g a n d ) This comprehensive approach encompasses multiple interaction-free energies, including van der Waals forces, electrostatic interactions, polar solvation effects, solvent-accessible surface area (SASA) contributions, and binding energies. In this study, we have employed g_mmpbsa package [ 40 , 41 ]. The trajectories from the protein-ligand MD simulation in explicit water from the Desmond module were used to generate the GROMACS trajectory file required for calculations using Schrödinger scripts. Additionally, topology files for protein and ligand were obtained separately by converting the *.cms files to *.gro and *.top files using the InterMol software [ 42 ].

(DOCX) (DOCX) (DOCX) (DOCX) (XLSX) (DOCX) (DOCX) (XLSX) (XLSX) (XLSX) (XLSX) (TIF)