key: cord-0823743-7mdxkax4 authors: Sagulkoo, Pakorn; Suratanee, Apichat; Plaimas, Kitiporn title: Immune-Related Protein Interaction Network in Severe COVID-19 Patients toward the Identification of Key Proteins and Drug Repurposing date: 2022-05-11 journal: Biomolecules DOI: 10.3390/biom12050690 sha: 31fef3695745a37781f8f84a1d671a98fd65c6e2 doc_id: 823743 cord_uid: 7mdxkax4 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. 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. 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. 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. 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). 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. 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. 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. 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 (), 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. 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. 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. 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. 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. 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 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. 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. 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. 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] . 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. 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. 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. 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. World Health Organization (WHO) Coronavirus (COVID-19) Dashboard. 2022 Severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) and coronavirus disease-2019 (COVID-19): The epidemic and the challenges The species Severe acute respiratory syndrome-related coronavirus: Classifying 2019-nCoV and naming it SARS-CoV-2 SARS: Clinical features and diagnosis Severe acute respiratory syndrome (SARS): A review Middle East Respiratory Syndrome (MERS) Host Determinants of MERS-CoV Transmission and Pathogenesis The Middle East Respiratory Syndrome Coronavirus-A Continuing Risk to Global Health Security Middle East Respiratory Syndrome Author Correction: A new coronavirus associated with human respiratory disease in China Addendum: A pneumonia outbreak associated with a new coronavirus of probable bat origin The first 12 months of COVID-19: A timeline of immunological insights International Public Health Responses to COVID-19 Outbreak: A Rapid Review. Iran Comparative sequence analysis of SARS-CoV-2 suggests its high transmissibility and pathogenicity The genetic sequence, origin, and diagnosis of SARS-CoV-2 Genetics and genomics of SARS-CoV-2: A review of the literature with the special focus on genetic diversity and SARS-CoV-2 genome detection Insights into SARS-CoV-2 genome, structure, evolution, pathogenesis and therapies: Structural genomics approach Molecular mechanism of interaction between SARS-CoV-2 and host cells and interventional therapy Role of Structural and Non-Structural Proteins and Therapeutic Targets of SARS-CoV-2 for COVID-19 SARS-CoV-2 Accessory Proteins in Viral Pathogenesis: Knowns and Unknowns Contributions of human ACE2 and TMPRSS2 in determining host-pathogen interaction of COVID-19 Novel Small-Molecule Inhibitors of the SARS-CoV-2 Spike Protein Binding to Neuropilin 1 Can SARS-CoV-2 Virus Use Multiple Receptors to Enter Host Cells? Coronavirus biology and replication: Implications for SARS-CoV-2 COVID-19 and the human innate immune system SARS-CoV-2 infection: The role of cytokines in COVID-19 disease Role of Inflammatory Cytokines in COVID-19 Patients: A Review on Molecular Mechanisms, Immune Functions, Immunopathology and Immunomodulatory Drugs to Counter Cytokine Storm. Vaccines 2021 The triumvirate of NF-κB, inflammation and cytokine storm in COVID-19 Transmission, Diagnosis, and Treatment of Coronavirus Disease 2019 (COVID-19): A Review COVID-19-associated acute respiratory distress syndrome (CARDS): Current knowledge on pathophysiology and ICU treatment-A narrative review Cytokine Storm in COVID-19-Immunopathological Mechanisms, Clinical Considerations, and Therapeutic Approaches: The REPROGRAM Consortium Position Paper. Front Severe COVID-19 Thromboembolic Complications in Severe COVID-19: Current Antithrombotic Strategies and Future Perspectives COVID-19-associated coagulopathy and disseminated intravascular coagulation COVID-19 and multiorgan failure: A narrative review on potential mechanisms Macrophage activation syndrome and COVID-19 Efficacy and safety of COVID-19 vaccines: A systematic review COVID-19 vaccine: A recent update in pipeline vaccines, their design and development strategies An update on drugs with therapeutic potential for SARS-CoV-2 (COVID-19) treatment Molnupiravir: A new candidate for COVID-19 treatment Corticosteroid use in COVID-19 patients: A systematic review and meta-analysis on clinical outcomes Systemic corticosteroids for the treatment of COVID-19 Corticosteroid therapy for COVID-19: A systematic review and meta-analysis of randomized controlled trials Efficacy of corticosteroid in patients with COVID-19: A multi-center retrospective study and meta-analysis A Review of the Evidence for Corticosteroids in COVID-19 The risks of systemic corticosteroid use A comprehensive review of the adverse effects of systemic corticosteroids Increased Secondary Infection in COVID-19 Patients Treated with Steroids in New York City Targeting hub genes and pathways of innate immune response in COVID-19: A network biology perspective Repurposing novel therapeutic candidate drugs for coronavirus disease-19 based on protein-protein interaction network analysis Identification of hub genes and molecular subtypes in COVID-19 based on WGCNA Identification of Hub Genes and Key Pathways Associated with Anti-VEGF Resistant Glioblastoma Using Gene Expression Data Analysis Network-based drug repurposing for novel coronavirus 2019-nCoV/SARS-CoV-2 Network analysis and molecular mapping for SARS-CoV-2 to reveal drug targets and repurposing of clinically developed drugs Multi-Data Aspects of Protein Similarity with a Learning Technique to Identify Drug-Disease Associations Meta-Path Based Gene Ontology Profiles for Predicting Drug-Disease Associations Repurposing of Drugs for SARS-CoV-2 Using Inverse Docking Fingerprints Identification of inflammatory bowel disease-related proteins using a reverse k-nearest neighbor search Hybrid Deep Learning Based on a Heterogeneous Network Profile for Functional Annotations of Plasmodium falciparum Genes Protein network analysis to prioritize key genes in amyotrophic lateral sclerosis Heterogeneous Network Model to Identify Potential Associations Between Plasmodium vivax and Human Proteins Prediction of Human-Plasmodium vivax Protein Associations From Heterogeneous Network Structures Based on Machine-Learning Approach NPF:network propagation for protein function prediction Reverse Nearest Neighbor Search on a Protein-Protein Interaction Network to Infer Protein-Disease Associations Network-based association analysis to infer new disease-gene relationships using large-scale protein interactions Identifying potential association on gene-disease network via dual hypergraph regularized least squares DDA: A Novel Network-Based Scoring Method to Identify Disease-Disease Associations Network-based characterization of disease-disease relationships in terms of drugs and therapeutic targets Network-Based Target Prioritization and Drug Candidate Identification for Multiple Sclerosis: From Analyzing "Omics Data" to Druggability Simulations Integrative Systems Biology Approaches to Identify Potential Biomarkers and Pathways of Cervical Cancer Network diffusion with centrality measures to identify disease-related genes Construction and Analysis of Protein-Protein Interaction Network to Identify the Molecular Mechanism in Laryngeal Cancer Construction and analysis of protein-protein interaction networks based on proteomics data of prostate cancer The Protein-Protein Interaction Network of Litopenaeus vannamei Haemocytes Construction and analysis of the protein-protein interaction network related to essential hypertension Construction and Analysis of Protein-Protein Interaction Network of Heroin Use Disorder Cascading from SARS-CoV-2 to Parkinson's Disease through Protein-Protein Interactions Network-based transcriptomic analysis identifies the genetic effect of COVID-19 to chronic kidney disease patients: A bioinformatics approach Detection of spreader nodes in human-SARS-CoV protein-protein interaction network Shortest-path network analysis is a useful approach toward identifying genetic determinants of longevity Network propagation: A universal amplifier of genetic associations STRING v11: Protein-protein association networks with increased coverage, supporting functional discovery in genome-wide experimental datasets Metascape provides a biologist-oriented resource for the analysis of systems-level datasets Dplyr: A Grammar of Data Manipulation The igraph software package for complex network research Transcriptional profiling of leukocytes in critically ill COVID-19 patients: Implications for interferon response and coagulation Use of Laplacian Heat Diffusion Algorithm to Infer Novel Genes With Functions Related to Uveitis Identification of proliferative diabetic retinopathy-associated genes on the protein-protein interaction network by using heat diffusion algorithm Network Diffusion Algorithms UpSetR: An R package for the visualization of intersecting sets and their properties Centrality analysis methods for biological networks and their application to gene regulatory networks The organisational structure of protein networks: Revisiting the centrality-lethality hypothesis Betweenness centrality of intracranial electroencephalography networks and surgical epilepsy outcome Variability of Betweenness Centrality and Its Effect on Identifying Essential Genes Centers of complex networks The connectivity structure, giant strong component and centrality of metabolic networks Eigenvector centrality mapping for analyzing connectivity patterns in fMRI data of the human brain Brain network alterations in Alzheimer's disease measured by eigenvector centrality in fMRI are related to cognition and CSF biomarkers The Association of Glucose Metabolism and Eigenvector Centrality in Alzheimer's Disease Altered eigenvector centrality is related to local resting-state network functional connectivity in patients with longstanding type 1 diabetes mellitus. Hum. Brain Mapp The Gene Ontology (GO) database and informatics resource New perspectives on genomes, pathways, diseases and drugs The reactome pathway knowledgebase Connecting communities Canonical pathways, networks and transcriptional factor regulation by clinical strains of Mycobacterium tuberculosis in pulmonary alveolar epithelial cells CORUM: The comprehensive resource of mammalian protein complexes-2019 An automated method for finding molecular complexes in large protein interaction networks Using MCL to extract clusters from networks Therapeutic target database 2020: Enriched resource for facilitating research and early development of targeted therapeutics Comparative Toxicogenomics Database (CTD): Update 2021 GeneCards Version 3: The human gene integrator Augmenting protein-chemical interaction networks with tissue and affinity data Network biology: Understanding the cell's functional organization Type I interferons in viral control and immune regulation Type I interferons: Expression and signalization Type I interferons in host defense Conditional expression of type I interferon-induced bovine Mx1 GTPase in a stable transgenic vero cell line interferes with replication of vesicular stomatitis virus Interferon-inducible protein Mx1 inhibits influenza virus by interfering with functional viral ribonucleopr.rotein complex assembly Type I interferons in infectious disease The Role of Type I Interferons in the Pathogenesis and Treatment of COVID-19 The signal pathways and treatment of cytokine storm in COVID-19 The Dual Nature of Type I and Type II Interferons COVID-19 in Children: Expressions of Type I/II/III Interferons, TRIM28, SETDB1, and Endogenous Retroviruses in Mild and Severe Cases Type I, II, and III interferon signatures correspond to COVID-19 disease severity Functional Exhaustion of Type I and II Interferons Production in Severe COVID-19 Patients Interferons alpha and beta as immune regulators-A new look Pathogenic potential of interferon αβ in acute influenza infection Chronic exposure to type-I IFN under lymphopenic conditions alters CD4 T cell homeostasis Lung epithelial apoptosis in influenza virus pneumonia: The role of macrophage-expressed TNF-related apoptosis-inducing ligand G1 phase cell cycle arrest induced by SARS-CoV 3a protein via the cyclin D3/pRb pathway SARS coronavirus 7a protein blocks cell cycle progression at G0/G1 phase via the cyclin D3/pRb pathway The nucleocapsid protein of severe acute respiratory syndrome-coronavirus inhibits the activity of cyclin-cyclin-dependent kinase complex and blocks S phase progression in mammalian cells Apoptosis-induced T-cell lymphopenia is related to COVID-19 severity T cell apoptosis characterizes severe Covid-19 disease Distinct Mitochondria-Mediated T-Cell Apoptosis Responses in Children and Adults With Coronavirus Disease Interferon signaling pathways ISG15 in antiviral immunity and beyond Interferon-induced Mx proteins in antiviral host defense Control of cell cycle transcription during G1 and S phases A unifying view of the broad-spectrum antiviral activity of RSAD2 (viperin) based on its radical-SAM chemistry A key anti-viral protein, RSAD2/VIPERIN, restricts the release of measles virus from infected cells RSAD2 and AIM2 Modulate Coxsackievirus A16 and Enterovirus A71 Replication in Neuronal Cells in Different Ways That May Be Associated with Their 5' Nontranslated Regions Rsad2 is necessary for mouse dendritic cell maturation via the IRF7-mediated signaling pathway Core genes involved in the regulation of acute lung injury and their association with COVID-19 and tumor progression: A bioinformatics and experimental study COVID-19 ARDS is characterized by a dysregulated host response that differs from cytokine storm and is modified by dexamethasone IFIT1 is an antiviral protein that recognizes 5'-triphosphate RNA Severe Acute Respiratory Syndrome Coronavirus 2 (SARS-CoV-2) Membrane (M) and Spike (S) Proteins Antagonize Host Type I Interferon Response IFP35 as a promising biomarker and therapeutic target for the syndromes induced by SARS-CoV-2 or influenza virus Interferon Inducer IFI35 regulates RIG-I-mediated innate antiviral response through mutual antagonism with Influenza protein NS1 Temporal dynamics of the host molecular responses underlying severe COVID-19 progression and disease resolution Effects of IFIH1 rs1990760 variants on systemic inflammation and outcome in critically ill COVID-19 patients in an observational translational study The role of interferon induced with helicase C domain 1 (IFIH1) in the development of type 1 diabetes mellitus SARS-CoV-2 sensing by RIG-I and MDA5 links epithelial infection to macrophage inflammation The RNA sensor MDA5 detects SARS-CoV-2 infection The transcriptomic profiling of SARS-CoV-2 compared to SARS, MERS, EBOV, and H1N1 A CRISPR screen identifies IFI6 as an ER-resident interferon effector that blocks flavivirus replication The Functional and Antiviral Activity of Interferon Alpha-Inducible IFI6 Against Hepatitis B Virus Replication and Gene Expression Host transcriptome-guided drug repurposing for COVID-19 treatment: A meta-analysis based approach Single-Cell Transcriptome Analysis Highlights a Role for Neutrophils and Inflammatory Macrophages in the Pathogenesis of Severe COVID-19 Weighted Gene Co-Expression Network Analysis Combined with Machine Learning Validation to Identify Key Modules and Hub Genes Associated with SARS-CoV-2 Infection In vivo antiviral host transcriptional response to SARS-CoV-2 by viral load, sex, and age IFITM proteins promote SARS-CoV-2 infection and are targets for virus inhibition in vitro Discriminating mild from critical COVID-19 by innate and adaptive immune single-cell profiling of bronchoalveolar lavages Inborn errors of type I IFN immunity in patients with life-threatening COVID-19 Role of Host Immune and Inflammatory Responses in COVID-19 Cases with Underlying Primary Immunodeficiency: A Review Altered ISGylation drives aberrant macrophage-dependent immune responses during SARS-CoV-2 infection SARS-CoV-2 Infection Boosts MX1 Antiviral Effector in COVID-19 Patients Stat1-dependent and -independent pathways in IFN-gamma-dependent signaling Transcription protein STAT1: Biology and relation to cancer Fagard, R. STAT1 and pathogens, not a friendly relationship A Putative Role of de-Mono-ADP-Ribosylation of STAT1 by the SARS-CoV-2 Nsp3 Protein in the Cytokine Storm Syndrome of COVID-19 An aberrant STAT pathway is central to COVID-19 Altered increase in STAT1 expression and phosphorylation in severe COVID-19 XAF1 directs apoptotic switch of p53 signaling through activation of HIPK2 and ZNF313 Down-regulation of the pro-apoptotic XIAP associated factor-1 (XAF1) during progression of clear-cell renal cancer Tumor suppressor XAF1 induces apoptosis, inhibits angiogenesis and inhibits tumor growth in hepatocellular carcinoma XAF1 forms a positive feedback loop with IRF-1 to drive apoptotic stress response and suppress tumorigenesis Combined effects of XAF1 and TRAIL on the apoptosis of lung adenocarcinoma cells Single-Cell Sequencing of Peripheral Mononuclear Cells Reveals Distinct Immune Response Landscapes of COVID-19 and Influenza Patients An integrative bioinformatics analysis for identifying hub genes associated with infection of lung samples in patients infected with SARS-CoV-2 Induction of cell cycle arrest via the p21, p27-cyclin E,A/Cdk2 pathway in SMMC-7721 hepatoma cells by clioquinol Cyclin A2/E1 activation defines a hepatocellular carcinoma subclass with a rearrangement signature of replication stress Expression and prognostic value of CDK1, CCNA2, and CCNB1 gene clusters in human breast cancer Comprehensive Analysis of Cyclin Family Gene Expression in Colon Cancer The role of Cdc25A in the regulation of cell proliferation and apoptosis Overexpression of CDC25A and CDC25B in head and neck cancers Overexpression of cdc25A and cdc25B is frequent in primary non-small cell lung cancer but is not associated with overexpression of c-myc Role of the Cdc25A phosphatase in human breast cancer Cell cyclerelated phosphatases CDC25A and B expression correlates with survival in ovarian cancer patients Phosphatase Cdc25A Negatively Regulates the Antiviral Immune Response by Inhibiting TBK1 Activity The regulation of Cdc20 proteolysis reveals a role for APC components Cdc23 and Cdc27 during S phase and early mitosis At the Crossroads between Chromosome Segregation and Mitotic Exit High Expression Levels of CDK1 and CDC20 in Patients with Lung Squamous Cell Carcinoma are Associated With Worse Prognosis Thrombocytopenia in COVID-19 and vaccineinduced thrombotic thrombocytopenia CMPK2 and BCL-G are associated with type 1 interferon-induced HIV restriction in humans Mitochondrial CMPK2 mediates immunomodulatory and antiviral activities through IFN-dependent and IFN-independent pathways Mitochondrial protein CMPK2 regulates IFN alpha-enhanced foam cell formation, potentially contributing to premature atherosclerosis in SLE Host metabolic reprogramming in response to SARS-CoV-2 infection: A systems biology approach Metformin inhibition of mitochondrial ATP and DNA synthesis abrogates NLRP3 inflammasome activation and pulmonary inflammation Decitabine Reactivation of FoxM1-Dependent Endothelial Regeneration and Vascular Repair for Potential Treatment of Elderly ARDS and COVID-19 Patients High expression of RRM2 as an independent predictive factor of poor prognosis in patients with lung adenocarcinoma Bioinformatics and machine learning approach identifies potential drug targets and pathways in COVID-19 Poly (I:C) transfection induces mitochondrial-mediated apoptosis in cervical cancer Exploiting poly(I:C) to induce cancer cell apoptosis Poly(I:C)-Mediated Death of Human Prostate Cancer Cell Lines Is Induced by Interleukin-27 Treatment Toll-like receptor 3 agonist poly I:C reinforces the potency of cytotoxic chemotherapy via the TLR3-UNC93B1-IFN-β signaling axis in paclitaxel-resistant colon cancer poly(I:C) costimulation induces a stronger antiviral chemokine and granzyme B release in human CD4 T cells than CD28 costimulation Intranasal treatment with poly(I·C) protects aged mice from lethal respiratory virus infections Induction of Innate Immune Response by TLR3 Agonist Protects Mice against SARS-CoV-2 Infection Interdependent and independent roles of type I interferons and IL-6 in innate immune, neuroinflammatory and sickness behaviour responses to systemic poly I:C Comprehensive evaluation of poly(I:C) induced inflammatory response in an airway epithelial model IL-6 expression promoted by Poly(I:C) in cervical cancer cells regulates cytokine expression and recruitment of macrophages Mitomycin C: Mechanism of action, usefulness and limitations Mitomycin C: A clinical update Mitomycin C plus S-1 as second-line therapy in patients with advanced gastric cancer: A noncomparative phase II study Mitomycin C: New strategies to improve efficacy of a well-known therapy Recent Results Cancer Res Role of gemcitabine in cancer therapy lycorine and oxysophoridine inhibit nov.vel coronavirus (SARS-CoV-2) in cell culture Influence of 4 -Substitution on the Activity of Gemcitabine and Its ProTide Against VZV and SARS-CoV-2 Comparison of Antiviral Activity of Gemcitabine with 2'-Fluoro-2'-Deoxycytidine and Combination Therapy with Remdesivir against SARS-CoV-2 Association of Antineoplastic Therapy With Decreased SARS-CoV-2 Infection Rates in Patients With Cancer Repurposing of Hydroxyurea against COVID-19: A Promising Immunomodulatory Role The Use of Hydroxyurea in the Treatment of COVID-19 The role of tamoxifen in the treatment and prevention of breast cancer Estrogen and Androgen Receptor Inhibitors: Unexpected Allies in the Fight Against COVID-19 Tamoxifen and clomiphene inhibit SARS-CoV-2 infection by suppressing viral entry Estrogen Receptor Modulators in Viral Infections Such as SARS-CoV-2: Therapeutic Consequences Anti-inflammatory properties of curcumin, a major constituent of Curcuma longa: A review of preclinical and clinical research Safety and anti-inflammatory activity of curcumin: A component of tumeric (Curcuma longa) Antioxidant and anti-inflammatory properties of curcumin The clinical applications of curcumin: Current state and the future Potential effects of curcumin in the treatment of COVID-19 infection Antiviral and immunomodulatory activity of curcumin: A case for prophylactic therapy for COVID-19 Curcumin as a Potential Treatment for COVID-19 Oral Curcumin With Piperine as Adjuvant Therapy for the Treatment of COVID-19: A Randomized Clinical Trial Effectiveness of Curcumin on Outcomes of Hospitalized COVID-19 Patients: A Systematic Review of Clinical Trials We thank the National e-Science Infrastructure Consortium (http://www.escience.in.th accessed on 7 April 2022) for providing computational resources that have contributed to the complete results in this article. We would like to thank William Lao, the University of New South Wales, Australia, for his constructive reviews. We also thank our families for supporting our efforts until this article is published. The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.