To predict multi-target drug candidates for COVID-19, we developed a drug scoring method S_PCN that was calculated as follows:

where n is the target number of a drug community; S_DisGi is the therapeutic target score of target protein i and defined in Equation (1); S_DN is CoNet score of a drug community and defined in Equation (3). Finally, we computed the S_PCN values of 2,303 FDA-approved drugs.

To infer disease genes, we developed a scoring function to calculate the disease gene score (S_DisG) of DEGs. The S_DisGi score of a DEG (i) consists of the gene expression significance p value (S_PVi) between disease and corresponding normal samples, its druggability (S_{Druggabilityi}) measuring the number of interacting drugs, and its involvement significance (S_HiSBiNi) in pathways and subsystems using HiSBiN (Figure a). S_DisGi is defined as follows:

where i is the DEG (i), and w_PV is the weight of S_PV, w_D is the weight of S_Druggability, and w_p is the weight of S_HiSBiN. Here, the weights were set to 0.5, 2, and 1 to w_PV, w_D, and w_p respectively. We statistically analyzed and normalized these three scores ranging from 0 to 1 across all the DEGs and 0 ≤ S_DisG ≤ 3.5. The detailed scoring method and scheme for identifying therapeutic target score for the gene i are as follows, and the example of PADI4 is shown in Figure.

In the first stage, we consider the gene importance for a disease by using S_PVi, which is the maximum p value of target i of gene expression in two SARS-CoV GEO datasets (GSE5972: precrisis, crisis and GSE1739) and was computed as follows:

The p value was calculated by the GEO2R analysis tool, s was number of datasets.

In the second stage, S_{Druggabilityi} is used to determine the draggability of a DEG i, which is considered as the number of targeting drugs derived from the Homopharma. For example, for the gene PADI4, which is one of the DEG of methotrexate community, has reported as being the target of four drugs (e.g., Azithromycin, Tetracycline, Streptomycin, and Citrulline) in the recorded drug-target pair database (i.e., Drugbank) and 42 predicted approved drugs derived from Homopharma (Figure b).

In the third stage, to evaluate how DEGs and drug targets are involved in pathways during SARS-CoV-2 infection, we first collected 342 human pathways containing 7,806 proteins (genes) from the Kyoto Encyclopedia of Genes and Genomes (KEGG) database³. We then modified the hierarchical systems biology model (HiSBiM)² from our previous work to measure the involvement of drug targets in the KEGG pathways, including DEG identification, pathway enrichment, and subsystem contribution by ‘hierarchical systems biology network (HiSBiN)’. The DEGs for COVID-19 were identified based on fold change and p value between normal and SARS-CoV (or SARS-CoV-2) infection samples. To yield KEGG pathway enrichment of the DEGs, the p value of hypergeometric distribution was calculated as follows:

where N is the number of all genes in the KEGG database, M is the number of DEGs, n is the number of genes and k is the number of DEGs in a specific KEGG pathway.

To understand how each DEG is involved in KEGG subsystems and pathways, we defined the S_HiSBiNi of a DEG i as follows:

where  and Z_Pathi are the meta-z-scores of the subsystems and pathways in which a DEG i is involved; Z_t and Z_r are the z-scores of subsystem t and pathway r, respectively; and m and n are numbers of a DEG-involved subsystems and pathways of 184 COVID-19 KEGG pathways, respectively. The z-score of each pathway was the p value transformation using a standard normal distribution. The detailed example of gene PADI4 was shown in Figure c.

For construction of the PharmaCoNets of each drug with its target proteins (drug community) in disease (e.g. coronavirus infectious disease), we evaluated the enrichment of co-expressed gene pairs between all genes of drug communities and disease-related pathways using the gene expression profiles of infected by SARS-CoV/SARS-CoV-2 human blood/ nasopharyngeal tissue samples. Based on our previous work, we have systematically study cancer membrane protein-regulated networks (CaMPNet) by the construction of membrane protein community and its regulated cancer-related pathways through this methodology⁴. We then modified from our previous work to measure the involvement of drug communities regulation of diseases-related in the KEGG pathways, including drug community identification, and pathway enrichment contribution by ‘co-expression network (CoNet)’.

For each coronavirus disease type, two DEG’s with Pearson’s correlation coefficient |Pearson’s r| ≥ h (here, h is set to 0.7 based on a large effect size) were considered a co-expressed pair. For each DEG as an involved gene in the drug community, we first used the co-expressed pairs between it and a disease-related pathway to determine its involvement (enrichment p value is converted to z-score) for this pathway in each cancer type based on hypergeometric distribution. Moreover, for each drug community, we measured the involvement between its involved genes and all the DEGs of regulated pathways in disease. To yield KEGG pathway enrichment of the drug communities, the p value of hypergeometric distribution was calculated as Equation (2), where i and n are the numbers of co-expressed gene pairs and all the combinational gene pairs, respectively; x is the observed co-expressed gene pairs with |Pearson’s| ≥ 0.7; M and N are the total numbers of all the co-expressed gene pairs and combinational gene pairs, respectively, between all the involved DEGs of the drug community and all the DEGs in 342 KEGG pathways.

To understand how each drug community is involved in KEGG pathways, we defined the CoNet score (S_DN) of a drug community as follows:

where S_DPj is the drug community-regulated pathway score (z-score converted by enrichment p value) of the pathway j in which drug community are involved; and m is numbers of disease DEGs enriched pathways (i.e., here is the 184 COVID-19 KEGG pathways identified by SARS-COV DEGs). The z-score of each pathway was the p value transformation using a standard normal distribution.