Rapid Identification of Druggable Targets and the Power of the PHENotype SIMulator for Effective Drug Repurposing in COVID-19

The current, rapidly diversifying pandemic has accelerated the need for efficient and effective identification of potential drug candidates for COVID-19. Knowledge on host-immune response to SARS-CoV-2 infection, however, remains limited with very few drugs approved to date. Viable strategies and tools are rapidly arising to address this, especially with repurposing of existing drugs offering significant promise. Here we introduce a systems biology tool, the PHENotype SIMulator, which – by leveraging available transcriptomic and proteomic databases – allows modeling of SARS-CoV-2 infection in host cells in silico to i) determine with high sensitivity and specificity (both>96%) the viral effects on cellular host-immune response, resulting in a specific cellular SARS-CoV-2 signature and ii) utilize this specific signature to narrow down promising repurposable therapeutic strategies. Powered by this tool, coupled with domain expertise, we have identified several potential COVID-19 drugs including methylprednisolone and metformin, and further discern key cellular SARS-CoV-2-affected pathways as potential new drugable targets in COVID-19 pathogenesis.


INTRODUCTION
The rapid emergence and spread of the virulent novel severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) has hijacked and largely disrupted human civilization as we know it, bringing about countless global challenges but also the urgent need for innovative vaccine and drug discovery 1 . Soon after its emergence in Wuhan (Eastern China) in late 2019, coronavirus disease 2019 (COVID- 19) was declared a pandemic by the World Health Organization 2 , as it continuously spreads and holds the world hostage. As of February 2021, over 100 million COVID-19 cases and 2 million deaths have been reported worldwide. It is apparent that our civilized, well-organized and hitherto functioning societies were not adequately prepared nor equipped to deal with the high infectivity, transmissibility, mortality and global impact of the COVID-19 pandemic. While vaccine development and deployment is well underway, widespread distribution remains challenging, and at present only an antiviral (Remdesivir) and glucocorticoids (Dexamethasone/ Methylprednisolone) have been approved for treatment of severe COVID-19. Recently, the virus-neutralizing antibody cocktail (Casirivimab and Imdevimab, termed REGN-COV2) also received emergency use authorization for treatment of mild to moderate COVID-19 in high-risk patients 3 . Otherwise, no established drug is available to prevent or adequately treat COVID-19 and in the absence of a clear etiological understanding, treatment has remained largely supportive and symptomatic 4,5 .
Making matters worse, the mutating virus is now posing additional challenges.
Therefore, next to in vitro studies, in silico studies are of great value for rapid and effective drug discovery. Indeed, computational structure-based drug design and immunoinformatics have recently resulted in identification of potential SARS-CoV-2 target proteins and drugs that are being selected for further testing 4,6,7 . Another promising avenue for obtaining effective and readily available therapeutic strategies is the repurposing of drugs already approved for other indications. Drug repurposing strategies provide an attractive and effective regulatory elements 18 . The algorithm thus computes, under-specified biological contexts, by iteratively propagating the effects and alterations of one or more biomolecules (differentially expressed genes (DEGs), proteins, microRNAs, or metabolites), thus making use of published virus-human interaction data 7 . Here we compare our results with available data from recently published in vitro studies based on transcriptomics and proteomics in different model systems 5,10 . Relevant and significantly affected pathways are further detailed on a protein interaction level. Finally, we show the potential of the PHENSIM in selecting promising hypothesis-driven COVID-19 drug candidates, which has applicability to other diseases and broader aspects of clinical practice, thereby outlining the potential power of PHENSIM in drug repurposing in COVID-19 and beyond.

PHENSIM model: from in vitro to in silico
Innovative approaches to rapidly elucidate a pathogens' mechanism of action have proven crucial for containing the global burden of communicable diseases. The PHENSIM approach, described here, is based on the definition of a newly introduced protocol for in silico simulation of novel emerging pathogens, such as SARS-CoV-2, and it aims at elucidating distinct hostresponses and molecular mechanisms triggered by that particular pathogen, all while defining possible candidate drugs for indication repositioning.
For our strategy to be viable, even when only limited direct knowledge is available on the host-pathogen interaction, we need direct infection (in vitro) data that can be exploited to predict such interactions. To acquire this knowledge, we therefore employ transcriptomic and proteomic experiments of in vitro infected vs. normal, pathogen-free cell lines. When available, we leverage Differentially Expressed Genes (DEGs) as a means to simulate the direct and indirect effect of the virus on a host without a priori knowledge regarding the mechanism of infection.
Using DEGs as input for our cell PHENotype SIMulator PHENSIM [ 15,19 we define a signature of pathogen predicted effects on human pathways (pathogen alterations profile; here termed the "viral signature"; see Fig. 1). To build the viral signature, we use pathway endpoints; an endpoint is a biological element in a pathway whose alteration, based on current knowledge, affects the phenotype in a specific way 20 .
By leveraging PHENSIM we aimed to determine the impact of such viral infection induced alterations on an array of human cell lines in silico. Simulation results are used to define a "viral signature", that can then be employed to identify candidate drugs. Once a cellspecific SARS-CoV-2 viral signature is defined, potential repositioning drugs can be identified by building a "drug signature" database queried by means of a similarity measure using pathway endpoints (Fig. 1). Given a candidate drug identified through a database (i.e. Drugbank or Pubchem) and literature (Pubmed) search, we define all known targets and alterations (up/down-regulations caused by the drug). Alterations are then provided as input to PHENSIM, together with the corresponding cell-specific viral signature. Next, distinct endpoint pathways 20 are identified and resulting drug signatures relating to a specific candidate can subsequently be compared with acquired viral signatures to evaluate the inhibitory potential of that candidate drug. Both viral and drug signatures are collected in a database, where a similarity search is performed using a Pearson correlation ρ(x,y) since the propagation algorithm is linear in time complexity 21 ; see methods section equation (1). All drugs whose correlation with the virus is negative (green) are considered possible repositioning candidates, since they predict inhibition of the viral signature, whereas a positive correlation (red) suggests exacerbation of the viral signature when introducing the candidate drug.

Validation of PHENSIM transcriptomic strategy in SARS-CoV-2-infected host cells
To validate our PHENSIM model on a transcriptomic level in the context of SARS-CoV2, we sought to replicate the in vitro experiments using publicly available data presented by Blanco-Melo et al 10 . The in-depth transcriptomic analysis of SARS-CoV-2 elicited host-response by Blanco-Melo et al. recently revealed an inappropriate inflammatory response driven by reduced innate antiviral defenses, with low or delayed type I and type III interferon (IFN) and exaggerated inflammatory cytokine response, with elevated chemokines and IL-6 10 .
As SARS-CoV-2 largely affects the lungs and respiratory tract, and because of its apparent affinity for lung tissue, the authors make use of several respiratory epithelial cell lines to assess the transcriptomic host-response. Here we use PHENSIM to reproduce transcriptomic effects in silico, as described in vitro for the following cell lines, namely undifferentiated normal human bronchial epithelial (NHBE) cells, cultured human airway epithelial cells (Calu-3) cells and A549 lung alveolar cells. The comparison of these results is depicted in Fig. 2. A549 cells are described to be relatively non-permissive to SARS-CoV-2 replication in comparison to Calu-3 cells, which is attributed to low expression of the viral entry receptor angiotensin-converting enzyme (ACE)2 10,22 . Thus, A549 cells were transduced with human ACE2 (A549-ACE2), which enabled apparent SARS-CoV-2 replication at low-MOI (multiplicity of infection of 0.2).
Furthermore, to induce significant IFN-I and -III expression, a high MOI of approximately 2-5 was necessary.
Here we leveraged the data published by Blanco-Melo et al. to run our PHENSIM simulation pipeline. In Fig. 2A  at least for the top DEGs (Fig. 2B). To quantify the overall predictive accuracy of PHENSIM, genome-wide transcriptomic data was assessed for all scenarios as described in Fig. 2 specificity for this in-depth SARS-CoV-2 transcriptomic analysis. Furthermore, the positive predictive value (PPV) and False negative rate (FNR) are shown for each tested scenario (see Table 1).
In order to further verify PHENSIM's robustness in whole genome pathway analysis, we next explored PHENSIM's ability to predict significantly affected signaling pathways in SARS-CoV-2 infection. In Fig. 2C we highlight PHENSIM's predicted perturbation of a select set of affected pathways during infection, as recently identified to be of importance by Catanzaro et al. For further verification of our PHENSIM pathway analysis prediction in silico, we compared our results with those obtained using our previously described MITHrIL (Mirna enrIched paTHway Impact anaLysis) tool 20 to analyze the Blanco-Melo et al. acquired in vitro data ( Fig. 2D-E). Given DEGs, MITHrIL first computes a perturbation for each gene in the metapathway (as described in Methods section). The perturbation can be considered as the predicted state that the node will have given the input DEGs. Next, we sum the perturbation of all nodes for each pathway to acquire the "accumulated perturbation," or the Accumulator. The accumulator is equivalent to a pathway expression and is a sum of all perturbations computed for that particular pathway. MITHrIL pathway analysis for A549-ACE2 at low viral load (MOI 0.2) revealed Chemokine, JAK-STAT, PI3K-Akt signaling and cytokine-cytokine interaction as a few of the top upregulated pathways, according to impact (circle size), significance (color-gradient for adjusted p-value) and accumulated perturbation computed for that particular pathway (accumulator).
For A549-ACE2 at high viral load (MOI 2.0; Fig. 2E), next to similar pathways at low viral MOI, Toll-like receptor (TLR) and NOD-like receptor signaling were among the top pathways observed, corresponding to the observation that high viral MOI was needed to induce significant type I IFN signaling 10 . Interestingly, both at low and high MOI various metabolic pathways were significantly affected with a negative accumulator. Overall, the MITHrIL analysis results show the most affected pathways to be similar to the PHENSIM in silico predicted results.

Modeling proteomics in SARS-CoV-2-infected host cells leveraging PHENSIM
Using combinatorial profiling of proteomics and translatomics to study host-infection on a cellular and molecular level give opportunity to study relevant viral pathogenicity in the search of potential drug targets 5

PHENSIM proteomic validation
To validate PHENSIM on a proteomic level, we used our in silico approach to replicate the in vitro SARS-CoV-2 infection of human Caco-2 cells 5 Table 2).

PHENSIM proteomics from in vitro to in silico
Next, to compare the Reactome-based in vitro functional pathway analysis to our PHENSIM in silico approach, a representative selection of significantly affected pathwayscorrectly predicted by PHENSIM -is depicted in Fig. 3D. Pathways were selected according to the cellular mechanisms highlighted by Bojkova et al. 5 . The centered heatmap shows an increasing activity score (top to bottom) as predicted by PHENSIM for each pathway. An indepth analysis of proteomic pathway at 24hrs revealed distinct upregulation of various pathways involving cellular metabolism such as fatty acid degradation, glycolysis and glyconeogenesis, carbon metabolism, inflammatory and immune signaling pathways and also cellular senescence signaling pathways (Fig. 3D). A selection of KEGG pathways similar to the Reactome protein interaction pathways are further highlighted in more detail, depicting the main protein-protein interactions for that particular pathway, with upregulated proteins depicted in red and downregulated proteins in blue. As authors used Reactome as their pathway knowledge-base, here we show only KEGG pathways matching their Reactome counterpart.

PHENSIM predicts a metabolic signature in SARS-CoV-2 infection in silico
As a metabolic signature was identified by PHENSIM's proteomic in silico simulation of SARS-CoV-2 infection (Fig. 3D), we next assessed the degree of intersection between the perturbed genes of these metabolic pathways in order to reject the hypothesis that a common set of altered proteins is driving the significant perturbation of these closely related metabolic pathways. All metabolic pathways considered essential for SARS-CoV-2 infection according to the acquired Bojkova et al. data (Fig. 3) were included in the analysis (FDR-adjusted p-value < 0.05) and a PHENSIM activity score was determined (see Supplementary Table S2

PHENSIM Drug repurposing strategy for COVID-19
The next step in our PHENSIM approach is the employment of our drug strategy in order to test candidate drugs for potential COVID-19 repurposing. This approach takes advantage of existing knowledge on drug-related pharmacology and toxicology for rapid therapeutic selection 8 . As schematically described in Fig. 1, once a cell-specific viral signature is defined, it can be exploited to search for possible repositioning candidates by leveraging our select drug signature database. We used a Pearson correlation p(x,y) to compare the viral and drug signatures, which gives rise to a correlation score specific to that candidate drug, computed for SARS-CoV-2 infection in a particular setting. Here we set out to test a selection of hypothesis-and datadriven candidate drugs as shown in Fig. 4. One such drug which regrettably failed to live up to its anticipated potential to effectively treat COVID-19 is the antimalarial drug hydroxychloroquine (HCQ), currently approved for rheumatologic implications, although associated with cardiac toxicity [26][27][28] .
Although the efficacy of corticosteroids in viral acute respiratory distress syndrome

PHENSIM Methylprednisolone treatment of SARS-CoV-2 infected host cells in silico
As a next step in our drug repurposing efforts, we simulate the simultaneous host-cell infection of SARS-CoV-2 and in silico treatment with Methylprednisolone (MP), hereby combining the drug action and pathogen infection on a host-cell, in order to further assess MP as top candidate. We simulate SARS-CoV-2 viral infection and simultaneous MP treatment in silico, in order to more closely resemble the in vivo situation (Fig. 5). The heatmap in Fig

Discussion
The current pandemic has accelerated the need for efficient and effective identification of potential drug candidates for COVD-19 pathogenesis. Knowledge on host-immune response to SARS-CoV-2 infection, however, remains limited with very few drugs approved to date. Various viable strategies and tools are rapidly arising to address this, where repurposing of existing drugs can offer a feasible mechanism of deployment (4). Here we introduce one such strategic approach, the PHENotype SIMulator, which allows the modeling of SARS-CoV-2 infection in silico and implementation to select promising candidates for further in vitro and in vivo analysis.
We show that PHENSIM can effectively be used to 1) determine the viral effects on cellular host-immune response and cellular pathways and 2) evaluate a myriad of therapeutic strategies.
As previously described, PHENSIM uses a probabilistic randomized algorithm to compute the effect of a particular biological scenario on gene regulation, protein expression, miRNA and metabolite involvement, with use of KEGG meta-pathway analysis 15 . Here, we simulate SARS-CoV-2 viral infection based on publicly available data to acquire a specific cellular SARS-CoV-2 signature. To verify a transcriptomic-based PHENSIM strategy, we compared the in silico simulation results to publicly available transcriptomic data of SARS-CoV2 infected cell lines 10 . We find that PHENSIM performs at a high overall accuracy with high PPV, sensitivity and specificity for all airway and lung-related cell lines evaluated. Key SARS-CoV2 infection-related signaling pathways could be discerned as such, comprising the viral signature.
PHENSIM predictive performance was further validated using our previously described MITHRIL transcriptomic pathway analysis 20 , showing similar results. Interestingly, key signaling pathways proposed to be crucial in SARS-CoV2 infection 4 , were shown to be significantly perturbed in all cell lines studied in silico using PHENSIM, simulation offering promising potential molecular drug targets.
The PHENSIM strategy was also suitable for a proteomics/translatome-data based approach. PHENSIM simulation was compared to published SARS-CoV2 infection specific proteomic effects in host cell lines 5 . Comparing simulation results to proteomic data 24 hours post infection showed, for the proteins available in KEGG, high accuracy with a PPV, sensitivity and specificity well above 97%. Inhibition of several of these protein associated pathways pathways was shown to prevent viral replication in human cells 5 .
We next used the transcriptomic-based PHENSIM approach to compare the viral signatures, computed with respect to model cell lines, to in-silico-derived drug signatures for a selection of drugs analyzed. Our overall correlation results show several potential drug repurposing candidates negatively correlating with SARS-CoV-2, varying from corticosteroids such as MP (already approved for treatment of COVID-19 patients) to biologics such as BTKinhibitors that are currently being studied in clinical trials 31 to metformin 32 . Individual signaling pathway contribution to the observed correlation score could be further delineated for each individual drug, providing specific targets for in depth analysis and potential for pathway-specific therapeutic targeting. As expected, the individual pathways most targeted by the in silico drug interventions (Fig. 4C&D) were similar to pathways found most perturbed by PHENSIM during host transcriptomic response to SARS-CoV-2 viral infection (Fig. 2C&D), emphasizing their potential therapeutic effects. HCQ, although hypothesized to be a good potential candidate to treat COVID-19, has not proven effective in vivo 33 . The exact reason why HCQ has failed in COVID-19 remains to be fully understood. Interestingly, COVID-19 is associated with a variety of hematologic complications 34 , and increased HCQ use during the COVID-19 pandemic has induced the emergence of methemoglobinemia, including tissue hypoxia and reduced oxygenation 35,36 amongst others. Evidently, evaluating the risk-benefit ratios -drug safety and efficacy -is crucial when selecting drugs to be repurposed for COVID-19 8 , which particularly holds true for HCQ 27,28,37 .
In Fig. 4, we depict our drug repurposing PHENSIM approach that functions as a screening tool for initial drug candidate screening, based on the anti-correlation of the viral and drug signatures, and gives a particular score for each candidate; a negative correlation constitutes higher potential for that particular drug. This broader correlation approach described in Fig. 4 and Supplementary Fig. S5 and S6 can be used to screen large sets of candidate drugs. Next, as depicted in Fig. 5, a more dynamic and extensive analysis can be simulated, in order to simulate the interaction between the SARS-CoV-2 host-cell infection (column B) and subsequent in silico treatment with a candidate drug, here Methylprednisolone (column C).
Although Methylprednisolone is a known broad-spectrum corticosteroid, with clear anti-  Supplemental Table S1).
As discussed, the in silico model presented here provides an interesting framework that could be further developed and expanded, achieving a more complete cell signature with input of available data on processes such as cell-cell communication through ligand-receptor complexes 50 or viral immune evasion e.g. 51 Fig. 4B). As shown, the sequence of top candidate drugs for repurposing is slightly different depending on the cell-type, viral load (MOI), and expression of the viral entry receptor ACE2 (see Fig. 4 and Supplementary Fig. S2). This speaks to the variability of how the virus might affect specific cell types and tissues, even within the same organ system such as the bronchial (NHBE), airway epithelial (Calu-3) and lung alveolar (A549) cells e.g., pointing to the difficulty in specifically targeting this viral infection therapeutically. It will also depend on the stage of infection and disease state, as to which course of treatment or combination of treatments will be optimal.
As demonstrated by our results, we believe that the PHENSIM system provides a multitude of powerful systems biology functions and implements them easily and efficiently.
PHENSIM is a simulation algorithm which follows the biological processes modeled by pathways. Therefore, PHENSIM is able to make a prediction of such processes and not only of the final effect, going beyond methods based on pathway enrichment. Furthermore, since pharmacological treatments may depend on the state of biological processes, PHENSIM may be of more appropriate use in this context. Comparison with other simulation algorithms such as BIONSI 13,14 has shown excellent performance by PHENSIM 15 . PHENSIM creates and builds on interpretable and intervenable mechanistic bio-chemical models, rather than combinatorial and statistical "black-box" models for joint stationary distribution of biological data, as in, say proteinprotein interaction (PPI) networks, Graphical or Deep-net models. potentially chronic disease such as inflamm-aging and diabetes. In conclusion, our PHENSIM approach will enable more rapidly initiated clinical trials and accelerated regulatory review of already pre-selected drugs with a high repurposable potential.

PHENotype SIMulator (PHENSIM) approach
PHENSIM is a systems biology approach, simulating the effect of the alteration of one or more biomolecules (genes, proteins, microRNAs, or metabolites) in a specific cellular context using KEGG (Kyoto Encyclopedia of Genes and Genomes) meta-pathways 15,20 . The meta-pathway concept, introduced by us previously 61  To start the simulation, PHENSIM requires a set of biomolecules as input, their direction of deregulation (activation/up-regulation or inhibition/down-regulation), and a set of inactive biomolecules in the cellular context (cell lines, tissue, e.g). The algorithm uses these details to compute synthetic Log Fold Changes (LogFC). Synthetic LogFCs are computed by sampling the normal distribution fitted to the actual LogFCs of a particular gene, as described previously 15 . Such values are then propagated within the meta-pathway, using the MITHrIL (Mirna enrIched paTHway Impact anaLysis) pathway perturbation analysis 20 . MITHrIL determines how local change can affect the cellular environment by computing a "perturbation". For each gene in the meta-pathway, the perturbation reflects its expected change of expression/activity (negative/positive for down-/up-regulation, respectively). Finally, these results are collected and synthesized using two values: the "Average Perturbation" and the "Activity Score" (AS). Given a node, the average perturbation is the mean for its perturbation values computed at each simulation step. It reproduces the expected change of expression for the entire process. The function of AS is twofold: 1) the sign gives the type of predicted effect (activation(+); inhibition(-)), 2) the value is the log-likelihood that this effect will occur. Together with AS, PHENSIM also calculates a p-value which determines how biologically relevant the predicted alteration is for the phenomena being simulated.
All p-values computed by PHENSIM are corrected for multiple hypotheses using the qvalue algorithm 65 . To determine this probability, PHENSIM randomly selects genes in the metapathway and runs the simulation on this random set. By repeating this procedure (n=1000 for our simulations), it is possible to empirically estimate the probability that a node has a higher activity score than the observed one. For this reason, we can employ such a value to determine which alterations are most specific for a particular infection, gaining novel hypotheses on the molecular action of the pathogen.

PHENSIM pathogen alterations profile
Our approach defines a protocol for the in silico simulation of emerging pathogen infection, aimed at defining candidate drugs for repositioning. First, we find a representation of the pathogen in the KEGG meta-pathway, which allows us to perform simulations. For a novel pathogen, such as SARS-CoV-2, interactions with the host genes might be unknown. Therefore, To build the pathogen signature, we use pathway endpoints; An endpoint is a biological element in a pathway whose alteration, based on current knowledge, affects the phenotype in a specific way 20 . Given the output of this simulation, we collect all Endpoint Activity Scores in a single signature, the 'pathogen alterations profile'. This profile can be exploited to search for possible repositioning candidates, by building a drug signature database queried by means of a similarity measure (Fig. 1). When building the simulation profile, we do not use any p-value.
Indeed, we need to consider not only the alterations, which are the most specific for a particular infection, but also the alterations caused by any cellular response to the infection. Since the pvalue represents the biological relevance for the phenomena that is being simulated, we can ignore its value to build the signature.

PHENSIM drug signature database
Given a particular drug identified through databases (i.e. Drugbank or Pubchem) and literature (Pubmed) searches, we define all known targets and their alterations (up/down-regulations caused by the drug), and these alterations are provided as input to PHENSIM together with the same cellular context specified for the viral simulation. The results are used to define a drug signature using pathway endpoints as described above, which are collected in a database used for repositioning.

PHENSIM drug repurposing approach
Our drug repurposing methodology is based on a similarity search performed on the drug signature database. Given a pathogen profile computed with PHENSIM, we use a correlation function to scan through each record in a drug profile database. This procedure yields a ranking on each drug in the range [-1;1], where negative values indicate that the virus alteration profile is opposite to the drug and positive values indicate the reverse; drugs with a negative correlation are considered possible candidates for repositioning. In our experiments, we employ a Pearson correlation function to run the similarity search. Since PHENSIM is based on MITHrIL pathway perturbation analysis, which computes results in a log-linear space 20 , we can assume a Pearson correlation is sufficient to determine similarity between the viral and drug signature.
A key characteristic of this approach is the capability to simulate both single and drug combinations. Furthermore, PHENSIM also provides a framework for extending pathways by adding new nodes and edges coming from results in the literature as well as other reputable sources.
Finally, to assess whether each drug candidate targets relevant infection processes, we decomposed the Pearson correlation in terms of KEGG pathways and reviewed the results.
More in detail, let D and V be drug and pathogen alteration profiles, respectively. That is, D[e] is the activity score of endpoint "e" computed by PHENSIM for a drug simulation, and V[e] is the activity score for the same endpoint in the pathogen simulation. Pearson correlation ( , ) can be written as equation (1) A significant feature of this partial correlation approach is that we obtain the total correlation by summing up all values for each pathway P. Therefore, we can determine which biological processes are impacted by the drug administration.

PHENSIM combined drug/pathogen simulation
To further evaluate whether the results of the correlation could be confirmed by PHENSIM, we devised a strategy to simultaneously simulate drug action and pathogen infection on a host cell line. First, we collected DEGs between pre-and post-infection samples as described in the previous section. Then, given a drug, we gather its known targets and their alterations (up/down-regulations caused by the drug) through databases (i.e. Drugbank or Pubchem) and literature (Pubmed) searches. Therefore, we extend the meta-pathway by adding two nodes, representing the virus and the drug, respectively. The virus node is connected to each DEG with an edge as described in the "PHENSIM pathogen alterations profile" section. Then, an activating (inhibiting) edge is added between the drug node and a target, for any up-regulated (/down-regulated) target. Finally, we can run a simulation by giving as input the simultaneous upregulation of both virus and drug nodes (results depicted in Fig. 5).

PHENSIM method validation
To determine the efficacy of our model we used several datasets obtained in the context of SARS-CoV-2 infection. For each dataset, we computed the genome-wide Log Fold Changes (FC). As PHENSIM does not require any quantitative information, DEGs were termed upregulated if LogFC>0.6, and downregulated if LogFC<-0.6.

PHENSIM CD147 gene and pathway extension
Prior to running all simulations, we verified if viral entry points were present in KEGG to better represent viral activity. SARS-CoV-2 is said to invade host cells via these two receptors: angiotensin-converting enzyme 2 (ACE2) and CD147 (also known as Basigin or EMMPRIN) 49 .
In KEGG the latter gene is missing and to extend our simulation model, a new node, representing CD147 was added and connected with its known interactions and downstream nodes according to literature (Supplementary Table S1). The incoming edges to CD147 represent the possible activators/inhibitors (upstream genes), the outgoing edges represent the actions performed by CD147 towards its downstream genes. In Supplementary Table S1, we report a list of up-and downstream genes added to our simulations and their respective references containing interaction details. CD147 is a transmembrane protein of the immunoglobulin super family, expressed in many tissues and cells, acting as the main upstream stimulator of matrix metalloproteinases (MMPs) and playing a crucial role in intercellular recognition 66 . Over the last decade, several groups have shown that CD147 acts as a key molecule in the pathogenesis of several human diseases including infectious diseases (HIV, HBV, HCV, KSHV) 66 , and it has now been posed to recognize and internalize/endocytose SARS-CoV-2 in certain cell types 49 .

PHENSIM transcriptomic reliability assessment
To assess the reliability of the results, we focused on two fronts: ii) the ability of PHENSIM to predict genes altered in the expression data, and ii) the ability to predict the correct direction of the alteration. In detail, we define as altered all genes having an absolute LogFC>0. 6 Following the same procedure used by Blanco-Melo et al., raw counts were normalized and analyzed for differential expression using the DESeq2 pipeline. All genes with an FDR-adjusted p-value<0.05 and absolute LogFC>0.6 were considered differentially expressed. Non-expressed genes for a specific cell-line were defined as genes that showed an average read count lower than 10.

PHENSIM proteomic approach
To determine if our methodology can also exploit proteomic data, we leveraged data from  (Fig. 3D). Finally, since we found many similar metabolic pathways significantly affected in silico as described in vitro by Bojkova et al., we aimed to determine if a core set of proteins was common between pathways; results are displayed in the VENN diagrams in Supplementary Fig. S4 and Table S2.

Performance Evaluation: PHENSIM Genome-wide and Proteome network analysis
We

Pathway Analysis
Pathway analysis was applied to the transcriptomics data to determine which biological processes were altered by the viral infection. We used 4 pathway analyses approaches to assess the most impacted pathways:  Fig. S3 and Supplementary material).

Statistical analysis
Statistical methods for transcriptomics and proteomics were applied as described by Blanco-Mello and Bojkova et al. respectively 5,10 . For transcriptomic data, raw counts were normalized and analyzed for differential expression using the DESeq2 pipeline as previously described 72 .
All genes with an FDR-adjusted p-value<0.05 and absolute LogFC>0.6 were considered differentially expressed. In addition, we considered all genes showing an average read count<10 as non-expressed. All p-values computed by PHENSIM are corrected for multiple hypotheses testing, using the q-value algorithm 65 . For proteomic data, Normalized LC-MS/MS data were downloaded and significance was tested using unpaired two-sided Student's t-tests with equal variance assumed. All values were aggregated on a pathway basis to compute an Accumulator and a p-value, and p-values were adjusted for multiple hypotheses using the Benjamini-Hochberg FDR correction. Results were filtered by an FDR-adjusted p-value of 0.05 and ranked using the Accumulator.