Immune-Related Protein Interaction Network in Severe COVID-19 Patients toward the Identification of Key Proteins and Drug Repurposing

Coronavirus disease 2019 (COVID-19) is still an active global public health issue. Although vaccines and therapeutic options are available, some patients experience severe conditions and need critical care support. Hence, identifying key genes or proteins involved in immune-related severe COVID-19 is necessary to find or develop the targeted therapies. This study proposed a novel construction of an immune-related protein interaction network (IPIN) in severe cases with the use of a network diffusion technique on a human interactome network and transcriptomic data. Enrichment analysis revealed that the IPIN was mainly associated with antiviral, innate immune, apoptosis, cell division, and cell cycle regulation signaling pathways. Twenty-three proteins were identified as key proteins to find associated drugs. Finally, poly (I:C), mitomycin C, decitabine, gemcitabine, hydroxyurea, tamoxifen, and curcumin were the potential drugs interacting with the key proteins to heal severe COVID-19. In conclusion, IPIN can be a good representative network for the immune system that integrates the protein interaction network and transcriptomic data. Thus, the key proteins and target drugs in IPIN help to find a new treatment with the use of existing drugs to treat the disease apart from vaccination and conventional antiviral therapy.


Introduction
Coronavirus disease 2019 (COVID-19) has been spreading worldwide, despite several developed vaccines, still causing numerous cases. Moreover, most causes of death are from severe complications from the disease. Currently, the global statistical data from the World Health Organization (WHO) indicate that 476,374,234 and 6,108,976 cases are infected and dead, respectively (26 March 2022) [1]. COVID-19 is an infectious disease caused by severe acute respiratory syndrome coronavirus-2 (SARS-CoV-2), a positivesense single-strand RNA (+ssRNA) virus [2]. SARS-CoV-2 is a betacoronavirus, classified in the Coronaviridae family [3]. Although SARS-CoV-2 incurs a well-known pandemic coronavirus infection in the present, for severe respiratory diseases resulting from the two coronavirus diseases, most known cases emerged before, in the last two decades.
networks [53,[70][71][72]. Nonetheless, analyses of importance and difference of the centralities in biological roles are still needed for further investigation.
There are two methods for extracting a specific subnetwork: the neighborhood and diffusion approaches. Neighborhood method is a traditional technique performed widely in many biological network studies [73][74][75][76][77]. Furthermore, there are several COVID-19 network studies constructing protein-protein interaction networks and selecting a group of related proteins using this method [50,51,54,55,[78][79][80]. This method extracts subnetworks by considering the nodes with the low shortest path with disease-associated genes or seed nodes [81]. The benefits of neighborhood-based subnetwork construction are that the approach is user-friendly and has fast running time. However, some neighbor nodes in the networks are not disease-associated genes, causing topological changes and missing identification of key nodes. Therefore, the diffusion-based method is preferred to build specific protein-protein interaction subnetworks, such as immune, inflammation, and viralhost interaction network, with lower false-positive disease-related nodes, although it is time consuming and requires coding skills. Network diffusion is the method used to predict novel disease-associated genes based on known disease-associated genes via considering the diffusion or probability scores from iterative running algorithms at time stable. Nodes with high diffusion scores are inferred as theoretical disease-related genes [82].
Moreover, there have been no COVID-19 studies involved in immune-related biological networks using the network diffusion method. As a result, we proposed this method to construct an immune-related protein interaction network (IPIN). A network diffusion method named Laplacian heat diffusion (LHD) algorithm on a human interactome network was conducted in this study to construct an IPIN in severe COVID-19 patients. Key immune-related proteins in the network were identified using several centralities and ranking score measurement. Additionally, drug repurposing was also performed to find target medication to those key proteins. This study aims to discover candidate target drugs to treat severe COVID-19 at clinical levels.

Materials and Methods
The summary of materials and methods used in our study is illustrated as a diagram in Figure 1. First, differentially expressed genes (DEGs) of the transcriptomic data downloaded from GEO DataSets were mapped with the human interactome data from STRING v11.0 database [83] for forming the seed nodes. Second, construction of a human protein interactome network obtained from STRING database [83]. Third, the Laplacian heat diffusion (LHD) algorithm operated network diffusion. Fourth, a permutation test filtered out false-positive high diffusion score nodes. Fifth, significantly high diffuse score nodes were used to construct the IPIN. Metascape [84] was also performed for the functional enrichment analysis of the network. Sixth, Molecular Complex Detection (MCODE) and the Markov Clustering (MCL) algorithm were conducted to find IPIN modules. Next, the centralities and the ranking score were calculated to identify the key immune-related proteins. Finally, candidate drugs targeting the key proteins were discovered by chemical and drug databases searching for drug-gene and drug-protein interactions.

Human Protein Interactome
The human protein interactome data, containing 19,566 proteins (or nodes) and 11,938,498 interactions, was obtained from STRING v11.0 database (https://string-db.org/ accessed on 20 November 2021) [83]. R package 'dplyr' [85] was used to manipulate the data. The interactions with a combined score between 900 and 999 were included. R package 'igraph' [86] was conducted to eliminate isolated nodes and multiple edges. As a result, the rest of the interactome had 11,334 proteins and 123,263 interactions. The combined scores between 900 and 999 were changed to a weighted score by rescaling into the range of [0.01, 1].
1 Figure 1. Summary of the process to identify the key proteins and drug repurposing in the severe COVID-19 based on an immune-related protein interaction network (IPIN). Circles represent protein nodes in a protein-protein interaction network. Dark blue circles are nodes of the DEGs' proteins. Light blue circles are nodes of proteins having diffuse scores. Orange circles represent high diffusion score proteins for IPIN. Red, green, and yellow circles are proteins in different modules. Pink circles are target proteins of an existing drug.

Transcriptomic Data
The collected data was from Gill et al.'s study (2020) in GEO DataSets (GSE154998) (https://www.ncbi.nlm.nih.gov/geo/ accessed on 20 November 2021) [87]. The researchers collected leukocytes samples from COVID-19 cases and controls in an intensive care unit (ICU) to perform transcriptomic profiles by RNA sequencing (RNA-seq) method. The total sample size was 14, with 7 samples being COVID-19 cases and the rest being controls. The SARS-CoV-2-positive cases were confirmed using reverse transcription-polymerase chain reaction (RT-PCR) method. In addition, genes from RNA-seq data having the false discovery rate (FDR) < 0.05 and log2fold change (log2 FC) > 1.5 were differentially expressed genes (DEGs). Hence, there were 224 genes meeting the criterion (Table S1 in Supplementary Materials). The DEGs then were mapped with the protein list in the protein-protein interaction network of the STRING v11.0 database by using Ensembl ID joining [83]. Thus, there were 189 Ensembl protein IDs (Table S2 in Supplementary Materials).
With the use of network diffusion, the seed nodes were prepared from mapping between these 189 Ensembl IDs (the immune-related proteins from DEGs) and the human interactome data (in Section 2.1.1). This resulted in 141 seed nodes (Table S3 in Supplementary Materials).

LHD Algorithm
This study used LHD algorithm to operate human interactome network propagation. In many studies, LHD is one of the most common network diffusion algorithms used to infer disease-associated genes or proteins [88,89]. Given a network called G, let W be a weight adjacency matrix of network G calculated from rescaling the combined scores and D be a diagonal matrix whose values are a degree of each node arranged diagonally. Laplacian matrix (L) was calculated from D−W. An initial diffusion vector (H 0 ) of all nodes is conducted normally by setting the initial heat diffusion score (h 0 ) to each node. In general, the initial heat scores (h 0 ) are set as 1/n where n is the number of seed nodes related to the subject of interest (immune-related proteins as seeds in our case) while the other nodes are set as 0. The diffusion vector at time t (H t ) was updated based on the previous diffusion vector at time t−1 (H t−1 ) according to this equation: where e is defined as Euler's number (≈2.1828). The network diffusion was iterated based on Equation (1) until two consecutive diffusion vectors were relatively similar. In our case, the relative similarity is met if H t − H t−1 < 10 −6 , then the diffusion becomes stable. The latest diffusion scores after the stability are used as an indicator of the relevant scores to the seed nodes. In another word, a node with a high diffusion score was strongly associated with seed nodes. In our cases, the initial diffusion vector (H 0 ) contained the initial scores of 11,334 protein nodes in the whole protein-protein interaction network. The initial scores for 141 seed nodes were set as 1/141 while other protein nodes are set as 0. After the network propagation was stable, the final scores were then used for further permutation test analysis. In this study, LHD algorithm was carried out by using R package 'diffusr' [90]. The parameters in the package were well-established using the default setting.

Permutation Test
When LHD algorithm finishes the network propagation, in theory, the nodes with the high diffusion score are strongly associated with the disease-associated proteins. However, some nodes receiving the high diffusion score can be false-positive. This false-positive result can occur because some factors such as the topological structure of network G can provide a high diffusion score apart from the actual association with the seed nodes. Therefore, a permutation test should be conducted to filter out nodes with the false-positive high diffusion score. The permutation test measures whether nodes have a high diffusion score due to statistical significance or chance. The test was operated by assigning 1000 different sets of the initial seeds into the human interactome network for LHD algorithm. Before running the algorithm in each set, 141 seed nodes were randomly assigned, independent of DEGs, from 11,334 nodes. Hence, the sets of seed nodes in the original set and in the 1000 sets were different. The z-score of the diffusion score of node n (Z(n)) was calculated according to this equation: where h(n) is defined as the diffusion score of node n in the original set, X(n) and SD(n) are the mean and standard deviation, respectively, of diffusion score of node n in the original set and 1000 permutation sets. A node with a z-score more than 1.96 (p-value < 0.05) had the true-positive high diffusion score and was selected for the IPIN construction. Table S4 in Supplementary Materials shows the original diffuse scores including with the mean and standard deviation calculated from the diffusion score in the original set and 1000 permutation sets in each node in the whole protein-protein interaction network.

IPIN Construction
After the network diffusion by LHD algorithm and the validation of the high diffuse score nodes by permutation testing were performed, 154 nodes with significant diffusion scores were obtained (Table S5 in Supplementary Materials). The filtered significant nodes from the permutation test were mapped with the STRING human interactome network. The largest component of the significantly high diffusion score nodes consisting of 97 nodes and 778 interactions was selected as our immune-related protein interaction network or IPIN for further functional enrichment analysis, topological analysis, and centrality measurement. The construction of the IPIN was based on leukocytes or white blood cells' transcriptome data. Leukocytes are the cells that play an important role in pathogen defense mechanisms or immunity. The gene expression of leukocytes during the infection indicates the host's immune status against the pathogens. Thus, this IPIN was constructed to represent a core host immune system against the severe COVID-19. The edge list information for this IPIN is also provided in Table S14 in Supplementary Materials.

Topological Analysis, Network Centrality, and Ranking Scores
R package 'igraph' [86] was performed to calculate global and local topological parameters and network centrality [86]. The global topological parameters such as the number of nodes (N), the number of edges (M), average degree (<k>), diameter (D), mean shortest path length (mspl), and average clustering coefficient (acc) were computed. Furthermore, degree distribution and a clustering coefficient versus degree curve were plotted to evaluate the scale-free network properties. Several local parameters and centrality measurements such as degree centrality (C D ), betweenness centrality (C B ), closeness centrality (C C ), and eigenvector centrality (C E ) were calculated to find the essential nodes in the network. A node with more than 90th percentile centrality value was listed and its functional importance in the network and the disease was discussed. These high-value nodes in the four centrality sets were plotted using an upset plot from R package 'UpSetR' [91] to find the node overlapping in each centrality.
Degree centrality is the number of adjacent nodes having the interaction with interested node i, according to this equation where A ij is a value of a non-weight adjacency matrix (A) of node i and j, respectively. In network biology, the high-degree proteins are hub nodes and play an important role in the network function [92]. Therefore, numerous medications are designed for targeting the hub nodes. Betweenness centrality is the summation of the ratio of the shortest path between node u and v that passes through node i. The betweenness centrality equation is where σ uv is a total number of the shortest path between node u and v and σ uv (i) is the number of the shortest path between node u and v that pass through node i. The highbetweenness proteins indicate the bottleneck property, forming the bridges controlling the flow of information in the network. Interruption of the bottlenecks by several targeted therapies can result in network function destruction in many diseases, improving the disease outcomes [92][93][94][95]. Closeness centrality is the summation of inverse shortest path distance between node i and all other nodes in the network. The closeness centrality equation is shown as where N is the total number of nodes in the network and d(i, j) is the shortest path length between node i and node j. Some studies used closeness centrality to find essential genes and proteins in several biological problems [96,97]. Eigenvector centrality is the extended form of degree centrality that focuses on the global high-degree nodes more than the local high-degree nodes. Eigenvector centrality has the assumption that the essential nodes should connect to the other important nodes. Therefore, nodes with high eigenvector centrality have high degree centrality, and their neighbored nodes also have a high degree value. The equation is demonstrated as where λ is the largest eigenvalue of A. Eigenvector centrality is applied to analyze many connectome studies in neurological diseases [98][99][100][101].
As mentioned above, all four centralities play a vital role in discovering essential genes or proteins. Moreover, numerous traditional network biology studies usually use degree and betweenness centrality to find essential genes and proteins. Closeness and eigenvector are occasionally applied to find vital genes and proteins. Therefore, a node with a high value of all centralities was a key protein in the network. The ranking score (S R ) was applied to rank the nodes based on the centralities. Let C be a set of all centralities and c represent a centrality measure in C. A ranking score of any node i (S R (i)) was calculated based on the reciprocal of the product of a ranking position of node i in each centrality c (k(i, c)), according to this equation Nodes with high-ranking scores greater than the 90th percentile were considered as key proteins in the IPIN network. In addition, nodes with degree and betweenness centrality greater than the 90th percentile were used to compare the ranking scoring nodes.

Functional Enrichment Analysis and Network Clustering
Metascape (https://metascape.org/gp/index.html#/main/step1 accessed on 10 December 2021) [84] was conducted for functional enrichment analysis in the largest component network according to the six terms: Gene ontology biological process (GO-BP) [102], Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways [103], REACTOME pathways [104], WikiPathways [105], Canonical pathways [106], and CORUM pathway [107]. Moreover, the Molecular Complex Detection (MCODE) algorithm [108] in Metascape was operated to cluster the enrichment terms into large groups to find the more common biological terms. The Markov Clustering (MCL) algorithm [109] was also used in the STRING v11.0 database to explore the network communities in the IPIN. The inflation parameter of MCL was set to 1.5. Functional enrichment analysis in each module was also conducted in STRING v11.0 database using GO-BP, REACTOME pathways, KEGG pathways, and WikiPathways term. In addition, the clusters were visualized by using STRING v11.0.

Detection of Potential Drugs for Drug Repurposing
The key proteins, having a ranking score, degree, or betweenness centrality above the 90th percentile, were used as the input to find drug-gene and drug-protein interactions from DrugBank database (https://go.drugbank.com/ accessed on 15 December 2021) [110], Therapeutic Target Database (TTD) (http://db.idrblab.net/ttd/ accessed on 15 December 2021) [111], Comparative Toxicogenomics Databases (CTD) (http://ctdbase.org/ accessed on 15 December 2021) [112], and GeneCards (https://www.genecards.org/ accessed on 15 December 2021) [113]. Drugs which have United States Food and Drug Administration (FDA) approval and evidence of interactions with the key genes or proteins were selected. The STITCH v5.0 database (http://stitch.embl.de/ accessed on 18 December 2021) [114], containing drug-protein interaction information, was performed to confirm the chosen drugs. A confidence score of the interaction in the STITCH database was used to find suitable drug-protein interactions. The confidence score is the probability value calculated based on both experimental and computational evidence such as text mining, high-throughput experiments, co-expression and gene fusion data, and information from other databases. A drug with a confidence score of more than 0.9 was considered a candidate drug having efficiency for severe COVID-19 treatment.

IPIN Construction and Topological Properties
The immune-related protein interaction network, known as IPIN, was obtained after the network diffusion and the validation of the high diffuse score nodes by permutation test. This network consisted of 97 nodes and 778 interactions as shown in Figure 2. In addition, the STRING reports showed that the average node degree and the expected edges are 16 and 50, respectively.

IPIN Construction and Topological Properties
The immune-related protein interaction network, known as IPIN, was obtained after the network diffusion and the validation of the high diffuse score nodes by permutation test. This network consisted of 97 nodes and 778 interactions as shown in Figure 2 In addition, the STRING reports showed that the average node degree and the expected edges are 16 and 50, respectively.  Table 1 summarizes the global topological parameters of the IPIN. The average degree and diameter of the IPIN are 16.04 and 7, respectively. The network has the small-world effect because it provided the low mean shortest path length (mspl = 3.01) but the high average clustering coefficient (acc = 0.74). These behaviors are concordant with other biological networks. Local topological parameters in each node in the network, for example, clustering coefficient and degree, are summarized in Table S6 in Supplementary Materials. To find the scale-free property, a degree distribution was plotted to prove the power law, as shown in Figure 3a. Furthermore, a clustering coefficient versus degree plot is illustrated in Figure 3b. The degree distribution plot reveals that it does not follow the power-law distribution because it provides a low correlation (R 2 = 0.1). This appearance is explained by the IPIN is a real-world network extracted from the human interactome network. Thus, the power-law properties can be disrupted due to the subnetwork construction. Nevertheless, the clustering coefficient versus degree plot shows the independence between the clustering coefficient and degree (R 2 = 0.061). The independent relation between clustering coefficient and degree is found in random and scale-free networks [115]. Analysis in the STRING revealed that the IPIN's protein-protein interaction (PPI) enrichment p-value is significant (p-value < 10 −16 ), indicating that the proteins have interactions with each other more than by chance. Thus, the interactions in the IPIN were more significant than random interactions. The IPIN was less likely to be a random network even though there is no suitable reason to explain the scale-free properties of the IPIN.

Functional Enrichment and Module Identification in IPIN
The functional enrichment analysis was performed using a hypergeometric test and Benjamini-Hochberg statistical correction algorithm in Metascape [84], which integrates six biological and pathway enrichment terms: GO-BP [102], KEGG pathways [103], RE-ACTOME pathways [104], WikiPathways [105], Canonical pathways [106], and CORUM pathway [107]. The results revealed that the terms were primarily involved in immune pathways, cell division, nucleotide metabolisms, and protein processing, as shown in Figure 4. The immune pathways were associated with innate immune response and antiviral signaling pathways such as type I and II IFN (IFN-I and IFN-II) and IFN-stimulated genes (ISGs). IFN-I mainly comprises of IFN-α and IFN-β, while IFN-γ is a component in IFN-II. The cell division terms were relevant to mitotic cell cycle process, mitotic metaphase, regulation of cell cycle, phase transition of cell cycle checkpoint at G1/S and G2/M, cytoskeletondependent cytokinesis, and regulation of sister chromatid separation. Moreover, nucleic metabolic pathways, such as pyrimidine metabolism, DNA metabolic process, and regulation of DNA replication, were identified. Apoptosis was also found from the enrichment analysis. Other enrichment terms were protein and enzyme processing, such as protein modification by small protein conjugation, negative regulation of catalytic activity, protein tetramerization, and protein localization to organelle.   Table 2 shows the enrichment analysis results from the In addition, MCODE algorithm clustered the functional enrichment terms into four modules, MCODE1, MCODE2, MCODE3, and MCODE4, as shown in Figure 5. For instance, MCODE1 (Figure 5a) represented the biological term related to innate immune and proinflammatory cytokine signaling pathways: IFN-α and IFN-β. MCODE2 (Figure 5b) was associated with mitotic cell division, cell damage detection, and chromosome segregation. MCODE3 (Figure 5c) was involved in cell cycle checkpoint and cell cycle signaling pathways, and MCODE4 (Figure 5d) was enriched with protein processing and antigen presentation. Table 2 shows the enrichment analysis results from the MCODE algorithm.    Table  S8 in Supplementary Materials. The enrichments result was aggregable with the result from the MCODE analysis.   detail about the clusters is described in Table S7 in Supplementary Materials. The result from functional enrichment analysis of the four modules revealed that the MCL1 (Figure 6a) was related to cell cycle regulation functions while MCL2 (Figure 6b) was mainly enriched in innate immune responses. Moreover, MCL3 (Figure 6c) had an association with nucleic acid metabolism. MCL4 (Figure 6d) mixed cell division and immune response terms. The details of enrichment results in each MCL module are listed in Table S8 in Supplementary Materials. The enrichments result was aggregable with the result from the MCODE analysis.  Table 2). The innate immune response found in the study mostly correlated with IFN signaling pathways. IFN can be separated into two groups: type I and II IFN (IFN-I and IFN-II). IFN-I cytokines, such as IFN-α and IFN-β, play an essential role in innate immunity to many viral infections [116]. They are released by the interaction between PAMPs or DAMPs and PRRs. The IFN-I activation controls innate immunity and increases viral clearance by stimulating various antiviral proteins, such as MX1, OASs, and ISGs [117][118][119][120][121]. IFN-I is also involved in NK cells, B lymphocyte, CD4 + , and CD8 + T lymphocyte stimulation [121]. Nonetheless, IFN-I signaling is suppressed and delayed in coronavirus infections, such as SARS-CoV, MERS-CoV, and SARS-CoV-2, causing viral clearance dysfunction. Persistent viral replication causes the release of uncontrolled proinflammatory cytokines, resulting in the cytokine storm [116,122,123]. This phenomenon is also called paradoxical hyperinflammation. Although IFN-II, consisting of IFN-γ, has an immune function overlapping with IFN-I to stimulate the antiviral and innate immune response, it causes enhancement predominantly in antigen-presenting cells (APCs) [124]. Several studies have revealed that IFN-II exhaustion is usually found in severe COVID-19 patients, suggesting a vital role of IFN-I in the immunopathology of the disease [125][126][127]. As a result, IFN signaling enrichment found in the peripheral white blood cells can be from the compensating mechanism of immune cells for IFN suppression and delay in severe cases. However, IFN signaling predominance in this study can indicate the ongoing activation. Persistent IFN stimulation induces apoptosis of CD4+ T lymphocytes and causes lymphopenia. It also increases proinflammatory cytokine production [128][129][130]. Other evidence supports the notion that increased IFN-I in influenza viral infection can release excessive proinflammatory cytokines, resulting in respiratory epithelial apoptosis and severe pneumonia [131]. Because the IFN function in severe COVID-19 is complicated, the study of IFN roles in the disease should be more investigated.
Leukocyte proliferation is the immune defense mechanism responding to infec-  Table 2). The innate immune response found in the study mostly correlated with IFN signaling pathways. IFN can be separated into two groups: type I and II IFN (IFN-I and IFN-II). IFN-I cytokines, such as IFN-α and IFN-β, play an essential role in innate immunity to many viral infections [116]. They are released by the interaction between PAMPs or DAMPs and PRRs. The IFN-I activation controls innate immunity and increases viral clearance by stimulating various antiviral proteins, such as MX1, OASs, and ISGs [117][118][119][120][121]. IFN-I is also involved in NK cells, B lymphocyte, CD4 + , and CD8 + T lymphocyte stimulation [121]. Nonetheless, IFN-I signaling is suppressed and delayed in coronavirus infections, such as SARS-CoV, MERS-CoV, and SARS-CoV-2, causing viral clearance dysfunction. Persistent viral replication causes the release of uncontrolled proinflammatory cytokines, resulting in the cytokine storm [116,122,123]. This phenomenon is also called paradoxical hyperinflammation. Although IFN-II, consisting of IFN-γ, has an immune function overlapping with IFN-I to stimulate the antiviral and innate immune response, it causes enhancement predominantly in antigen-presenting cells (APCs) [124]. Several studies have revealed that IFN-II exhaustion is usually found in severe COVID-19 patients, suggesting a vital role of IFN-I in the immunopathology of the disease [125][126][127]. As a result, IFN signaling enrichment found in the peripheral white blood cells can be from the compensating mechanism of immune cells for IFN suppression and delay in severe cases. However, IFN signaling predominance in this study can indicate the ongoing activation. Persistent IFN stimulation induces apoptosis of CD4+ T lymphocytes and causes lymphopenia. It also increases proinflammatory cytokine production [128][129][130]. Other evidence supports the notion that increased IFN-I in influenza viral infection can release excessive proinflammatory cytokines, resulting in respiratory epithelial apoptosis and severe pneumonia [131]. Because the IFN function in severe COVID-19 is complicated, the study of IFN roles in the disease should be more investigated.
Leukocyte proliferation is the immune defense mechanism responding to infections. However, cell division and cell cycle regulation can also participate in the pathophysiology of COVID-19-associated cytokine storms. In the RNA-seq study used to construct the IPIN, the patient data revealed that severe COVID-19 cases had elevated neutrophils more than lymphocytes [87]. An excessive number of neutrophils can cause increased production of proinflammatory cytokines. Hence, finding the key immune-related proteins in the cell cycle regulation is necessary to modulate the immune response in severe cases. Nucleic acid metabolism was also found in the enrichment analysis. Increased cell division also stimulates nucleic acid metabolism to produce DNA materials for chromosome segregation. Although cell division is found in leukocytes during the infection, several studies revealed that SARS-CoV could induce cell cycle arrest in host cell lines by using viral proteins interacting with cyclin and the cyclin-dependent kinase (CDK) complex [132][133][134]. Nevertheless, there is still no research about the relation between SARS-CoV-2 and leukocyte cell cycle regulation. Therefore, studying the cell cycle regulation of host cells in COVID-19 requires further investigation.
Another term found in the enrichment analysis was apoptosis or programmed cell death. Apoptosis in COVID-19 causes deleterious effects, leading to severe complications. For instance, lymphocyte apoptosis results in lymphopenia and delays adaptive immune response, increasing the cytokine storm risk. Many studies have found that B and T lymphocyte apoptosis is associated with severe COVID-19 [135][136][137]. Apoptosis is also found in respiratory epithelial and endothelial cells, causing the blood-air barrier defect. This event leads to ARDS progression. Therefore, apoptosis regulation in severe COVID-19 is essential to improve adaptive immunity and reduce fatality rates.
Furthermore, there was protein processing found in the IPIN enrichment analysis. Protein processing is post-translational modification occurring in the endoplasmic reticulum (ER) and Golgi apparatus. The immune cells during the infection have more active functions, such as cell proliferation and cytokine production. Hence, protein processing is highly expressed to play a role in these activities.

Important Immune-Related Proteins with Node Centralities
Network centrality aspects such as degree, betweenness, closeness, and eigenvector centralities of each node in the IPIN are illustrated in Table S6 in Supplementary Materials. Nodes with degree, betweenness, closeness, or eigenvector centralities above the 90th percentile are shown in Tables S9-S12 in Supplementary Materials, respectively. From the tables, 15 nodes have large degree values, and 10 nodes have high betweenness values. In addition, 11 and 10 nodes have high closeness and eigenvector scores, respectively. These nodes are proteins playing a role in antiviral, innate immune, apoptosis, and cell cycle regulation signaling pathways. The function of each protein including high centrality predominance is displayed in Table 3.
From the tables, we found that nodes with high betweenness values were primarily involved in cell cycle regulation (6 of 10 nodes), while nodes with high degree and eigenvector values were mainly related to antiviral and innate immune signaling pathways (13 of 15 nodes and 9 of 10 nodes with high degree and eigenvector values, respectively). In addition, we discovered that the high closeness nodes had a similar proportion between immune signaling and cell proliferation, explaining that 6 of 11 nodes were involved in the innate immune response. At the same time, the rest of the nodes had a function associated with cell cycle regulation. The upset plot In Figure 7 shows that all nodes with high eigenvector values are in the high degree nodes. That means the high-eigenvector nodes form a subset of the high-degree nodes. Furthermore, the high-betweenness nodes share the node members mainly with the high-closeness nodes (the intersection size is calculated from 5 + 1 = 6). There is one node shared in all centralities (RSAD2). The degree, betweenness, and closeness centrality provide unique nodes that do not intersect with other centralities. The numbers of unique nodes in degree, betweenness, and closeness centrality are 3, 2, and 1, respectively.

Important Immune-Related Proteins with the Ranking Scores
As shown in Table 4, there are 10 nodes with a ranking score above the 90th percentile. Among them, eight nodes have functions relevant to innate immune response and antiviral activity: IFIT1, IFIT2, IFIT3, IRF7, ISG15, MX1, RSAD2, and STAT1. Interestingly, these immune nodes participate in IFN signaling pathways. The signaling pathway is usually activated when viral infections invade the hosts [138]. Stimulated IFNs increase the production of antiviral proteins, for example, ISG15 and MX1 [139,140]. In addition, the rest of the nodes, such as CDC25A and CCNA2, are involved in cell proliferation and cell cycle regulation. During infection, the immune system increases leukocyte proliferation to eradicate pathogens. As a result, cell division and cell cycle regulators can be found in the analysis. Both CDC25A and CCNA2 regulate cell cycle transition in G1/S phase [141]. Cell cycle control during G1/S phase is a critical point for cell division. Hence, drug repurposing targeting these regulators could improve the excessive leukocyte proliferation in the cytokine storm, leading to decreased morbidity and mortality rate in severe COVID-19.

Important Immune-Related Proteins with the Ranking Scores
As shown in Table 4, there are 10 nodes with a ranking score above the 90th percentile. Among them, eight nodes have functions relevant to innate immune response and antiviral activity: IFIT1, IFIT2, IFIT3, IRF7, ISG15, MX1, RSAD2, and STAT1. Interestingly, these immune nodes participate in IFN signaling pathways. The signaling pathway is usually activated when viral infections invade the hosts [138]. Stimulated IFNs increase the production of antiviral proteins, for example, ISG15 and MX1 [139,140]. In addition, the rest of the nodes, such as CDC25A and CCNA2, are involved in cell proliferation and To compare with the conventional methods usually used for identifying key proteins in IPINs such as degree and betweenness centrality, Figure 8 displays a Venn diagram of nodes found in the ranking score, degree, and betweenness centrality. The figure shows that the score covers mostly proteins in degree centrality (9 of 15 nodes in the degree set). Conversely, the score captured a few nodes in betweenness centrality (3 of 10 nodes in the betweenness set). Noticeably, the nodes merging between degree and betweenness set cover almost the high-value nodes analyzed from the four centralities (22 of 23 nodes). Thus, a combination of degree and betweenness centrality can provide the best result for identifying key proteins in IPINs. Although the ranking score could not capture other key proteins different from these two centralities, RSAD2 and IFIT1 (the top two genes in Table 4) were detected as the most important immune-related proteins. To compare with the conventional methods usually used for identifying ke teins in IPINs such as degree and betweenness centrality, Figure 8 displays a Ven gram of nodes found in the ranking score, degree, and betweenness centrality. T ure shows that the score covers mostly proteins in degree centrality (9 of 15 nodes degree set). Conversely, the score captured a few nodes in betweenness centralit 10 nodes in the betweenness set). Noticeably, the nodes merging between degree a tweenness set cover almost the high-value nodes analyzed from the four centralit of 23 nodes). Thus, a combination of degree and betweenness centrality can prov best result for identifying key proteins in IPINs. Although the ranking score cou capture other key proteins different from these two centralities, RSAD2 and IFIT top two genes in Table 4) were detected as the most important immune-related pro RSAD2 or viperin is a broad-spectrum antiviral protein in several viruses s measles, coxsackievirus A16, and enterovirus A71 [142][143][144]. Moreover, an animal revealed that RSAD2 was necessary for dendritic cell development [145]. The analysis of the overlapping genes in postmortem lung tissue from COVID-19 cas acute lung injury (ALI) murine model found that RSAD2 had a high degree centra a COVID-19-associated regulatory network [146]. Moreover, a lower respirator transcriptomic study revealed that RSAD2 expression correlated with the viral l mild and severe COVID-19 [147]. IFIT1 is an antiviral protein interacting with oth family proteins to form an IFN-dependent multiprotein complex. The complex pl important role to increase innate immunity against RNA viruses via binding be IFIT1 and 5′-triphosphate RNA (PPP-RNA) [148]. Nevertheless, an experimental f revealed that several SARS-CoV-2 proteins, such as nsp7, nsp15, M, 3CLpro, he and N proteins, suppressed IFT1 expression in HEK293T cells [149]. RSAD2 or viperin is a broad-spectrum antiviral protein in several viruses such as measles, coxsackievirus A16, and enterovirus A71 [142][143][144]. Moreover, an animal study revealed that RSAD2 was necessary for dendritic cell development [145]. The DEGs analysis of the overlapping genes in postmortem lung tissue from COVID-19 cases and acute lung injury (ALI) murine model found that RSAD2 had a high degree centrality in a COVID-19-associated regulatory network [146]. Moreover, a lower respiratory tract transcriptomic study revealed that RSAD2 expression correlated with the viral load in mild and severe COVID-19 [147]. IFIT1 is an antiviral protein interacting with other IFIT family proteins to form an IFN-dependent multiprotein complex. The complex plays an important role to increase innate immunity against RNA viruses via binding between IFIT1 and 5 -triphosphate RNA (PPP-RNA) [148]. Nevertheless, an experimental finding revealed that several SARS-CoV-2 proteins, such as nsp7, nsp15, M, 3CLpro, helicase, and N proteins, suppressed IFT1 expression in HEK293T cells [149].
Although the ranking score was computed to find the key proteins by considering the four centralities, it captured a few proteins when compared with the combination between degree and betweenness centrality. Moreover, the score covered the proteins mostly found in degree centrality. The reason to describe the result is that the proteins with high eigenvector centrality were the subset in degree centrality. Meanwhile, some proteins with high closeness centrality were found in both eigenvector and degree centrality. In addition, proteins found in betweenness centrality rarely overlapped with degree centrality. Therefore, the score calculation is weighted to degree centrality more than betweenness centrality. Noticeably, the combination of degree and betweenness centrality covered almost high-value centrality proteins rather than the ranking score, although it lost one protein with a high closeness value (CCNE1).
As described earlier, the 23 key proteins identified using centrality measurement were involved in innate immune response, cell cycle regulation, and apoptosis. IFI35 is an IFN signaling regulator response to viral infections. Many studies have found that IFI35 plays an essential role in cytokine storms and severity in COVID-19 and influenza [150][151][152][153]. IFIH1 is an intracellular viral sensing protein that stimulates the IFN signaling pathway when viral particles are detected in the cell [154]. IFIH1 was reported to participate in SARS-CoV-2 sensing and was associated with proinflammatory cytokine overproduction in COVID-19 [155][156][157]. DDX58 also plays a role as a cytoplasmic viral sensor. High DDX58 expression in COVID-19 was associated with cytokine responses [158]. IFI6, IFIT1-3, IRF7, ISG15, MX1, OAS1, OAS2, OASL, and RSAD2 have antiviral activity functions. Many studies in IFI6 have revealed that IFI6 plays an essential role against hepatitis B virus (HBV) and flavivirus replication [159,160]. Several systems biological and transcriptomic studies have also shown that IFI6 is a hub gene in the gene co-expression networks and transcriptomic profiles in COVID-19 [52,[161][162][163]. IFIT1-3, OAS1-3, and OASL were upregulated in SARS-CoV-2 and other coronavirus infections [164][165][166][167]. Furthermore, the IRF7 mutation causing loss of function was reported to be associated with severe pneumonia progression in COVID-19 [168,169]. The infected macrophage cell line study showed that extracellular ISG15 stimulated proinflammatory cytokine production, leading to hyperinflammation [170]. The result of a COVID-19 case-control study revealed that MX1 expression was increased depending on elevated viral load, and the expression was decreased in elderly patients [171]. Older patients have a high risk for COVID-19-associated cytokine storms, suggesting that low MX1 expression could play a vital role in the cytokine storm. STAT1 is a signal transduction protein related to various signaling pathways such as IFN, IL-6, epidermal growth factor (EGF), and platelet-derived growth factor (PDGF) pathways [172][173][174]. Several studies have indicated that phosphorylated STAT1 increases in severe COVID-19 patients, causing STAT1 signaling dysfunction and failed IFN activation [175][176][177]. XAF1 is a tumor suppressor protein playing as a positive feedback regulator in the p53-induced apoptotic signaling pathway [178]. Numerous studies have reported that XAF1 dysfunction plays a vital role in tumor progression [179][180][181][182]. A single-cell transcriptomic study in peripheral blood mononuclear cells showed that COVID-19 caused XAF1-induced T lymphocyte apoptosis, leading to adaptive immune impairment [183]. In addition, IPIN analysis from COVID-19 patient lung tissue revealed that IFIH1, DDX58, ISG15, OASL, and XAF1 were hub genes in the network [184].
CCNA2 and CCNE1 are CDK kinase regulators during G1/S and G2/M phases in the cell cycle. Numerous studies have indicated CCNA2 and CCNE1 play a central role in various types of malignancy such as hepatocellular carcinoma, breast, and colon cancer [185][186][187][188]. CDC25A is a protein required in the cell cycle by activating CDKs [189]. CDC25A overexpression was found in head and neck, breast, ovarian, and non-small cell lung cancer [190][191][192][193]. An immune study revealed that CDC25A had activities decreasing IFN-β transcription and DDX58-mediated antiviral signaling pathway [194]. CDC20 has a function involved in chromosome segregation and is the target for spindle assembly checkpoint (SAC) [195,196]. High-expressed CDC20 was related to the worst prognosis in lung squamous cell carcinoma [197]. An IPIN analysis in COVID-19-induced thrombocytopenia also reported that CDC20 was highly expressed in COVID-19 with thrombocytopenia [198]. CMPK2 is an enzyme associated with the nucleotide salvage pathway. Many studies have shown that CMPK2 participates in IFN-I activation and antiviral immune response [199][200][201]. In COVID-19 studies, CMPK2 was highly upregulated in severe cases related to ARDS [202,203]. Moreover, FOXM1-dependent tissue regeneration is impaired in severe COVID-19 cases, causing sustained lung injury and a high fatality rate [204]. RRM2, a cell cycle regulator, had an increased expression in lung adenocarcinoma with a poor prognosis [205]. A gene co-expression network and functional enrichment study also showed that RRM2 was a component in a module involved in p53 signaling pathway, a cell cycle and apoptosis pathway, in COVID-19 [206].

Potential Drugs to Cure Severe COVID-19 Patients
We used the key immune-related proteins from the ranking score, degree, and betweenness centrality to find the drug-gene or drug-protein interactions from four well-known public databases: DrugBank database [110], TTD [111], CTD [112], and GeneCards [113]. In addition, the protein not found in both degree and betweenness centrality (CCNE1) was used to discover the interaction. The result from database searching revealed 115 FDAapproved drugs interacting with the key genes or proteins. Table S13 in Supplementary Materials shows the drug-gene and drug-protein interactions in detail among the 23 key immune-related proteins.
STITCH v5.0 [114], a drug-protein interaction database, was conducted to confirm the result from the databases by the confidence score cut-off value of 0.9. The STITCH result is demonstrated in Figure 9. The key immune-related proteins are classified into two groups. The former is an antiviral, innate immune response, and apoptosis signaling pathway group, and the latter is a cell division and cell cycle group. Figure 9a displays the drug-protein interaction network in the innate immune response, and Figure 9b illustrates the drug-protein interaction network involved in cell cycle regulation. There are seven candidate drugs interacting with these seven key proteins. Three drugs are associated with the key proteins related to innate immune response and apoptosis. In contrast, the rest interacts with the proteins involved in cell cycle regulation. In the innate immune and apoptosis network, polyinosinic:polycytidylic acid (poly(I:C)) interacts with IFIH1 and DDX58, while mitomycin C interacts with MX1. Decitabine also binds to XAF1 in the network. There are four drug-protein interactions in the cell cycle regulation network. RRM2 is interacted with gemcitabine and hydroxyurea, respectively. Tamoxifen binds to FOXM1, and curcumin has an interaction with CCNE1.
From the STITCH v5.0 result, the seven candidate drugs interacted with the seven key proteins. For instance, Poly(I:C) interacted with IFIH1 and DDX58. Poly(I:C) is an immune stimulant used to activate innate immunity such as IFN by the TL3 agonist effect. It also induces cancer cell apoptosis in various types of malignancy: cervical, prostate, and colon cancer [207][208][209][210][211]. Poly(I:C) also increases cytotoxic activity in CD4+ T lymphocytes in viral infections, promoting adaptive immune response [212]. A study in influenza A virus (IAV) and a SARS-CoV-infected mice model revealed that poly(I:C) had a protective effect in fatal respiratory infections [213]. Interestingly, intranasal poly(I:C) in mice with SARS-CoV-2 infection showed a decreased viral load, suggesting that poly(I:C) can be an effective drug for treating the disease [214]. However, some studies reported that poly(I:C) increased proinflammatory production [215][216][217]. Therefore, further studies about poly(I:C) in COVID-19 treatment should consider the drug dosage and administration route to prevent the cytokine storm due to the medication. Mitomycin C, a chemotherapeutic agent, interacted with MX1. It is used to treat many cancer types [218][219][220][221]. A systems biological study revealed that mitomycin C interacted with MX1, a key protein in an IPIN, suggesting further studies in the role of the drug in antiviral stimulation [50]. In addition, there was a drug-protein interaction between decitabine and XAF1. Decitabine is a pyrimidine nucleoside antimetabolite used to treat myelodysplastic syndrome (MDS) and acute myeloid leukemia (AML) [222,223]. A study in the mice model showed that decitabine improved FOXM1-dependent endothelial regeneration and vascular repair [204]. As mentioned above, lung tissue degeneration can cause worsening lung injury. Decitabine then could play a role in decreasing lung injury in severe COVID-19. Interestingly, there is a clinical trial that has been studying decitabine treatment in critical ill COVID-19 patients. The estimated research completion date August 2022. Gemcitabine and hydroxyurea also interacted with RRM2. Gemcitabine, a pyrimidine nucleoside analog and chemotherapeutic agent, is used to treat solid tumors such as bladder, pancreatic, breast, and non-small cell lung cancer [224]. Several studies in cell lines have found that decitabine decreases SARS-CoV-2 replication [225][226][227]. In addition, in a cohort study, gemcitabine reduced SARS-CoV-2 infection in cancer patients [228]. Hydroxyurea, an antimetabolite treating sickle cell anemia, has anti-inflammatory and immunomodulatory effects and was expected to apply well in COVID-19 [229]. From the STITCH v5.0 result, the seven candidate drugs interacted with the seven key proteins. For instance, Poly(I:C) interacted with IFIH1 and DDX58. Poly(I:C) is an immune stimulant used to activate innate immunity such as IFN by the TL3 agonist effect. It also induces cancer cell apoptosis in various types of malignancy: cervical, prostate, and colon cancer [207][208][209][210][211]. Poly(I:C) also increases cytotoxic activity in CD4+ T lymphocytes in viral infections, promoting adaptive immune response [212]. A study in influenza A virus (IAV) and a SARS-CoV-infected mice model revealed that poly(I:C) had a protective effect in fatal respiratory infections [213]. Interestingly, intranasal poly(I:C) in mice with SARS-CoV-2 infection showed a decreased viral load, suggesting that poly(I:C) can be an effective drug for treating the disease [214]. However, some studies reported that poly(I:C) increased proinflammatory production [215][216][217]. Therefore, further studies about poly(I:C) in COVID-19 treatment should consider the drug dosage and administration route to prevent the cytokine storm due to the medication. Mitomycin C, a chemotherapeutic agent, interacted with MX1. It is used to treat many cancer types [218][219][220][221]. A systems biological study revealed that mitomycin C interacted with MX1, a key protein in an IPIN, suggesting further studies in the role of the drug in antiviral stimulation [50]. In addition, there was a drug-protein interaction between decitabine and XAF1. Decitabine is a pyrimidine nucleoside antimetabolite used to treat myelodysplastic syndrome (MDS) and acute myeloid leukemia (AML) [222,223]. A study in the mice model showed that decitabine improved FOXM1-dependent endothelial regeneration and vascular repair [204]. As mentioned above, lung tissue degeneration can cause worsening lung injury. Decitabine then could play a role in decreasing lung injury in severe COVID-19. Interestingly, there is a clinical trial that has been studying decitabine treatment in critical ill COVID-19 patients. The estimated research completion date August 2022. Gemcitabine and hydroxyurea also interacted with RRM2. Gemcitabine, a pyrimidine nucleoside analog and chemotherapeutic agent, is used to A clinical study also revealed that the mortality rate was low in COVID-19 patients receiving hydroxyurea, suggesting a vital role of hydroxyurea in COVID-19 treatment [230]. Tamoxifen, a selective estrogen receptor modulator, had an interaction with CCNE1. The drug is used to treat estrogen receptor-positive breast cancer [231]. Tamoxifen downregulates TMPRSS2, preventing SARS-CoV-2 entry into host cells [232]. A preclinical study showed that tamoxifen reduced SARS-CoV-2 in vitro and in vivo [233]. Moreover, a review article in clinical studies found that tamoxifen decreased COVID-19 susceptibility rates in breast cancer patients [234]. Tamoxifen also inhibited viral replication and virus entry in many virus types such as EBOLA, MERS, and SARS-CoV-2 [234]. The drug repurposing result also revealed that curcumin interacted with CCNE1. It is worth noting that curcumin is a promiscuous molecule acting on many receptors. Hence, the effect of curcumin in COVID-19 can be from other mechanisms. Our study only proposed one of the possible mechanisms of curcumin in severe COVID-19 treatment. Curcumin is a natural product found in turmeric (Curcuma longa). Many studies have indicated that curcumin has anti-inflammatory and antioxidant effects [235][236][237]. Moreover, it provides effects in apoptosis promotion, cell proliferation inhibition, anti-cell adhesion and invasion, decreased angiogenesis, and anti-microbial activity [238]. Therefore, its clinical application is related to numerous diseases, such as rheumatoid arthritis, inflammatory bowel diseases, osteoarthritis, and various types of cancer. In COVID-19, several review articles have revealed that curcumin inhibits viral entry and replication [239][240][241]. It also promotes IFN and antiviral signaling pathway and decreases proinflammatory cytokine production. Curcumin has protective effects on ARDS by reducing NF-κB, inflammasome, and IL-8 pathway. Furthermore, a randomized control trial study showed that mild, moderate, and severe COVID-19 patients taking oral curcumin with piperine had better clinical outcomes and lower hospitalization duration than the controls [242]. A systematic review reported that curcumin reduced the proinflammatory cytokines, such as IL1β and IL6. It also increased the anti-inflammatory cytokines, for example, IL-10, IL-35 and transforming growth factor α (TGF-α) [243]. Therefore, further study should focus on the effective dose and administration routes of curcumin.

Discussion
This study constructed the immune-related protein interactions network, IPIN, for severe COVID-19 based on the leukocyte transcriptomic profile of critically ill patients using network propagation on the human interactome network. Functional enrichment analysis and network clustering were operated to discover the underlying molecular mechanisms of immune-related severe COVID-19. Topological analysis, centrality, and ranking score measurements were calculated to identify the key immune-related proteins. Finally, the drug-protein interactions were searched to find the candidate drugs to treat the severe COVID-19 patients.
Diffusion-based IPIN construction and permutation testing provided the highly connected immune-related proteins in IPIN. IPIN was a network with a small-world effect in relation to other biological networks. The small-world effect was proved by the low average shortest path length and high average clustering coefficient. The scale-free property cannot be explained in the network due to a lack of a relationship between degree and its probability. However, the IPIN was less likely to be the random network because it provided the significant PPI enrichment p-value from the STRING database. We performed the four network centralities (degree, betweenness, closeness, and eigenvector) and ranking score to find the key immune-related proteins. The results showed that the combination of degree and betweenness centrality covered a wide range of the key proteins. However, the ranking score can detect the main key proteins: RSAD2 and IFTI1.
We identified 23 key immune-related proteins, such as CCNA2, CCNE1, CDC20,  CDC25A, CMPK2, DDX58, FOXM1, IFI6, IFI35, IFIH1, IFIT1, IFIT2, IFIT3, IRF7, ISG15,  MX1, OAS1, OAS2, OASL, RRM2, RSAD2, STAT1, and XAF1, using the four centralities and ranking score measurement. These proteins all play an important role in severe COVID-19, validated by several computational, experimental, and clinical studies. The functional enrichment analysis from the whole network and the modules obtained from both MCODE and MCL methods produced similar results. The enrichment terms were divided into four main categories: cell cycle regulation, antiviral and innate immune response, apoptosis, and nucleotide metabolism. These terms were in accordance with leukocytes during viral infections. Furthermore, the main terms were accepted with the functional classification in the key proteins. For instance, the key proteins involved in cell cycle regulation consisted of CCNA2, CCNE1, CDC20, CDC25A, and RRM2, while the others associated with antiviral and innate immune response were DDX58, FOXM1, IFI6, IFI35, IFIH1, IFIT1, IFIT2, IFIT3, IRF7, ISG15, MX1, OAS1, OAS2, OASL, RSAD2, and STAT1. In addition, the remaining key proteins such as CMPK2 and XAF1 play a crucial role in nucleotide metabolism and apoptosis, respectively.
Drug repurposing based on drug-gene and drug-protein interaction database searching provided the seven potential candidate drugs, poly(I:C), mitomycin C, decitabine, gemcitabine, hydroxyurea, tamoxifen, and curcumin. There were three drugs interacting with the key proteins related to antiviral and innate immune response: poly(I:C), mitomycin C, and decitabine. Other drugs interacted with the key proteins involved in cell cycle regulation and apoptosis. Among the candidate drugs, we recommend that the drugs interacting with proteins involved in IFN and antiviral signaling should be used carefully in clinical application because they promote IFN stimulation. IFN overactivation can result in excessive proinflammatory cytokine production in some studies that we mentioned [128][129][130][131]. Moreover, chemotherapeutic agents cause many adverse side effects and should be performed as the second choice. Therefore, we suggest using curcumin and tamoxifen as the first choices for clinical application. They have fewer side effects than other chemotherapeutic agents because curcumin is a natural product and tamoxifen is a targeted drug. In addition, both curcumin and tamoxifen have several clinical studies to support their effectiveness in COVID-19 treatment.

Limitations and Future Study Suggestions
Although this study provides novel knowledge and candidate targeted drugs in COVID-19, it has some limitations to be explained. First, the network diffusion method is an algorithm consuming computational time and memory space which depends on the number of nodes, interactions, and permutation tests. We conducted high-performance computing (HPC) for running the LHD algorithm in the original and 1000 random sets to perform the permutation test. Second, many proinflammatory cytokine signaling pathways, such as IL-1β, IL-6, IL-12, IL-18, IL-33, and TNF-α, were rarely detected in this IPIN analysis. Our reason for explaining the issue is that the data came from comparing controls and COVID-19 cases in an intensive care unit. Typically, critically ill patients have stress and inflammatory responses, leading to increased proinflammatory cytokines. Therefore, there was no difference in the proinflammatory cytokine gene expression between cases and controls. Moreover, proinflammatory cytokines are usually released from respiratory epithelial and immune cells in the lung parenchymal tissue. Studies in peripheral white blood cells can lose this proinflammatory cytokine information. Our suggestion for future research is to perform lung transcriptomic profiles comparing severe COVID-19 cases and mild illness or healthy cases for IPIN construction and analysis. Furthermore, a single-cell approach should be conducted to identify an IPIN in each cell type. In COVID-19, there are differences between each cell type such as cell count, behavior, function, and pathogenesis. Therefore, identifying key proteins in these cells can help to treat the disease more precisely.

Conclusions
This study proposed LHD algorithms to perform network diffusion on the human interactome network to construct the immune-related IPIN in severe COVID-19 based on the transcriptomic data. Functional enrichment analysis found that the network contained the proteins involved in antiviral and innate immune response signaling pathways, cell cycle regulation, apoptosis, and protein processing. The degree and betweenness centrality combination covered almost the key proteins from the four centrality measurements. These key proteins play a significant role in cell proliferation, antiviral activity, and innate immunity responding to the SARS-CoV-2 infection. Moreover, the candidate drugs targeting the key proteins were found from database searching. Most of them have experimental data supporting their effectiveness in COVID-19 treatment. However, other key proteins and candidate drugs were not found in our method and need further investigation. Therefore, a combination of advanced experimental and computational tools should be conducted for further efficient treatment discovery relating to COVID-19.
Supplementary Materials: The following supporting information can be found and downloaded at: https://www.mdpi.com/article/10.3390/biom12050690/s1. Table S1: 224 differentially expressed genes (DEGs) (FDR<0.05 and log2 FC > 1.5); Table S2: The DEGs in the protein-protein interaction network from STRING with Ensembl protein IDs; Table S3: The immune-related proteins from DEGs in the human interactome data as seed nodes for network diffusion; Table S4: Diffusion scores of the original set and statistic data of each node in the whole PPI network; Table S5: 154 protein nodes with significant diffusion scores by the permutation test; Table S6: Degree, betweenness, closeness, and eigenvector centrality values of each node in the IPIN; Table S7: Network clustering of the IPIN by MCL algorithm provided 4 modules: MCL1, MCL2, MCL3, and MCL4. Table S8: The detail of enrichment results in each MCL module; Table S9: List of nodes with high degree centrality; Table S10: List of nodes with high betweenness centrality; Table S11: List of nodes with high closeness centrality; Table S12: List of nodes with high eigenvector centrality; Table S13: Drug-gene interactions of the identified key immune-related proteins; Table S14: Edgelist of the IPIN.