Identifying potential drug targets and candidate drugs for COVID-19: biological networks and structural modeling approaches

Data collection

The coronavirus infected host gene expression profiles were extracted from the GEO database¹⁹. Agilent-014850 whole human genome microarray 4x44K G4112F platform was employed to GSE100509, GSE86529, GSE47962, GSE47961, GSE47960, GSE45042, GSE37827 and GSE33267 studies and one study GSE56677 which was performed based on Agilent-039494 SurePrint G3 human GE v2 8x60K microarray 039381 platform was included^20–24. All the selected studies contain control and CoV infected human lung epithelial cells samples. The samples infected with other viruses were excluded. The sample details of different datasets are indicated in Table 1.

Differential expression of CoV infected host genes

Using NetworkAnalyst (NA) web interface, the tab-delimited text (.txt) files of selected datasets were given as an input for preprocessing, normalization, and probe identification²⁵. Variance stabilizing normalization (VSN) and quantile normalization was applied to reduce false-positive errors and equal distributions of datasets for statistical analysis^26,27. R package “LIMMA” (Linear models for microarray analysis) of from the Bioconductor project were used to perform differential expression of CoV infected host genes²⁸. Log2 transformation, Benjamini and Hochberg and t-test were used to perform normalization and calculate false discovery rate (FDR; p<0.05) of samples²⁹.

Construction of DE-CoV-HGs PPI network, modules and PDTs identification

We constructed the PPI network from identified DE-CoV-HGs using Search tool for retrieval of interacting genes/proteins (STRING) interactome³⁰. The highest confidence interaction score was set to 0.9, which reduces false positive interactions³¹. The random walks (R package “igraph”) was used to extract modules based on the Walktrap algorithm from the DE-CoV-HGs interaction network³². It runs several short random walks within a group of nodes that are highly connected to detect small modules. From the modules, the hub genes (nodes) were identified using two different centrality measures “degree” and “betweenness”³². The degree of the gene is the many connections it has to other genes. Genes with a high degree act as hubs within the network. The betweenness of a gene is the number of paths that pass through it when considering the pair-wise shortest paths between all genes in the network. A node that occurs between two dense clusters will have a high betweenness. In the present study, the PDTs of DE-CoV-HGs were identified based on the degree centrality and betweenness centrality measures.

Pathway enrichment analysis

We have used ClueGO v2.5.3, which is a Cytoscape plugin for function and pathway enrichment analysis of DE-CoV-HGs³³. A list of PDTs were provided as input into ClueGO with select specific parameters like species such as-Homo sapiens; ID type-Entrez gene ID; enrichment functions -KEGG pathways for the analysis. Each enrichment was calculated based on the Bonferroni method (kappa score 0.96; p > 0.005). ImageGP was employed to visualize the results of functional enrichment analysis.

Structural modeling approach to study host-SARS-CoV-2 proteins interaction and drug-target interaction

Preparation of host and SARS-CoV-2 viral receptor. The FASTA sequence of selected proteins namely 3C-like protease (3cLpro), angiotensin I converting enzyme 2 (ACE2), B-cell lymphoma 3 (BCL3), casitas b-lineage lymphoma (CBL), CCAAT enhancer binding protein beta (CEBPB), CRK like proto-oncogene (CRKL), histone deacetylase 9 (HDAC9), Kelch-like protein 12 (KLHL12), v-myc myelocytomatosis viral oncogene (MYC), mothers against decapentaplegic homolog 3 (SMAD3), nucleocapsid protein (Npro), nuclear receptor coactivator 3 (NCOA3), receptor tyrosine-protein kinase erbB-4 (ERBB4), spike protein-receptor binding domain (SRBD), SMAD-specific E3 ubiquitin protein ligase 1 (SMURF1), ubiquitin C (UBC), and vascular endothelial growth factor A (VEGFA) were obtained from Uniprot for protein-protein docking studies (Table 3)³⁴.

Protein-protein docking. PPD was performed in the HDOCK web-interface using ab initio docking with a fast Fourier transform (FFT) based hierarchical algorithm^35,36. The FASTA sequence of the selected protein (protein 1 and protein 2) was given as input. The server obtained and chosen one of the highest sequence coverage and similarity and with the highest resolution template automatically from PDB and then, the homology modeling was performed using MODELLER. Then, the corresponding structures constructed by superimposing the modeled protein structure on to the template, where 5,000steps of AMBER MD minimization was also performed to remove severe atomic clashes at the interfaces. The putative binding modes between the protein-1 and protein-2 are systematically samples based on the FFT algorithm. In total 4,392 binding modes are ranked depending on their cluster and binding energy scores. The final complex was selected with RMSD of two binding modes was below a cut off 5Å. RMSD was calculated based on the backbone atoms of the protein 2³⁶.

Whereas, “a” indicates as iterative step. “x” and “y” indicates the types of a pair of atoms in the protein 1 and protein 2. ${\mathit{h}}_{\mathit{x}\mathit{y}}^{\mathit{o}\mathit{b}\mathit{s}}\mathit{(}\mathit{z}\mathit{)}$ is the pair distribution function for atom pair xy calculated from the ensemble of experimentally measured protein 1- protein 2 complex structures. ${\mathit{h}}_{\mathit{x}\mathit{y}}^{\mathit{a}}\mathit{(}\mathit{z}\mathit{)}$ is the pair distribution function calculated from the ensemble of the binding modes of the protein 1 – protein 2 complexes predicted with the trial potentials {${\mathit{h}}_{\mathit{x}\mathit{y}}^{\mathit{a}}\mathit{(}\mathit{z}\mathit{)}$} at the a-th step. ${\mathit{u}}_{\mathit{x}\mathit{y}}^{\mathit{a}\mathit{+}\mathit{1}}\mathit{(}\mathit{z}\mathit{)}$ are the improved potentials from ${\mathit{u}}_{\mathit{x}\mathit{y}}^{\mathit{a}}\mathit{(}\mathit{z}\mathit{)}.{\mathit{a}}_{\mathit{B}}\mathit{T}$. is the generality value set as 1³⁷.

Molecular dynamics simulations. The results of the PPD complex with the highest binding affinity were used to prepare input files for MD simulations. Solution builder in the CHARMM-GUI web interface was employed to solvate the complex with TIP3P and neutralized by potassium chloride (KCl) ions at the 0.15 mol-1 concentration under CHARMM36 force field^38,39. The initial configuration of KCl ions is then estimated using short Monte Carlo simulations (2000 steps) through Coulombic and the van der Waals (vdW) interactions. The solvated complex was simulated by using Nanoscale molecular dynamics (NAMD)⁴⁰. Full system long-range Coulombic interactions were determined by the particle mesh Ewald (PME) method⁴¹. The integrator parameters time was set to 2 fs/step. The simulations were executed in the constant pressure and temperature (NpT ensemble) with pressure bar (1) and temperature (303 K) through the Langevin dynamics with a damping coefficient of 1/ps. After 9×10⁷ steps of minimization, water molecules and ions were equilibrated for 2 ns around proteins, which were restrained using harmonic forces. MD simulation of selected drug target (ERBB4-3cLpro) was started from the last frame of restrained equilibration. The production was carried out for 150 ns.

Identification of possible candidate drug for PDT and molecular docking

Drug Gene Budger a web-interface was employed to identify possible candidate drugs or small molecules in LINCS L1000 data that significantly regulated selected drug target (ERBB4) expression based on LIMMA algorithm⁴². While entering the gene symbol of drug target, we could retrieve the following information such as log-transformed fold change, p-value and q-value for each potential small molecule. Small molecules/drugs that are liable for considerable overexpression and under expression of downregulated and upregulated drug target can increase the therapeutic action against COVID-19. The threshold q value was derived from the following equation

Where X is the number of false discoveries and Z is the number of true discoveries. Y=X+Z is the number of rejected null hypotheses. The false discovery rate (FDR) is as follows

Where E [Q] is the expected value of Q. The drugs/small molecule have the highest log2FC value was selected for further protein-docking studies. The drug target ERBB4 sequence was retrieved from the Uniprot database (ID: Q15303) for three structure modeling. The protein structure modeling was performed in a Swiss modeling server. The structure of the candidate drug “wortmannin” was retrieved from the ZINC database (ID: ZINC1619592). Then, protein-ligand docking was performed to calculate the binding energy of the docked complex in the Swiss dock web-interface, which was developed based on the EADock dihedral space sampling (DSS) algorithm. The energy functions were calculated using the CHARMM force field. Each docking cluster was the output of 250 different sequential runs. The binding models (BM) having the most favorable energies were calculated by fast analytical continuum treatment of solvation (FACTS) and clustered. Binding modes were estimated through Full fitness and clustered. Then, the Full fitness of the clusters was estimated by averaging 30% of the most favorable effective energy of their elements.

The effective energy calculated by the following formula

Whereas, E_intra^Ligand, E_intra^receptor are the internal energies of ligand and receptor; E_inter is interaction energy between the receptor and ligand; σ value (0.0072 kcal/mol); SASA is the solvent accessible surface area. The results were visualized in Maestro Visualizer. All the supporting data associated with this study freely available (http://doi.org/10.5281/zenodo.4458252).

Differentially expressed -CoV-HGs

In total, 1,968 DE-CoV-HGs in human lung epithelial cells were obtained from the initial analysis. To obtain better results of biomarker identification, we have selected 315 overexpressed and 112 under expressed DE-CoV-HGs based on the log-fold changes -1>log2FC>1. Table 2 illustrates the top 10 DE-CoV-HGs, their description, and their fold change and p value.

Modules and potential drug targets identification

PPI network was constructed using 427 DE-CoV-HGs (Figure 1a). The major subnetwork contains 531 edges with an average local clustering coefficient 0.353 and the PPI enrichment p-value was 2.36 ×10-3. Then, using the Walktrap algorithm, 47 modules were retrieved from the subnetwork. Among them, highly interconnected modules includes module 1 (No. of nodes=202; p=5.20×10^-32); module 2 (No. of nodes=126; p=1.90×10^-28) and module 3 (No. of nodes=121; p=3.13×10^-26) was employed to identify the PDTs (Figure 2b–d). Table 3 demonstrated that list of PDTs identified from different modules based on the degree and betweenness centrality measures. BCL3, CBL, CEBPB, CRKL, HDAC9, KLHL12, MYC, SMAD3, NCOA3, ERBB4, SMURF1, UBC, VEGFA were identified as PDTs which are highly expressed in human cells after CoV infection.

Figure 2. PPI network and Walk trap modules from DE-CoV-HGs.

(A) PPI network with highly interconnected nodes; (B) Module-1; (C) Module-2; (D) Module-3. Overexpressed genes are represented in green color, under-expressed genes in red color, interconnected nodes in gray color.

Functional and pathway enrichments of PPI network and PDTs

Many important gene ontology terms and pathways were identified from the PPI network namely, cellular response to growth factor stimulus, wound healing as a biological process term; growth factor receptor binding as a molecular function term; and TNF signaling pathway and transcriptional mis-regulation in cancer as a pathway enrichment terms (Figure 3). Moreover, the biological process term that is directly or indirectly related to the viral life cycle and human immune system includes RNA polymerase II promoter in response to hypoxia, peptidyl-tyrosine auto-phosphorylation, ERBB2 signaling pathway, anaphase-promoting complex dependent catabolic process, interleukin-1 mediated signaling pathway, TAP dependent antigen processing and presentation of peptide antigen via MHC class I, regulation of cellular amino acid and metabolic process, negative regulation of protein acetylation, negative regulation of gene silencing by miRNA, white fat cell differentiation and histone H4 deacetylation enriched in three different modules. Besides, stimulatory c-type lectin receptor signaling pathway (module 1), regulation of hematopoietic stem cell differentiation, T cell receptor signaling, innate immune response signaling and Fc-epsilon receptor signaling pathway (module 2), lymphocyte co-stimulation (module 3) were highly enriched immunological terms of three different modules (Figure 4). Further, pathways enrichment of modules demonstrated that numerous signal transduction pathways and cytokine immune signaling pathways were illustrated in Figure 5.

Figure 3. GO and Pathway enrichments of the PPI network.

(A) Biological process term; (B) Molecular function term; (C) Pathway enrichment terms of PPI network

Figure 4. Functional enrichments of the Walk trap modules.

(A) Biological process terms of three different modules; Immunological terms of (B) Module-1; (C) Module-2; (D) Module-3.

Figure 5. Pathways enrichment of Walktrap modules.

(A) Signal transduction pathways; (B). Cytokine immune signaling pathways. The figures are retrieved from Reactome pathway analyzer using PTGs.

Molecular interaction of host drug targets with SARS-CoV-2 protein

Protein-protein docking interaction. The potential target genes encoding proteins are defined as drug targets to control the pathogenesis of SARS-CoV-2. The drug targets regulated (overexpressed or under-expressed) by SARS-CoV-2 during pathogenesis in the host cytoplasm and different cellular components. Accordingly, PPD of SRBD, 3cLpro, and N-protein of SARS-CoV-2 with proposed drug targets happened by different biological signaling pathways. Thus, we hypothesized, proposed drug targets may control the pathogenesis of SARS-CoV-2 through the respective signaling pathways. The possible molecular interaction of host drug targets with SRBD, 3cLpro, and N-protein of SARS-CoV-2 was determined by PPD method. In the present study, we performed protein-protein rigid-body docking, which treats both proteins as rigid and discovers only six degrees of translational and rotational freedom. It excludes any kind of flexibility. Drug targets include HDAC9, SRBD, SMURF1, SRBD, CBL, ERBB4, BCL3, CEBPB, CRKL, KLHL12, MYC, NCOA3, SMAD3, UBC, and VEGFA were in the human respiratory system, which was possibly able to interact with the S-RBD, 3cLpro, and N - protein of SARS-CoV-2 and regulate the viral pathogenesis through replication and development. PPD binding affinities were illustrated in Table 3. Among the drug targets, ERBB4 demonstrated showed high binding affinity -260.46, -314.57 and -291.72 kcal/mol with SRBD, 3cLpro and N-protein of SARS-CoV-2 respectively. Thus “ERBB4” was selected as potential drug target and used for candidate drug selection study. Hydrogen bond forming residues between the ERBB4 and SRBD, 3cLpro, and N-protein of SARS-CoV-2 complexes represented in Figure 6.

Figure 6. Possible action and PPD of drug target ERBB4 with SARS-CoV-2 proteins.

(A) Possible pathological action of SRBD, 3cLpro, and N-protein of SARS-CoV-2 with “ERBB4”. Results of PPD of (B) ERBB4-SRBD; (C) ERBB4-3cLpro; and (D) ERBB4-N-protein complex by using LigPlot. The interacting residues of ERBB4 (magenta) are labelled in green color and SARS-CoV-2 protein structures (brown) are labelled in blue color. The interacting hydrogen bonds indicated in green color.

MD simulation of ERBB4-3cLpro complex. MD simulations demonstrated the stability of the ERBB4-3cLpro complex. The structure and dynamic properties of the complex were followed by monitoring the backbone Cα root mean square deviation (RMSD) from the initial structure during the simulation period (150ns). The RMSD of the ERBB4-3cLpro complex gradually increased for the first 95 ns before steadily oscillating around a value of ~0.85 nm, indicative of the ERBB4-3cLpro complex stability over the simulation timescale (Figure 7a). The root mean square fluctuation (RMSF) reflects the fluctuations in the positions of the ERBB4-3cLpro complex residues. Briefly, notable fluctuations have been observed in the multiple furin-like cysteine-rich domains, a tyrosine kinase domain, and a phosphatidylinositol-3 kinase binding site of ERBB4. However, the other amino acid residues in the PDZ binding domain (amino acid range 400–600) were not deviating more (Figure 7b). This indicates the reason for the stable RMSD plot of the ERBB4-3cLpro complex. The strong hydrogen bonding interaction (average ~1.89) between ERBB4 and 3cLpro observed throughout the 150ns of MD simulation (Figure 7c).

Figure 7. MD simulation of ERBB4-3cLpro complex.

(A) NAMD trajectories of RMSD; (B) RMSF; (C) Total number of intermolecular hydrogen bond formation trajectories of ERBB4-3cLpro complex.

Drug candidate for proposed drug target

The results of Drug Gene Budger demonstrated that Wortmannin (q- value = 3.31×10^-21; log2 Fold change = 2.997; Specificity = 1.09×10^-04) was identified as the candidate drug for the proposed target ERBB4. The ERBB4 was downregulated after the coronavirus infection that was then upregulated by the possible use of the Wortmannin as candidate drugs, which may inhibit the pathological cycle and development of SARS-CoV-2 in the human hosts. 30 different protein-ligand cluster of binding modes were obtained after molecular docking. Among them, the binding mode of cluster 1 indicated that the ERBB4-Wortmannin complex has the energy in term of full fitness -2563.6987 kcal/mol and Gibbs free energy ∆G was -7.93 kcal/mol (Figure 8).

Figure 8. Molecular interactions of Wortmannin within the catalytic pocket residues of receptor ERBB4.

(A) Two dimensional (hydrogen bond formation between Cys523 residue of ERBB4 and OH atom of Wortmannin showed in magenta line) and (B) three-dimensional interaction (hydrogen bond formation between Cys523 residue of protein-ligand complex showed in yellow color dotted line).

Precision based therapeutic countermeasures for COVID-19 patients

In the present study, the highly downregulated gene ERBB4 was selected as a drug target to control SARS-CoV-2 infection. Figure 9 illustrates the hypothesized mechanism action of downregulated ERBB4 in SARS-CoV-2 infected cells. A key receptor tyrosine kinase ERBB4 contains four domains includes multiple furin-like cysteine-rich domains, a tyrosine kinase domain, a phosphatidylinositol-3 kinase binding site, and a PDZ binding domain. Ligands are binding to the ERBB4 which induces different cellular response including alveolarization, alveolar epithelial type 2 cells (AEC2) differentiation, anchoring of cell surface receptor, and morphogenesis via extracellular-signal-regulated kinase pathway^81,82. Downregulation of ERBB4 inhibits regular action namely alveolarization, morphogenesis, and AEC2 differentiation. Many proteolytic events allow for the release of ERBB4 cytoplasmic fragments. The downregulation of ERBB4 by the SARS-CoV-2 in ACE2 cells may induce receptor tyrosine kinase-mediated macropinocytosis for the host cell entry. Besides, activation of macropinocytosis can increase non-specific fluid uptake, the pH in the vesicle is not decreased and, in that scenario, SARS-CoV-2 may be recycled to the cell surface. Interestingly, fluid-phase uptake (macropinocytosis and cell-to-cell spread) was reported in murine coronavirus and SARS-CoV infection^83–85. Both macropinocytosis and PI3K signaling are directly associated with the progression of the number of tumor cells including pancreatic and lung cancer^86,87. Further, growing evidence suggested that PI3K signaling involves the generation of the second messenger lipids and activates Rac-Pak signaling which plays a critical function namely macropinocytosis, mTOR signaling clathrin-mediated endocytosis, and cytoskeletal rearrangements^88,89. A recent systematic review by 16 and his colleagues showed the increased level of biomarkers such as C-reactive protein, serum amyloid A, interleukin-6, lactate dehydrogenase, D-dimer, cardiac troponin and renal biomarkers were higher in severe complicated COVID-19 patients’ plasma and infected lung tissues. Another study reported that low levels of type I and III interferons reduced innate antiviral defenses, which act as driving features of COVID-19⁹⁰. In addition, proteome analysis of SARS-CoV-2 infected hosts demonstrated that proteins involved in translation, splicing, carbon metabolism, proteostasis and nucleic acid metabolism could act as a drug target for SARS-CoV-2⁹¹. Moreover, Gordon et al. reported 67 druggable targets and repurposed 69 FDA-approved drugs for COVID-19⁹². This information supports the technical concept of the present study which enhances drug repurposing strategy and advances on precision-based therapeutics to COVID-19.

Figure 9. Hypothesized mechanism of action of downregulated ERBB4 in SARS-CoV-2 infected cells.

There are three possible mechanisms (1) SARS-CoV-2 an enveloped single-strand RNA virus, virions are internalized probably by receptor tyrosine kinase (EGFR) family receptor ERBB4 dependent macropinocytosis like SARS-CoV-1 (the process of macropinocytosis indicated in magenta circle). The cellular entry of SARS-CoV-2 involves binding to ERBB4 expressing AEC2 in the plasma membrane, macropinocytosis, and cathepsin based cleavage of the spike glycoproteins. The host transmembrane serine protease 2 (TMPRSS2) for spike priming and cathepsin cleaves spike glycoprotein S1 domain, permits fusion with endosomal membranes and release virions genome into the host cytoplasm; (2) Single-strand RNA translated polypeptides cleaved by the 3cLpro of SARS-CoV-2 possibly attached with mitochondrial ERBB4 and involve viral replication through MAPK signaling pathway. Then, the cell to cell spread action, the replicated SARS-CoV-2 attached to the uninfected AEC2 receptor and then restart macropinocytosis; (3) The part of S2 domain presents peptides into antigen-presenting cells (APC) to modulate host immune system through NFkB- signaling.

In the present study, we repurposed Wortmannin as a candidate drug that can inhibit PI3K signaling and receptor tyrosine kinase (ERBB4) mediated macropinocytosis, to control virus-cell proliferation and may enhance ERBB4 functions namely alveolarization, morphogenesis, and AEC2 differentiation. Wortmannin is a covalent inhibitor of phosphatidylinositol-3-kinase (PI3K) and related enzymes include the mammalian target of rapamycin (mTOR), DNA-dependent protein kinase catalytic subunit (DNA-PKcs), and mitogen-activated protein kinase (MAPK)^93,94. In humans, PI3K activates cell signaling pathways in the head and neck, urinary tract, cervical, ovarian, and lung cancer. Studies reported that the inhibition of PI3K signaling enhances tumor suppression and anti-tumor activity⁹⁵. Besides, Wortmannin is an extensively used cell biology research includes to inhibit cell proliferation, DNA repair, and receptor-mediated endocytosis^96,97. This report supported that the repurposed Wortmannin may have the potential to control SARS-CoV-2 infection in COVID-19 patients with downregulated expression of receptor tyrosine kinase ERBB4.

Underlying data

Zenodo: Datasets for SARS-CoV-2 drug target and candidate drug identification, http://doi.org/10.5281/zenodo.4458252⁹⁸

This project contains the following underlying data:

Data are available under the terms of the Creative Commons Attribution 4.0 International license (CC-BY 4.0).