Introduction

As of the end of October 2020, the novel Coronavirus Disease 2019 (COVID-19) caused by severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) has crossed over 50 million cases globally with more than 1 million deaths. While SARS-CoV-2 is known to cause substantial complications in respiratory and pulmonary systems, causing pneumonia and acute respiratory distress syndrome (ARDS), several COVID-19 cases present various extrapulmonary manifestations of COVID-19. These include manifestations of cardiovascular, renal, hematologic, gastrointestinal, hepatobiliary, endocrinological, ophthalmological dermatological, and neurological systems [1,2,3,4]. To date, increasing evidences are indicating the neurological manifestations of the central and peripheral nervous system in COVID-19 [5]. These neurological complications include headaches, dizziness, nausea, loss of consciousness, seizures, anorexia, anosmia, ageusia, encephalopathy, and meningo-encephalitis [6,7,8,9,10,11]. In addition, ischemic stroke [12], acute necrotizing encephalopathy (ANE) [8], and acute inflammatory demyelinating polyneuropathy (Guillain-Barré syndrome, GBS) [13] have also been reported to be associated with COVID-19. Moreover, COVID-19 patient autopsy studies have also identified viral RNA transcripts and viral proteins in brain tissues [14, 15].

The potential mechanisms underpinning the various neurological syndromes include direct or indirect viral neuronal injury [16], a secondary hyperinflammation syndrome related to cytokine storm [17], post-infectious immune-mediated disorders, or the effects of a severe systemic disorder with the neurological consequences of sepsis, hyperpyrexia, hypoxia, vasculopathy, and/or coagulopathy.

Recently, two studies utilized the 3D human brain organoids to study the neuroinvasive potentials of SARS-CoV-2. Ramani et al. [18] have revealed that SARS-CoV-2 readily targets cortical neuronal cells and induces tauopathies and neuronal cell death. The other study by Song et al. [19] models the SARS-CoV-2 infection of neuronal cells in hiPSC-derived brain 3D organoids. These comorbidities thus lead to higher risk of disease development and higher mortality associated with COVID-19. It is largely assumed that other neurological manifestations will likely surface in the near future. Thus, understanding the infection mechanisms of SARS-CoV-2 and associated comorbidity is a high priority to contain its short- and long-term effects on human health.

SARS-CoV-2 causes the disease by hijacking the host cell machinery and perturbs the highly organized cellular networks. Moreover, the highly coordinated interactions between molecules observed in a healthy cell are also gets altered in a disease condition [20]. This suggests that SARS-CoV-2 interaction with host human cells will be different in healthy and diseased cells, and thus could lead to different impacts of COVID-19 infection. Therefore, pre-existing clinical conditions can facilitate the appearance of another disease if they share the same or functionally related genes [21]. As SARS-CoV-2 has been shown to be associated with several neurological manifestations, we here predicted the risk of COVID-19 infection in patients with various neurological disorders.

Over the last 20 years, numerous human-viral interactomes have been constructed to understand the mechanisms of viral entry, infection, and disease progression [22,23,24,25,26]. Investigating such interactomes has led to the discovery of shared and distinct molecular pathways associated with viral pathogenicity. In the present work, we have utilized a network-based system biological framework (Fig. 1) to investigate the molecular interplay between COVID-19 and various human neurological disorders. Here, we have constructed a brain-specific protein-protein interaction network of the 332 genes of human’s targeted by SARS-CoV-2 reported by Gordon et al. [27], with their neighboring genes and named the network as COVID-19 target network (CTN). The genes in the network were then further used to identify the neurological disease associated with them. Based on the shared genes, we have integrated the CTN with brain diseases and generated a disease-gene network of the brain (BDGN). This human brain disease-gene network consists of a total of 123 various brain disorders including COVID-19 interacting with 653 genes of CTN. Out of these 123 diseases, 28 diseases were directly linked with the COVID-19, indicating the comorbidity and complexity of COVID-19. Next, we have identified the functional modules and hub genes of the CTN network that can be considered as the hotspot for comorbidity and could be the target for drug repurposing. We therefore emphasized that targeting these functional modules will inhibit the viral replication and growth and thus will improve the medical conditions in comorbidity associated with COVID-19.

Fig. 1
figure 1

A strategic workflow adopted in this study with self-explained legends

Material and Methods

Brain-Specific Protein-Protein Interaction (PPI) Network for COVID-19 Target Genes

Interactome data of the human brain was retrieved from the TissueNet v.2 database with 165,240 interactions [28]. Quantitative tissue association for human PPI is provided by TissueNet. TissueNet gathers RNA-Seq and protein-based assay profiles from the genotype-tissue expression project (GTEX) and human protein atlas (HPA) for preparing extensive interaction networks. Experimentally validated protein interaction information extracted from four major databases (DIP, BioGrid, MINT, and IntAct) was also included for preparing PPI networks. A list of 332 genes of humans known to interact with COVID-19 was obtained from Gordon et.al [27] and was used for creating a subnetwork having interaction between the COVID-19 target genes and their neighboring genes. The subnetwork was called the COVID-19 target network (CTN). Out of 332 genes, 4 genes (GDF15, INHBE, SBNO1, and CEP43) did not show any interactions in the brain.

Assembling Brain-Specific Disease-Gene and Disease-Disease Interaction Network

After obtaining the list of COVID-19 target genes and their neighboring genes from CTN, Gene ORGANizer [29] tool was used to identify the brain-related disorders linked with these genes. Gene ORGANizer database allows us to analyze the relationship between the query genes and the organs affected by them. The database provides organ-specific disease-gene relation from highly curated DisGeNET and human phenotype ontology (HPO) tools. Disease-gene interactions having valid HPO identifiers were considered for further study. A total of 2002 disease-gene interactions related to brain disorder was obtained, which included 127 various brain disorders interacting with 653 genes. Forty-three out of 653 genes were the direct target of the COVID-19. Finally, a disease-gene interaction network was prepared which consists of 780 nodes (i.e., diseases/genes) and 2046 interactions. Disease-Disease interaction network was then prepared using the disease-gene interaction network in which two diseases are considered to be associated only if they share a common gene known to cause those diseases. A network of 123 diseases connected with each other was prepared, in which 28 diseases were directly connected to the COVID-19. The calculation of the topological coefficient of the disease-gene interaction network also indicated the tendency of nodes in the network to have shared neighbors [30]. The values of topological coefficient also suggest that CTN may have functional modules within the network, as the genes in the modules are more connected in-between rather than the genes present outside the module [31].

Identification of Functional Protein Modules of CTN

Organization of biological function is believed to be in a modular and hierarchal manner [32]. It is believed that cellular function modular organization is reflected as modular structure in molecular network [33]. Protein modules are groups of highly connected proteins that tend to be co-localized and co-expressed. Proteins in protein-modules form a single multimolecular machine by interacting with each other at the same time and place [34, 35]. The MCODE plugin of the Cytoscape tool was used for identifying protein modules in CTN with degree cut-off value equals to 2, node score cut-off value equals to 0.2 and k-core value equals to 2. MCODE algorithm operates on three stages, i.e., vertex weighing, complex prediction, and optionally post-processing step, for identifying the locally highly connected region in the interaction network [36]. A total of 32 modules were obtained, from which the five largest modules in terms of their sizes have been selected for further studies. These selected five modules have the maximum number of genes and also have multiple COVID-19 target genes in it. These five functional modules thus cover a major part of the COVID-19 target network and most likely will play an important role in the progression of the disease. The functional modules identified in CTN are large, so it is important to identify which genes in each module best explains its behavior. A widely used approach is to identify highly connected genes (also known as hub genes) in the modules [37].

Identification of Hub Genes of CTN

Hub genes in the network are the genes that are highly connected with other genes in the network on a direct basis. Any change in the expression or activity of the hub gene has the potential to influence the working of the network. Hubs are frequently more relevant to the functionality of the network than other genes. We calculated the topological properties like the degree of connectivity (K), betweenness centrality value, and closeness centrality for all the 5 functional modules, similar to our previous study [26]. All these network topology parameters were calculated using the network analyzer plugin of the Cytoscape tool [38]. Briefly, degree (k) signifies the number of interactions made by nodes in a network and is expressed as:

$$ \mathrm{Degree}\ \mathrm{of}\ \mathrm{connectivity}\ \left(\mathrm{k}\right)={\sum}_{v\varepsilon {K}_u}w\left(u,v\right) $$

where Ku is the node-set containing all the neighbors of node u, and w(u,v) is the edge weight connecting node u with node v.

Betweenness centrality (Cb) represents the degree to which nodes stand between each other based on the shortest paths. A node with higher betweenness centrality represents more control over the network. It is expressed as:

$$ {C}_b\left(\mathrm{u}\right)={\sum}_{k\ne u\ne f}\frac{p\left(k,u,f\right)}{p\left(k,f\right)} $$

where p(k,u,f) is the number of interactions from k to f that passes through u, and p(k,f) denotes the total number of shortest interactions between node k and f.

Closeness centrality (Cc) is a measure of how fast information is traveled from one node to other nodes in the network. Closeness centrality value ranges from 0 to 1, and isolated genes have closeness centrality value equal to zero.

$$ {C}_c(z)=\frac{1}{avg\left(L\left(z,m\right)\right)} $$

where z is the node for which the closeness value is calculated and L(z,m) is the length of the shortest path between two nodes z and m. It has been seen that genes having a high degree of connectivity also have high closeness centrality score.

Eccentricity (eG) is defined as shortest distance between one node from all the other nodes present in the network. The shortest the distance between the nodes, the faster the travel of information among them. If a node is an isolated node, its eccentricity value will be zero.

$$ {e}_G(v)=\max \left\{{dist}_G\left(v,u\right):u\ \epsilon\ V(G)\right\} $$

where v is the central node in a connected network G for which we are calculating the eccentricity and u are the nodes presented in the network.

Topological coefficients (Tf) indicate the tendency of the nodes in the network to have shared neighbors. Nodes having no or one neighbor are assigned a topological coefficient of zero. The topological coefficient of a node, n with kf neighbors is computed as:

$$ {T}_f=\frac{avg\left(j\left(f,p\right)\right)}{k_f} $$

where j(f,p) is the number of shared neighbors between f and p, plus 1 if there is an edge between f and p.

The top 35 genes having a high degree of connectivity and betweenness values were considered as hub genes of the CTN. Apart from topological parameters, we have also calculated the Dyadicity, which measures the connectedness of the nodes belonging to the same groups, and the Jaccard similarity coefficient to check the extent of molecular overlapping among COVID-19 and other neurological disorders.

Identification of Key Genes and Key Regulators of the Protein Modules

After identifying the hub genes of the CTN network, we have identified one key gene in each of the functional modules. A key gene is among the hub genes of the modules which also plays role in connecting the functional modules with each other and provides the idea about how these modules shares information with each other. Furthermore, for identifying the key regulator or driver genes which control the regulation of the modules, tracing of hub genes of each module up to motif level was performed. Key regulators are the genes that are highly connected in the network not only when seen from the top of the network but its connection reaches the lower levels in the network making the gene an important and influential part of the network. Apart from identifying the important genes in the modules, it is equally important to identify the processes and the pathways in which these modules play role, to better understand how the alterations in the function of these modules affects the host condition.

Gene Ontology and Pathway Enrichment Analysis

For enrichment analysis of protein-modules, Database for Annotation Visualization and Integrated Discovery (David v.6.7) tool was used [39]. Gene Ontology (GO) enrichment analysis includes annotation at biological, cellular, and at molecular levels. DAVID uses the GO and Kyoto encyclopaedia of genes and genomes (KEGG) database for the enrichment analysis of the genes. Pathways and functions having P value less than 0.05 were considered significantly enriched.

Network-Based Drug Repurposing

For identifying the drug targets for hub genes, key genes, and key regulators, databases such as DrugBank, Clue.io [40], ChemblInteraction [41], and DGIdb [42] were screened out. A final drug-gene interaction network was prepared using the STITCH database [43]. STITCH is a database known to predict the physical and functional interaction between the query genes and drugs. The interactions in the STITCH database are derived from five sources namely, automated text mining, high throughput lab experiment data, co-expression interaction data, interaction prediction by genomic context, and by previous knowledge from other databases. For each interaction, STITCH calculated a combined score. A combined score is calculated by combining the corrected probability of observing an interaction randomly and probabilities of interaction from different evidence channels. The drug-gene interactions having a combined score value > 0.7 were considered high confidence interactions.

Apart from chemical molecules as drugs, we also identified the microRNAs (miRNAs) as a potential drug candidate, the miRNAs interacting with the 35 hub genes were retrieved from the miRTarBase database [44]. mirTarBase database provides miRNA interaction with several hosts organism including humans. This database uses several validation methods such as reporter assay, western blot, qPCR, microarray, NGS, and pSILAC to validate the interaction between host genes and miRNAs. For 35 hub genes, 1997 miRNA-gene interactions were retrieved. Interactions having at least one strong evidence of validation methods were used for further analysis.

Results

Protein-Protein Interaction Network Between COVID-19 and Human Host in the Brain

For constructing COVID-19 and its human interaction network, we retrieved the protein-protein interaction network of the brain from TissuevNet2.0 [28] database. A network consisting of 12,968 number of nodes and 165,241 number of edges was constructed using the brain’s PPI data (Supplementary Fig. 1). A list of 332 human target genes of COVID-19 was retrieved from Gordan et.al [27]. A subnetwork of these 332 genes with their neighboring genes was constructed from the brain’s PPI network and was named as the COVID-19 target network (CTN) (Fig. 2a, Supplementary Table 1). Out of the 332 COVID-19 target genes, 327 genes were part of the subnetwork. The subnetwork has 5061 number of nodes and 95,802 number of edges. More than 50% of interactions from the main network become the part of a subnetwork which clearly depicts how deeply the COVID-19 target genes are connected in the brain. The degree distribution of CTN indicates that the network has a scale-free property (Fig. 2b). A scale-free network is one that follows a power-law distribution and is independent of the size of the network (i.e., the number of nodes in the network), which means the basic structural foundation of the network remains the same even when the network grows [45, 46]. Next, we also calculated the dyadicity (D) between the COVID-19 target genes in CTN and obtain a dyadicity value as 24.2. Dyadicity refers to the connectedness of the nodes belonging to the same group in a network. In a dyadic network, nodes belonging to the same group are more connected with each other than in a random network. A network is called dyadic when D > 1 [47, 48]. A value of 24.2 in our analysis indicates that COVID-19 target genes are very well connected and making a complex in the network, which can hijack the host cellular machinery. The complexes in a network are likely to contribute to the progression of disease and comorbidity at the molecular level [49]. Proteins present in a community have a higher probability to play role in comorbidity as compared to the other proteins which are not part of the complex. Therefore, to understand the high risks of comorbidities associated with covid19, we have prepared and briefly analyzed the disease-gene and disease-disease interaction network of CTN.

Fig. 2
figure 2

a COVID-19 target genes (in red) interaction network in the brain with neighboring genes (in green). b Scatter plot showing the distribution of degree (k) in the COVID-19 target network (CTN)

Disease-Gene and Disease-Disease Interaction Network in the Brain

A disease-gene interaction network and disease-disease interaction network from CTN were prepared to understand the links of COVID-19 with neurological comorbidities. A disease-gene interaction map of CTN was constructed from the disease-gene interaction data of the brain obtained from the Gene ORGANizer database. A total of 313 brain diseases, 715 genes, and 2965 disease-gene pairs were considered for network construction. The disease-gene association map in CTN will then be constructed by connecting the associated node (genes) and brain disorder. The resulting network revealed 653 genes linked to a total of 127 disorders, which also includes COVID-19 (Fig. 3a, Supplementary Table 2). Figure 3a represents the disease-gene network map comprising of 780 nodes and 2002 disease-gene interactions and is termed as the brain disease-gene network (BDGN).

Fig. 3
figure 3

a Disease gene interaction network: the figure represents the interaction of CTN network genes with their related brain disorder. COVID-19 target genes are represented in red, neighboring genes of COVID-19 target genes are represented in green, diseases are represented by blue color, and COVID-19 disease is represented by yellow color. b Dot plot of highly connected diseases and the number of genes associated with disease in the brain’s gene-disease interaction network. c Bar plot of genes highly connected to multiple diseases in the brain’s gene-disease interaction network. d Brain’s disease-disease interaction network. Red color nodes represent the disease directly connected to COVID-19

The disease-gene interaction network showed that several disorders were connected with more than one gene in the network such as ataxia (k = 245), dementia (k = 91), autism (k = 69), and COVID-19 (k = 43) (Fig. 3b). Similarly, the disease-gene interaction network also showed that many of the disorder share the common genotype. For example, SARS-CoV-2 targets, TBK1 (k = 16), BCS1L (k = 7), DNMT1 (k = 7), FBN1 (k = 7), WFS1 (k = 7) and neighborhood nodes, CDKL5 (k = 17), MECP2 (k = 16), VCP (k = 16), TARDBP (k = 14), and ATXN2 (k = 12) are linked to multiple disorders (Fig. 3c, Supplementary Table 3). Also, the calculation of the topological coefficient of the BDGN network indicated that several nodes in the network have shared neighbors (Supplementary Fig. 2). Thus, the resulting BDGN reveals the molecular connection of COVID-19 with various brain disorders and also the close association of COVID-19 targets with the genes causing the brain disorders (Supplementary Fig. 3).

To further demonstrate the association between COVID-19 and brain disorders, a disease-disease association network was then constructed, where two diseases were considered to be related if they share one common gene. The resulting disease-disease interaction network contains a total of 123 diseases (nodes) and 436 edges, representing a higher clustering between diseases (Fig. 3d and Supplementary Table 4). Out of 123 disease nodes in the network, 28 nodes (in red) were directly linked with COVID-19 (yellow node) (Fig. 3d). Jackard similarity coefficient was also calculated to check the extent of molecular overlapping between COVID-19 and other brain-related disorders. Several diseases like ataxia, dysarthria, spasticity, cerebral atrophy, autism, dementia, and stroke were closely related to the COVID-19 and are also connected with multiple genes in the brain (Supplementary Fig. 4A and 4B). Thus, because of these molecular overlapping, patients having these diseases are more prone to COVID-19 and vice-versa. The higher molecular similarities between brain diseases create a highly complex high-density comorbidity cluster which contributes to higher mortality in COVID-19 patients.

To further strengthen the above statement of the high-density comorbidity, we calculated the eccentricity value of the network. The eccentricity of a node in a biological network can be interpreted as the easiness of the node to be functionally reached by all other nodes in the network [50]. A node with a high eccentricity value compared to the average eccentricity value of the network will be more easily influenced by the activity of the other nodes and conversely could also easily influence the other nodes in the network by its activity. We observed that more than 80% of the nodes in the disease-disease network shows the eccentricity value greater than or equal to the network’s average eccentricity value (i.e., 4), which represents that the functionality of the nodes in the network is highly linked to each other making the network a highly complex cluster (Supplementary Table 5).

Recent studies evidently showed the extra pulmonary associations in COVID-19 that contribute to higher mortality in COVID-19 patients. Thus, it is of utmost necessity to develop effective drugs to target the patient-specific risks of comorbidity during COVID-19 infection. However, presence of several overlapping molecular connection makes it difficult to identify and prioritize the targets for the treatment of the COVID-19 infection. We therefore focused next to target the host functional protein modules linked with diverse brain diseases.

A Broad Range of Disorders Are Linked to Functional Protein-Modules of Host Interaction Network

It is generally accepted that biological networks are not randomly connected and follow a structural pattern which gives rise to the modular structure and hierarchical organization. Modularity suggests nodes that are highly connected in a community are most likely to have the same biological functions and play a role in similar pathways [51]. Many complex networks exhibit modular structures, where in-between modules interactions are less dense as compared to interaction within modules. These modules reflect the organization of the functional unit in a network with relative independence [52]. Generally, diseases sharing similar genes are more predisposed to form disease modules and comorbidity. Similarly, genes related to similar diseases are likely to highly interact with each other [53]. For identifying protein modules, we used the MCODE module of the Cytoscape tool. A total of 32 protein modules in CTN were identified by MCODE. The top 5 modules having a larger number of nodes and edges were selected for further studies (Supplementary Table 6). Module-1 (Fig. 4a) had 257 nodes and 1542 edges including 15 COVID-19 target genes (in red). Module-2 (Fig. 5a) had 253 nodes and 921 edges including 17 COVID-19 target genes. Module-3 (Fig. 6a) had 183 nodes and 619 edges including 14 COVID-19 target genes. Module-4 (Fig. 7a) had 175 nodes and 899 edges including 10 COVID-19 target genes and module-5 (Fig. 8a) had 155 nodes and 340 edges including 5 COVID-19 target genes. The presence of a high number of COVID-19 target genes in these protein modules indicated that during the infection, these modules might be hijacked and strongly altered as compared to other modules in the network and will eventually disrupt the majority of the network function. DAVID tool was used for gene ontology analysis and KEGG pathway analysis.

Fig. 4
figure 4

a Module1 with COVID-19 target genes (in red). b Biological functions and KEGG pathways related to module-1. c Dot plot of module-1 related disease with their gene counts

Fig. 5
figure 5

a Module-2 with COVID-19 target genes (in red). b Biological functions and KEGG pathways related to module-2. c Dot plot of module-2 related disease with their gene counts

Fig. 6
figure 6

a Module-3 with COVID-19 target genes (in red). b Biological functions and KEGG pathways related to module-3. c Dot plot of module-3-related disease with their gene counts

Fig. 7
figure 7

a Module4 with COVID-19 target genes (in red). b Biological functions and KEGG pathways related to module-4. c Dot plot of module-4-related disease with their gene counts

Fig. 8
figure 8

a Module-5 with COVID-19 target genes (in red). b Biological functions and KEGG pathways related to module-5. c Dot plot of module-5-related disease with their gene counts

Biological functional enrichment analysis of module-1 revealed its role in RNA splicing, mRNA processing, protein complex biogenesis and assembly, regulation of programmed cell death and pathway analysis reveal module-1 role in proteolysis, ribosome pathway, spliceosome pathway, in Huntington’s disease pathway, and cancer pathways (Fig. 4b). The nodes of module1 were associated with disorders like ataxia, dysarthria, spasticity, encephalopathy, coma, and delayed speech and language development along with many other disorders (Fig. 4c).

Similarly, biological functional analysis indicates that module-2 plays a role in macromolecular complex assembly, negative regulation of gene expression, RNA transport, and localization, and play role in pathways like spliceosome, ribosome, cell cycle, gliomas, pathways in cancer, and RIG-I like receptor signaling pathways (Fig. 5b). The nodes are associated with disorders like hydrocephalus, dementia, autism, muscular dystrophy, language impairment, and frontotemporal dementia (Fig. 5c).

Enrichment analysis of module-3 reveals its role in RNA splicing, cell division, mRNA processing, spindle assembly, and nuclear division-related biological process, whereas pathway analysis reveals roles in Notch signaling pathways, Toll-like receptor signaling pathways, and SNARE interaction in vesicular transport (Fig. 6b). Nodes of the module-3 are associated with diseases like ataxia, sleep apnea, leukodystrophy, cerebral ataxia, and Joubert syndrome (Fig. 6c).

Module-4 plays a role in biological functions like RNA processing, translational elongation, ncRNA processing, and in the viral infection cycle, and also plays a role in pathways like NOD-like receptor signaling pathways, Parkinson’s disease, Huntington’s disease, and Cytosolic DNA sensing pathways (Fig. 7b). Module-4 is associated with diseases like ataxia, spasticity, amyotrophic lateral sclerosis, chorea, and auditory neuropathy (Fig. 7c).

Module-5 biological functional analysis reveals its role in the regulation of transcription from RNA polymerase II promoter, chromatin modification, regulation of transcription, and in cellular protein catabolic process, pathways enrichment analysis reveal a role in Wnt signaling pathways, viral reproduction, MAPK signaling pathway, neurotrophin signaling pathway, TGF-beta signaling pathway, and ErbB signaling pathway (Fig. 8b). Module 5 is associated with diseases such as ataxia, myopathy, Chiari malformation, GAIT ataxia, limb ataxia, and febrile seizures (Fig. 8c).

Interestingly, during the early stage of infection, SARS-CoV-2 hijack these modules and the concerned biological processes and synthesize its RNA. Many clinical conditions such as ataxia, spasticity, encephalopathy, and dementia are associated with all modules, revealing a greater risk of the severe illness of COVID-19 patients. The presence of common disease, the pathway, and the biological functions in the modules indicated that these modules are connected with each other to some extent (Supplementary Table 7). Besides, we also observed a different spectrum of disorders like cancer, myopia, dysarthria, sleep apnea, and thromboembolic stroke that were associated with these modules. Therefore, besides neurologic comorbidity, disorders in various other organs can also be a potential threat for COVID-19 patients as also demonstrated by Gysi et al. [54]. Our network-based results will thus add to strengthen these observations.

Core Regulatory Hubs of the Protein Network

Several virus-host network studies have indicated that viral proteins target the hub proteins which have a high degree of connectivity in the network [22, 24,25,26, 55]. Here, we have identified 35 candidate hub proteins that exhibit a high degree of connectivity and betweenness value, using the Cytoscape’s Network analyzer tool (Supplementary Table 8). Out of these 35 hub proteins, 16 proteins bind to SARS-CoV-2 [27]. Biological process enrichment analysis of these 35 hub proteins using GeneMania webserver and GO biological process analysis indicated their role in proteasome-mediated ubiquitin-dependent protein catabolic process, G1/S transition of the mitotic cell cycle, Notch signaling pathway, ERBB signaling pathway, and response to hypoxia (Fig. 9 and Supplementary Table 9). The GO function analysis reveals the roles of hub proteins in nitric-oxide synthase regulator activity, ubiquitin-protein ligase binding, transcription regulator activity, RNA binding, and cytoskeletal protein binding (Supplementary Table 9).

Fig. 9
figure 9

Protein-protein interaction network of 35 hub genes derived from GeneMANIA along with functional enrichment

Strikingly, 26 of the SARS-COV-2 viral proteins interact with the above-identified 16 hub genes (Table 1). As can be seen, SARS-CoV-2’s NSP8, NSP13, ORF3a, and ORF9c target the most of these hub proteins (i.e., eight in total), SARS-CoV-2 NSP2 and ORF7a targets seven hub proteins, while SARS-CoV-2 N, M, NSP6, and NSP7 have six hub protein targets (Table S3). Other SARS-CoV-2 proteins like NSP’s (4, 9, 12) and SARS-CoV-2 ORFs’ (3b, 8, 10) have five targets. Intriguingly, all the 16 hub genes are targets of more than one SARS-CoV-2 protein (Table 1), out of which NPM1, SNW1, GRB2, HSP90AA1, and ELAVL1 are the target of several SARS-CoV-2 proteins. These findings are consistent with previous findings which indicates that an individual viral protein can target multiple host proteins and several viral proteins can interact with the same host protein [25, 26, 55,56,57].

Table 1 Interaction of 16 hub genes of the CTN with SARS-CoV-2 proteins

Moreover, the spatio-temporal expression analysis of these 35 hub genes was studied through the BEST (Brain Expression Spatio-Temporal) webserver (http://best.psych.ac.cn/#) [58]. This server utilizes eight human brain expression datasets obtained from BrainSpan Atlas, Allen brainmap, GTEx, and other sources. The expression pattern of selected genes in different brain regions (spatial pattern) and age stages (temporal pattern) were analyzed using RNA-seq Data from Brainspan and RNA-seq Data from GTEx and has been shown in Supplementary Fig. 5. As can be seen in the figure that most of the genes were upregulated from the neonatal stages to adulthood, these genes were downregulated in the hypothalamus, olfactory bulb, substantia nigra, insula, and parahippocampal region (Supplementary Fig. 5A). However, in older age, the genes were moderately upregulated in the hypothalamus and substantia nigra, whereas they largely get downregulated in the olfactory bulb, thalamus, and in the cortical region (Supplementary Fig. 5B). The downregulation of these genes in the olfactory bulb, thalamus, substantia nigra, and the cortical regions of the brain indicates the impairment of sensory systems, memory, and cognition.

Network-Based Drug Repurposing

Because of high connectivity and the ability to rapidly transfer information in the network, hub genes could be the most appropriate target in any network for drug identification [26]. We propose to target these hub proteins, hijacked by SARS-CoV-2, through drug repositioning. Initially, we identified drugs for the 35 hub genes by screening multiple databases related to drug-genes interaction [42, 59, 60]. A total of 3286 drug-gene interactions for 35 hub genes were identified (Supplementary Table 10). For 16 hub genes showing interaction with COVID-19 viral protein, 1477 drug-gene interactions were identified. Out of these 16 genes, 12 genes showed significant interactions with 41 drugs having a combined score > 0.7 in the STITCH database (Supplementary Fig. 6).

Considering the severity and complexity of COVID-19, we also target key genes connecting the functional modules using drug repurposing. Targeting these key genes makes much more sense, as information transmitted by these key genes will further provide instructions to other modules controlling the network. We have identified five key genes (ESR1, TP53, UBC, HSP90AA1, and MYC) that were linking the modules in-between and can act as a suitable candidate for drug repurposing (Supplementary Fig. 7A). A total of 138 number of FDA approved drugs were identified for the 5 key genes. A drug-gene interaction network was prepared, and the STITCH database was used for the final categorization of interactions. Interactions having a combined score value > 0.7 were considered as high confidence interactions (Supplementary Fig. 7B). Moreover, earlier clinical studies indicated that many of the drugs used for the treatment of viral CNS infections have poor CNS penetration abilities and thus were unable to inhibit viral RNA replication [61]. These drugs were shown to have lower log P values (CNS penetration ability) and hence were unable to cross the blood-brain barrier. Ideally, according to Lipinski’s rule of 5, a log P value less than 5 has been considered as a pharmacologically effective drug [62]. The log P value of the identified drugs has been shown in parentheses and this information would be highly useful for selecting the drugs for drug repurposing against COVID-19 induced CNS infection.

ESR1 gene showed high confidence interactions with FDA approved drugs, fulvestrant (logP, 6.54), tamoxifen (logP, 5.93), raloxifene (logP, 5.45), estradiol (logP, 3.57), and ethinylestradiol (logP, 3.63). TP53 gene showed interaction with etoposide (logP, 0.73), 5-fluorouracil (logP, − 0.89), mitomycin-c (logP, − 0.40), bortezomib (logP, 0.89), and doxorubicin (logP, 1.41) drugs. MYC gene showed interaction with imatinib (logP, 3.47), calcitriol (logP, 5.51), and amifostine (logP, − 1.4) drugs and HSP90AA1 gene showed interaction with rifabutin (logP, 4.25) drug.

Key regulators in a network are the proteins that are deeply rooted in the network/modules and serve as the backbone of the network. Since networks serve as hierarchical characteristics, removal of key regulators from the network may cause alarming changes in the whole network and especially in the deeper levels of the network [63]. Thus, these key regulators will be proved as a good target gene in the modules. A total of 12 key regulators in 5 modules were identified by tracing the hub proteins of each module up to the motif level (Fig. 10 a). For module 1, YWHAQ, CUL4B, CUL2, and HSPA8 genes were identified as key regulators. Similarly, for module 2, the HDAC1 gene was identified, for module 3, the KAT2B gene, for module 4, the MYH9 gene, and for module 5, SMARCA4, SMARCC2, MYC, H2AX, and GAN genes were identified as key regulators. Drug targets against these key regulators were identified. A total of 107 FDA approved drugs were identified for these key regulators. Out of these 12 key regulators, 6 key genes showed significant interaction with drugs (combined value of > 0.7) using the STITCH database (Fig. 10b). For the HSPA8 protein of module-1, significant interactions were obtained with compounds such as hydrogen peroxide (logP, − 0.45), ibuprofen (logP, 3.97), arsenite, and DB07045 (logP, 1.39). Similarly, the HDAC1 protein of module-2 showed interactions with trapoxin B (logP, 3.1), trichostatin A (logP, 2.36), and vorinostat (logP, 1.88) drug molecules. Module-3 protein, KAT2B showed interaction with R4368 (N-(3-AMINOPROPYL)-N-(2-NITROPHENYL)) (logP, 1.6) drug molecule. MYH9 protein of module-4 showed interaction with blebbistatin (logP, 2) and module 5 protein, MYC showed interactions with trichostatin A (logP, 2.36), troglitazone, tamoxifen, etoposide, doxorubicin, and arsenite. Another protein of module-5, H2AX showed interaction with vorinostat, etoposide, and doxorubicin. Many of the drugs targeting these key regulators can be used either individually or in combination.

Fig. 10
figure 10

Identification of key regulators of CTN. a Tracing of hub genes through different levels in the modules. The red arrow represents the transfer of hub genes to the next level. b Interaction network of key regulators with FDA-approved drugs. The green color interaction between the drug and genes was a significant interaction with a combined score > 0.7

Finally, we searched the miRNAs target for the 35 hub proteins that may be considered as a potential drug to target the comorbidity. Among these 35 proteins, 19 proteins were targeted by 157 miRNAs. Of which, 11 miRNAs including, miR-125a-5p, miR-125b-5p, let-7a-5p, miR-130a-3p, miR-34a-5p, miR-29b-3p, miR-27a-3p, miR-24-3p, miR-145-5p, miR-200a-3p, and miR-200c-3p were noticed to target more than three interconnected proteins (Fig. 11). These miRNAs thus may show its utility as a drug for therapeutic intervention.

Fig. 11
figure 11

miRNA network. The network shows the miRNA (red node) targeted 19 hub genes (green node). The miRNA (yellow highlighted node) that binds to three or more than three hub genes are shown only

Discussion

The most alarming situation in COVID-19 infection is the associated comorbidity in the aged patients which increases the severe health risks worldwide. The present work has shown the risk of COVID-19 infection on the onset of various neurological disorders and the molecular basis of this comorbidity through the principle of system-based network biology. We here showed the complexity of COVID-19 and wide range of SARS-CoV-2 targets in the host cell, which establishes the molecular connection with various brain-related disorders. The disease-gene and disease-disease network map has revealed the overlapping molecular connections and high clustering of diseases within the vicinity of the same network, thus demonstrating a close pathobiological similarity. Moreover, the formation of a scale-free network among brain disorders and COVID-19 indicates the presence of a high-density comorbidity cluster. These results increase our understanding of the molecular basis of the neurological comorbidity associated with COVID-19 patients.

Another significant finding is the existence of functional modules that exhibit greater connectivity among nodes within a module (Figs. 4, 5, 6, 7, 8). Therefore, targeting the functional protein modules which are primarily hijacked by SARS-CoV-2 and are the origin of many brain disorders will prevent the virus replication and growth. The most highly connected module relates to the ubiquitin-proteasome system (UPS), ribosomal proteins, cell signaling, and cellular export proteins that can be considered as critical targets for reducing the SARS-CoV-2 infection. Furthermore, in-depth network analyses revealed 35 hub proteins present in different functional modules and are involved in more than one function (Supplementary Table 8). To provide the system-wide importance of these hub proteins in COVID-19 and associated comorbidity, we categorized these hub proteins into three groups based on their functionality. The proteins in the first group have been largely relevant to the UPS, cell cycle, and cell death. Similarly, the proteins in the other two groups are involved in transcription regulator activity, RNA metabolic pathway, signaling pathway, nuclear transport, and cytoskeletal protein binding.

The group 1 proteins are generally multifunctional and contain genes responsible for highly represented UPS. These are CUL1, CUL2, CUL3, CUL7, CAND1, OBSL1, VCP, FBXO6, UBC, and RNF2. Four of these genes (CUL1, CUL2, CUL3, and CUL7) belong to the core component of E3 ubiquitin-protein ligase complexes, which mediate the ubiquitination of target proteins. Three other UPS-related proteins, VCP, FBXO6, and UBC are the key mediators of the ER-associated degradation (ERAD) pathway. The UPS is largely involved in the maintenance of cellular homeostasis and plays a fundamental role in viral replication [64, 65]. Several viruses, including coronaviruses, often modulate ubiquitin and ubiquitin-related pathways for their survival [66, 67]. In this study, we find that module 1 and module 2 includes several members of ubiquitin-ligase and ER-associated proteins (Fig. 4a and 5a). Since the loss of ubiquitin-proteasome activity results in aggregation of proteins and disturbs cellular functions, we suggested that SARS-CoV-2 targets these modules to hinder the UPS-mediated cellular responses. The other genes of this group are involved in the regulation of cell cycle and cell death. It was recently suggested that COVID-19 infection may cause neurodegenerative disease and leads to neuronal death. The study revealed that SARS-CoV-2 infection-induced pathological effects closely resembling tauopathies and neuronal cell death [18]. In our study, we found genes like APP, TP53, MYC1, VCP, and UBC as hub genes that are involved in protein misfolding and aggregation, which ultimately leads to cell death.

Moreover, the second largest group of hub genes is involved in transcription regulation and RNA binding. The genes in this group include MOV10, XPO1, ESR1, SIRT7, NPM1, ELAVL1, NXF1, TP53, CDC5L, SNW1, RNF2, and COPS5. A recent network study identified XPO1, NPM1, HNRNPA1, and JUN, as “hub proteins” with the highest number of functional connections within 119 host proteins that interact with the human SARS-CoV-2 [56]. Of note, both NPM1 and HNRNPA1 interact directly with XPO1 in normal cells, and JUN is involved in inflammatory responses, including those associated with viral infections.

A recent single-cell RNA sequencing of the SARS-CoV-2 infected human brain organoid showed that SARS-CoV-2 infection induces a locally hypoxic environment in neurons [19]. Also, post-mortem studies of the brain of COVID-19 patients indicated the acute hypoxic-ischemic damage in the brain [14]. Consistent with these results, our findings indicate the enrichment of the hypoxia-inducing genes in the BDGN, including CUL2, TP53, UBC, and MDM2, suggesting a possible mechanism of SARS-CoV-2 induced neuropathology.

We then mapped the known drug-target network to search for druggable, cellular targets. We identified 41 approved targets for the 16 hub genes. For example, TP53, CDK2, HSP90AA, HDAC1, NPM1, and ESR1 were the most targetable proteins. One of the important drug targets is the export protein, XPO1, which enables the transport of viral proteins from the nucleus to the cytoplasm. The FDA-approved drug, Selinexor (logP, 2.85), (Fig. 12a) is a selective inhibitor of nuclear export (SINE) compound which blocks the cellular protein XPO1 and is now entered into a randomized clinical trial in COVID-19 patients as announced by Karyopharm Therapeutics Inc. [68]. Selinexor is FDA-approved drug against multiple myeloma. SINE compounds have been used to block the replication of several viruses and act as anti-inflammatory drugs [69]. In a study, verdinexor (logP, 3.7), which is highly similar in structure and biological activity to selinexor, was identified as a candidate drug inhibiting the interactions of several of the SARS-CoV-2 proteins to their human targets and was recommended for further evaluation to target COVID-19 [27].

Fig. 12
figure 12

Drug-gene interaction network of two important key genes and its drugs using the STITCH database. a The network shows selinexor, verdinexor, and guanosine triphosphate as drug candidates for XPO1. b The network shows blebbistatin as a candidate drug that binds to MYH9 and could be considered for drug repurposing

The other key drug target gene is MYH9 which has been traced to a candidate gene of module- 4 largely involved in viral infection. The MYH9 gene is involved in synthesizing a protein called myosin-9 which is one subunit of the myosin IIA protein. The MYH9-related disorder includes thrombocytopenia, hearing loss, kidney disease, and cataracts [70]. Mutations in the MYH9 gene lead to a genetic condition known as MYH9-related thrombocytopenia (MYH9RD) which is characterized by the presence of large but lesser platelets, results in bleeding in the body or the skin [71]. Thrombocytopenia is now considered as a risk factor for increased morbidity and mortality in COVID-19 patients. According to a recent study, COVID-19 patients showed blood clots in small vessels, clots in the lungs, and stroke-causing clots in cerebral arteries [72]. In some cases, a stroke may also lead to a ruptured blood vessel, and cause bleeding in the brain [73]. We here found that the chemical compounds bind to MYH9 includes guanosine triphosphate, MgADP, MgATP, and blebbistatin (logP, 2) (Fig. 12 b). The drug blebbistatin is a myosin-2 inhibitor and has been commonly used drugs in several research disciplines including neuroscience, muscle physiology, cell migration, cell differentiation, cell death, and cancer research [74]. It has also been shown to reduce the viral infection in many of the viruses by inhibiting its entry [75] and reducing the viral loads [76].Here, we also suggests eleven micro-RNAs (miRNAs) which shows the binding with hub proteins and also regulates hub proteins. The miRNA has been considered as drug molecules for targeting the SARS-CoV-2 proteins [77]. It is largely accepted that host cellular miRNA can directly target both the coding region of the viral genome and 3’UTR to induce the antiviral effect. Very recently, Fulzele et al. [78] identified over 800 miRNAs that target the COVID-19 genome. Interestingly, many of the miRNAs identified in this study (e.g., miR-130a-3p, miR-29b-3p, miR-27a-3p, miR-145-5p, and miR-200a-3p) are found to be unique to the SARS-CoV-2 genome according to the results of the Fulzele et al. [78]. Moreover, miRNAs including miR-34a-5p, miR-200a-3p, and let-7a-5p have been shown to strongly interact with the SARS-CoV-2 [79]. Besides, polyphenols are also considered as a potential and valuable natural compound for designing new drugs that could be used effectively in the combat against COVID-19 [80,81,82]. Overall, our BDGN analysis opens up avenues for multiple candidate repurposable drugs that target diverse cellular pathways involved in COVID-19 and associated neurologic comorbidity.

Conclusions

In conclusion, we constructed a human-SARS-CoV-2 BDGN interactome, demonstrating the risk of COVID-19 infection on multiple brain-related disorders. The study of disease-gene and disease-disease association networks reveals COVID-19-related structural and functional modules that establish the molecular connection with brain disorders. Also, the highly clustered network of the diseases suggests close pathobiological similarity and large comorbidity. We also identified crucial functional hubs in the network that showed close links with brain-related disorders including dementia, ataxia, encephalopathy, stroke, and many other disorders that indicate a high degree of comorbidity and the severity of COVID-19. We here also propose candidate drugs and miRNAs, targeting the hub proteins interconnecting the molecular pathways linked to virus infection that can be used to treat SARS-CoV-2 infection. We believe our results will help to understand the molecular pathobiology of neurological comorbidities linked to COVID-19. Given the severity and complexity of SARS-CoV-2 infection, we further suggest that drugs in combinations or drugs targeting multiple proteins will be the choice to improve clinical outcomes.