A Functional Network Model of the Metastasis Suppressor PEBP1/RKIP and Its Regulators in Breast Cancer Cells

Simple Summary We extracted known interactions between metastasis suppressors and their regulators from scientific studies. Next, we used publically available data of drug treatments to identify which of them potentially perturbs these interactions. Finally, we studied the effect of several of these drugs on a particular metastasis suppressor called RKIP and found that our model accurately predicted its regulation in breast cancer cells. Our approach can discover alternative mechanisms of existing cancer drugs and repurpose them in different disease types. Abstract Drug screening strategies focus on quantifying the phenotypic effects of different compounds on biological systems. High-throughput technologies have the potential to understand further the mechanisms by which these drugs produce the desired outcome. Reverse causal reasoning integrates existing biological knowledge and measurements of gene and protein abundances to infer their function. This approach can be employed to appraise the existing biological knowledge and data to prioritize targets for cancer therapies. We applied text mining and a manual literature search to extract known interactions between several metastasis suppressors and their regulators. We then identified the relevant interactions in the breast cancer cell line MCF7 using a knockdown dataset. We finally adopted a reverse causal reasoning approach to evaluate and prioritize pathways that are most consistent and responsive to drugs that inhibit cell growth. We evaluated this model in terms of agreement with the observations under treatment of several drugs that produced growth inhibition of cancer cell lines. In particular, we suggested that the metastasis suppressor PEBP1/RKIP is on the receiving end of two significant regulatory mechanisms. One involves RELA (transcription factor p65) and SNAI1, which were previously reported to inhibit PEBP1. The other involves the estrogen receptor (ESR1), which induces PEBP1 through the kinase NME1. Our model was derived in the specific context of breast cancer, but the observed responses to drug treatments were consistent in other cell lines. We further validated some of the predicted regulatory links in the breast cancer cell line MCF7 experimentally and highlighted the points of uncertainty in our model. To summarize, our model was consistent with the observed changes in activity with drug perturbations. In particular, two pathways, including PEBP1, were highly responsive and would be likely targets for intervention.


Introduction
Several approaches for screening and identifying bioactive compounds exist and are being used to find effective remedies for cancer. Most strategies focus on quantifying the phenotypic effects of different biological models [1]. High-throughput technologies like massive parallel sequencing and mass spectrometry have the potential to augment these strategies by explaining the effect of panels of drugs on biological entities such as gene expression and protein abundance. Reverse causal reasoning enables the use

Knockdown of Metastasis Suppressors and Transcription Factors in MCF7
We used two datasets where each MSG and TF was knocked down in the MCF7 breast cancer cell line. The first is a dataset that consists of 77 microarray samples [28]. Seventeen MSGs were knocked down (siRNA) separately in MCF7 ( Table 1). The control set consisted of wild-type (0 h) and cells treated with scramble siRNA at 48 or 96 h. We obtained the processed data using GEOquery [29]. We mapped the probe intensities to gene symbols and applied differential expression between knockdown and control using limma [30]. We compiled a second dataset from the KnockTF database for the eleven TFs knockdown in the same cell line (Table 2) [31]. KnockTF curates gene expression profiles of TFs knockdown from the gene expression omnibus (GEO) and presents the results in the form of differential expression output (fold-change and p-values).

Pharmacological Perturbations of Breast Cancer Cell Lines
Two sources of publicly available drug perturbation data were used in this study. The first is growth inhibition (GRmax) data of drug perturbations (n = 35) in breast cancer cell lines (n = 71, including MCF7) [32]. Adenosine triphosphate content was measured as a proxy for the cell count in drug treated cell cultures relative to DSMO-treated controls. We included in the analysis only the drugs with at least one effective dose. This dataset was used as evidence for the efficacy of these drugs in producing phenotypic effects on cancer cells. The second data source is the library of integrated network-based cellular signatures (LINCS), which contains gene expression profiles of multiple cells under different types of perturbations, including compounds, gene overexpression, knockdown, or knockouts [33]. Several cell lines were profiled for gene expression using the L1000 technology, which measures the abundance of 1000 landmark genes and infers the expression of the rest of the genes. We used a subset of drug perturbations (n = 35, same set of drugs as above) in cancer cell lines (n = 67) that were also surveyed for drug sensitivity. We obtained the expression profiles of all genes (level 3 of the LINCS data) using Slinky [34].

Network Perturbation Amplitude (NPA)
Network perturbation amplitude (NPA) models the perturbation of a causal network as the changes in the expression of connected downstream nodes [4]. This method assumes that the downstream gene expression reflects the function of the entities in a biological pathway (e.g., the activity of a given TF). The input to this method is a two-layer network model of the pathway of interest and expression profiles of two or more conditions. First, the functional layer (backbone) encodes the causal relations between the entities of interest. The second is a transcription layer (terminal) with all the nodes downstream from each node in the backbone. We used the gene expression profiles from LINCS pharmacological perturbations to score the model. NPA is an R package that efficiently implements the method with the same name [35].

Measures of Agreement
We defined two notions of agreement between expectations and observations of the effects of drug perturbations on the metastasis network model. For every node (u), we created a graph of all n nodes connected to it by an edge (e) (subnetwork). Given a binarized perturbation coefficient of drug treatment, we formed an expectation by multiplying +1 or −1 with the direction of interaction +1 or −1 for every edge, u × e. We compared the expected sign to the observed perturbation coefficient of the downstream nodes (x ). The average for every subnetwork we termed a concordance rate. We similarly computed an agreement estimate for all possible directed graphs that lead to a particular node (path). The difference from above is that the upstream node is not fixed, but changes as we move one edge down the path. The average agreement between the observed and expected effect for every path we termed a coherence rate. The following formula shows how the agreement rates were calculated.
where x is a component connecting a node u to a downstream node by an edge e and x is the observed effect of drug perturbation on that downstream node in the component of the subnetworks (concordance) or paths (coherence).
In both cases, the agreement was expressed in terms of Cohen's κ [87]. This is a value between −1 (worse) and 1 (better) while taking into account the expected agreement by chance. The following shows how to calculate it.
Finally, we compared the probability distribution of the agreement measures to randomly generated values using the Kolmogorov-Smirnov (KS) test. D + is the maximum distance between the cumulative distribution functions (ECDF).

Reagents and Drugs
Reagents and drugs utilized in this study were purchased as follows: RPMI-1640 media (

Cell Culture and Drugs Treatment
MCF-7 breast cancer cells were cultured in RPMI-1640 supplemented with 10% Fetal Bovine Serum (GIBCO) and 100 µg/mL streptomycin and incubated in a 37 • C humidified atmosphere containing 5% CO 2 . For drug treatment, cells (3 × 10 5 cells) were plated in a 6-well plate and further incubated for 24 h. Then, cells were treated with indicated drugs: Epirubicin, Vorinostat, Methotrexate, Cisplatin, Sorafenib, and Imatinib at 1.5 µM concentration and collected for assays after 24 h of treatment.

RNA Extraction, and RT-qPCR
Total RNAs were extracted from MCF7 breast cancer cells using Trizol reagent according to the manufacture's manual. Prior to cDNA synthesis, RNA samples were treated with DNase I solution to remove trace amounts of DNA. The prepared RNA samples were then used as templates for reverse transcriptase to generate First-strand cDNA using Thermo Scientific RevertAid First Strand cDNA Synthesis Kit. The first-strand cDNA synthesis products were used directly in qPCR using the amfiSure qGreen Q-PCR Master Mix kit. Primers of target genes used in the assay are presented in Table 4. The expression of the five genes was quantified in the treated samples relative to the control gene GAPDH and the control condition of the DMSO treatment. ∆∆C t model was applied using the pcr R package [88]. Student t-test was used to compare the relative expression in each treatment to the control DMSO treated cells. p-values < 0.05 were considered significant. Experiments were performed in five or more replicates. Table 4. RT-qPCR primers.

Gene
Forward Primer Reverse Primer

Software Environment and Reproducibility
This analysis was performed in R and using Bioconductor packages [89,90]. The software environment was packaged and distributed as a Docker image (https://hub. docker.com/r/bcmslab/antimetastatic, accessed on 29 November 2021). The code to run the analysis and reproduce the figures and tables in this manuscript is available as open-source (GPL-3) (https://github.com/BCMSLab/antimetastatic, accessed on 29 November 2021).

A Workflow for Building a Network of the Metastasis Suppressors and Their Regulators
In this study, we built and evaluated a network model of MSGs and their transcriptional regulators in breast cancer ( Figure 1A). First, we identified PPIs of seventeen MSG and eleven TFs in the String database. Then, we vetted and supplemented every link with information from the relevant literature to determine the type of entity, direction of interaction, and evidence. We added a few interactions that were not present in the PPIs.
Since these causal links originated from experiments in different organisms and conditions, we further filtered the model to keep only the interactions relevant to the context of breast cancer. We obtained and analyzed the knockdown datasets of each of the input genes and regulators in the MCF7 cell line. We resolved the conflicts between the reported interactions and fold-change from the knockdown data by prioritizing the context of the interactions, the strength of the evidence, and the effect size ( Figure 1B). We used differential expression to identify possible links between the nodes that were not previously reported and add them to the network. We compiled the results were in a context-specific metastasis model using BEL. The full list of included regulatory links is reported in Table 3 into three categories: within MSGs, within TFs, and in between the two groups of gene products.

Figure 1.
A workflow for building, refining, and evaluating a breast cancer-specific network model of metastasis. (A) We queried the String database for protein-protein interactions for metastasis suppressor genes (MSG) and transcription factors (TF). Next, we surveyed the relevant literature to vet and supplement the interactions and express them as possible causal links in the biological expression language (BEL). Finally, we filtered and evaluated the network model using knockdown and drug perturbation datasets. (B) Conflicts between the data-driven and curated interactions were resolved by prioritizing the context, the strength of the evidence, and the effect sizes. We scored the metastasis model on drug perturbations using the network perturbation amplitude (NPA). We divided the network into smaller subgraphs and measured the agreement between the expected and the observed direction activity. We considered the consistency in (C) subnetworks of nodes connected to an upstream by one edge (concordance) and (D) paths connecting a particular node to its upstream by a sequence of edges (coherence).
Finally, we evaluated the context-specific model using gene expression and growth rate data in MCF7 treated with different compounds. The context-specific interactions served as a backbone to a transcript layer derived from the knockdown data to form a two-layer network. We scored the model using NPA under different conditions. We then compared the effect of the drugs on the entire network to their impact on the growth rate. Next, and given the expected effect of drugs on each node, we compared this predicted change on the downstream nodes to the observed using two different measures of agreement ( Figure 1C,D). Moreover, we constructed a reliable model of highly consistent paths for further exploration.

Identifying Possible Interactions of Metastasis Suppressors Using Text Mining of the Literature
We took as a starting point the known protein-protein interactions between seventeen MSGs (Table 1) and their eleven TFs ( Table 2). The initial query on the String database resulted in 93 interactions. Most interactions were based on text mining or the proximity of terms (gene products) in the text of one or more published study. We manually examined the provided evidence for each interaction. We also filtered, and supplemented with information from related studies. The interactions were coded in BEL and were used as the functional layer for the NPA analysis (Table 3). Each pair of interacting gene products were considered as a subject and an object. The type and direction of the relation between the two entities was represented according to the BEL standards.

Contextualizing Metastasis Suppressor Interactions in Breast Cancer Cells
Inevitably, conflicts arise between the relations suggested in the literature and the manifest changes of knocking down some genes on the expression of others. For example, not all interactions are relevant in the biology of breast cancer because they come from studies of different conditions. Additionally, not all interactions are reported in previous studies. We used a data-driven approach to specify the associations that are most relevant and augment the network with interactions based on changes in gene expression. The efficiency of the knockdown was verified by checking the expression of the corresponding genes (Figure 2, bottom). Most knocked down genes had a significant fold-change (log 2 FC < −1 and p-value < 0.05) when compared to the control cells. Similarly, a wide-ranging change (between −4 and 2 log 2 FC) in expression in the metastasis suppressors and regulators resulted from the knockdown, in most cases (Figure 2, top). We obtained expression profiles of seventeen metastasis suppressor genes (MSG) and eleven transcription factor (TF) knockdowns in the MCF7 cell line. We calculated the differential expression between the knockdown and control samples. Bottom, bars represent the fold-change (log 2 ) of the target gene (Self) between the control and knockdown. Top, box plots represent the distribution of the fold-change of (Other) MSGs and TFs in the corresponding knockdown condition.
We compared the observed fold-change to the directions of previously curated interactions. We kept the ones that had no evidence to the contrary as a result of the knock-down of the subject. We removed the links with a significant (absolute log 2 FC > 1 and p-value < 0.05) change in the opposite direction of change when the upstream node was knocked down. Finally, we coded the significant effects of the knockdown of one node (subject) on others (object) as an interaction between the mRNA to directly increase or decrease. The resulting network contained 78 edges between the input 35 unique nodes ( Table 3). The recorded interactions included relations between the coding RNA, protein, and protein activities. That direction of interaction was either positive (increase/directlyIncrease) or negative (decrease/directlyDecrease) ( Figure 3A).
TFs seem to have extensive interactions, reflecting the hierarchical nature of the regulator networks reported in other studies. Only a few nodes in the network had many edges, and most nodes had, on average, only two edges ( Figure 3B, top). The highly connected nodes were either TFs that received signals from other regulators, key nodes in signaling pathways, or effector nodes that are tightly regulated. Examples include FOS and ESR1 as TFs, MAP Kinase subunits as signaling nodes, or CDH1 as an effector. In addition to filtering in the relevant interactions, we used the knockdown dataset to build a transcript layer for the functional model. The numbers of nodes connected to those in the functional layer were more normally distributed. Similarly, the distribution of up and down-regulation is normal but down-regulated nodes were more frequent ( Figure 3B, bottom).

Evaluating the Metastasis Model Using Drug Perturbation Data
To evaluate the relevance of the metastasis network, we used a dataset of growth inhibition rates and gene expression profiles under drug perturbations. In this dataset, the growth rate was quantified in breast cancer cell lines (n = 71), including MCF7. Cell cultures were treated with several compounds (n = 35) and compared to the baseline replication rate. We included in the analysis only the drugs with at least one effective dose in terms of significantly (GRmax < 1) inhibiting cell growth ( Figure 4A). Drug responses in most cell lines were correlated with the response in MCF7 ( Figure 4B). NPA infers the activity of the nodes in the functional layer of the metastasis network using a transcript layer derived from the knockdown dataset. This method returns a perturbation coefficient that indicates the magnitude and direction of changes (−1, inhibition, and 1, activation) for the network as a whole and for each node in the graph ( Figure 4C). The perturbation coefficients ranged from −0.1 to 0.2, with most nodes being significantly (confidence intervals not including zero) perturbed at least once. The perturbation amplitudes of 67 cancer cell lines correlated well with the NPA value in MCF7 ( Figure 4D). Next, we divided the metastasis graph into smaller modules and evaluated the agreement between the expected direction of change and the observed coefficients for each perturbation. We considered every node as an upstream connected to one or more nodes by one edge (subnetwork). The effect of the drugs on the upstream node was binarized (activation > 0, repression < 0) and used to infer the expected impact on the downstream nodes. We compared the expectations to the observed perturbation coefficients of the downstream node after similar binarization (concordance; see Section 2).
Higher concordance supports the consistency of the suggested model. Cohen's κ was employed to evaluate the agreement between expectations and observation (κ > 0, better and κ < 0, worse), taking into account the chance agreement (50%). The average concordance of each drug was higher (D + = 2.9 and p-value < 0.0001 in KS test) than expected by chance alone (Figure 5A). The network accurately predicted the direction of drug perturbations. Averaging for each node, a high concordance (D + = 2.6 and p-value < 0.0001) was also observed ( Figure 5B). The correlation between the concordance of the subgraph and the effect of the drug on its source node, a potential source of bias, was very weak (r < −0.02 and p-value > 0.9). We observed higher than chance (D + = 2.47 and p-value < 0.0001) of consistency in other cell lines treated with the same set of drugs ( Figure 5C).

Constructing a Model of PEBP1 and Its Interaction with Other Metastasis Suppressors
To demonstrate the utility of this model, we used the same perturbation data to prioritize key nodes and consistent edges between them. We isolated all paths in the graph leading to PEPB1 (n = 11). Paths are all interactions in the directed graph that end with the node corresponding to a protein. We multiplied the binarized perturbation coefficients of the upstream nodes by the sign of the edges between the nodes in every path. We compared these expectations were the observed perturbation coefficients on the downstream nodes of the paths (coherence; see Section 2). Not all drugs produced coherent effects on PEBP1 pathways, but the agreement was bigger than to be by chance (D + = 2.44 and p-value < 0.0001 in KS test) ( Figure 6A). Most paths were coherent (D + = 1.45 and p-value < 0.0001) indicating the strong relevance of these pathways/interactions ( Figure 6B). The main sources of bias (upstream coefficients) were not strongly relevant (r < −0.1 and p-value > 0.7). The observed directions of interaction were also consistent (D + = 3.01 and p-value < 0.0001) with predictions across different cancer cell lines when treated with the same drugs ( Figure 6C).
We used the top five coherent (coherence > 0.6) PEBP1 paths to highlight interactions that are strongly relevant to breast cancer cell metastasis and are highly responsive to treatment with different drugs ( Figure 7A). We constructed a model of regulatory interactions between PEBP1 and its upstream regulators with the top five coherent paths represented by dash-lines ( Figure 7B). Interestingly, we found that the pathways flowed in two main directions to enhance or suppress PEBP1 activity-the upstream proteins of PEBP1 function as TFs/regulators. One involves RELA (p65) and SNAI1, which were previously reported to inhibit PEBP1 transcription. The other involves the estrogen receptor (ESR1), which induces PEBP1 through the kinase NME1, which to our knowledge has not been reported before. Figure 6. Coherence of expectations and observations in the paths to PEBP1. We calculated and binarized (1 for activation or −1 for repression) the perturbation coefficients of every node in the network. We then evaluated the agreement between the expected and observed direction of change in the paths connecting a PEBP1 to its upstream by a sequence of edges (coherence). First, we multiplied the drug's effect on the upstream nodes by the sign of the edge connecting it to the next node to form expectations. Next, we compared the expectations with the actual perturbation coefficients of the corresponding nodes. Negative Cohen's κ indicates worse and positive better agreement than expected by chance. Finally, we compared the probability distribution of the coherence to randomly generated values using Kolmogorov-Smirnov (KS) test. D + is the maximum distance between the cumulative distribution functions (ECDF). (A) A histogram of the average coherence for every drug.

Validating the Model Predictions
To test the validity of the predictive model of PEBP1/RKIP regulation, we selected several drugs that activate or repress PEBP1 and tested their effect on gene expression in MCF7. Indeed, four out of 6 drugs conformed to the expected direction of regulation. Epirubicin induced the expression of PEBP1 while Sorafenib, Cisplatin and Imatinib repressed it ( Figure 8A). One repressor (Sorafenib) seems to produce its effect on PEBP1 through either of the two suggested pathways ( Figure 8B). Sorafenib treatment induced RELA and SNAI1, which lowered PEBP1 and repressed ESR1 and NME1, which activated it. The effect size on NME1 was small and not significant at a cuttoff of p = 0.05. By contrast, Cisplatin and Imatinib repressed PEBP1 through one or the other regulatory pathways. Cisplatin induced RELA and SNAI1 which inhibited PEBP1 ( Figure 8D), while Imatinib repressed PEPB1 through ESR1 and NME1 inhibition ( Figure 8D). The activator, Epirubicin, induced the expression of PEBP1 through activating NME1 ( Figure 8E). Our model predicted the direction of regulation of PEBP1 and the pathway by which this regulation was achieved ( Figure 8F).

Discussion
We used text mining datasets and a manual literature search to extract evidence for possible interactions between several metastasis suppressors and their regulators. We filtered these interactions in the context of breast cancer using a knockdown dataset in the cell line MCF7. The resulting interactions were coded in BEL to build a functional metastasis model. Finally, we used a reverse causal reasoning approach to evaluate and prioritize these interactions and extract pathways that are most consistent with drug treatments that inhibit cell growth.
High-throughput technologies such as massive parallel sequencing and mass spectrometry produce simultaneous measurements on many biological entities. Interpreting these measurements requires expert biological knowledge. An approach to formalize the use of knowledge and data is to use causal reverse reasoning [4]. Here, we employ such an approach to infer the function of metastatic suppressors using curated causal links and the perturbation data. The implicit advantage of this approach is also to make claims about the role of gene products from measurements of gene expression. Castro and colleagues used a similar approach to construct a model of inflammatory bowel disease (IBD). They were able show that three pathways that contribute to the disease were regulated differently in the two subtypes of IBD [91]. Another study identified a mechanism the mediates the difference between the subtypes in response to environmental exposure such as cigarette smoke [92].
Studies based on text mining or manual curation alone can suffer from several drawbacks. Some of the identified interactions can only be relevant in a particular biological or pathological context. In addition, inconsistencies arise with changing the context of the interactions. Finally, inevitable biases exist in the literature toward studying specific pathways or models. We attempted to address these issues in our study with a number of strategies, which we will discuss next.
We took all possible interactions as defined by text mining of relevant studies and manual curation as a prior ( Figure 1A). Next, we evaluated the relevance and consistency of these links. Interactions that were not contradicted by the knockdown of their corresponding nodes are likely to be present in this context ( Figure 1B). This filtering step increases the likelihood of the reported links being relevant. Similarly, consistency in the direction of change with drug perturbations increases our confidence in the relevance of some interactions which we initially included. We implemented intuitive consistency measures for all possible subnetworks ( Figure 1C) and for selected paths ( Figure 1D) that include PEBP1 as the target node.
In addition to filtering out irrelevant and insignificant interactions, we used the knockdown dataset to identify regulatory links that were not previously reported in the literature. We inferred potential edges from significant changes in gene expression as a result of the knockdown of a given node (Figure 2). These links would augment the model by increasing the number of edges and mitigating biases arising from tendencies in the literature toward studying certain pathways and specific disease contexts. We suspected that larger effect sizes (perturbation amplitudes) could bias the agreement measures. However, neither in the case of subnetwork concordance nor the path coherence was this bias big enough to influence the results (Figures 5C and 6C).
The value of any predictive model is to make testable hypotheses. Here, we suggested two key pathways regulating PEBP1 (Figure 7). First, the inhibition of the zinc transcriptional repressor (SNAI1) results in PEBP1 suppression [75]. In contrast to PEBP1's role in metastasis, SNAI1 induced breast cancer EMT and metastasis by directly repressing the epithelial markers E-cadherin (CDH1) [74]. Martin and colleagues found SNAI1 to be overexpressed in invasive metastatic breast cancer compared to normal breast tissue [93]. Moreover, RELA (P65), which is one of five subunits of the Nuclear factor-kappaB family, positively regulated SNAI1 transcription. NF-κB binds directly to the SNAI1 promoter, resulting in increased SNAI1 transcription [94]. Overall, NF-κB, SNAI1, and PEBP1 form a feedback loop that regulates the development of metastasis and resistance to apoptosis [95]. MTA3 is another upstream node of SNAI1 and suppresses its transcription. Similarly, MTA3 represses the transcription of a series of EMT-promoting genes such as ZEB2 and N-cadherin [81]. Our model is consistent with the above-mentioned studies.
Another pathway to regulate PEBP1 in our model is through ESR1 and NME1. NME1 dramatically inhibits the NF-κB signaling pathway and alters the transcription of metastasisrelated genes, reducing metastatic pulmonary cells [96]. Additionally, NME1 regulates gene expression in breast cancer cells. Moreover, NME1-dependent genes have a prognostic value as they predicted survival in breast cancer patients [50]. Here, we propose for the first time an alternative link between the ESR1 and PEBP1 induction via the kinase NME1.
We experimentally tested the effect of several drugs on PEBP1/RKIP and its regulators. In addition to verifying the predicted interactions, we wanted to highlight the uncertainty associated with this model. The model had three points of uncertainty. The first is whether a given drug produces the predicted effect on the terminal node of interest. Four out of the six tested drugs conformed to the expectations. The second uncertainty is whether that effect is achieved through both, either or a third unspecified pathway. Sorafenib inhibited PEBP1 by activating the inhibitory pathway of RELA/SNAI1 and inhibiting the other activation pathway of ESR1/NME. By contrast, Cisplatin and Imatinib repressed PEBP1 through only one of the two pathways. The third uncertainty stems from the fact that the farther back in the pathway we travel, the less accurate the predictions are. For example, Epirubicin activated PEBP1 by inducing NME1, but the effect on ESR1 was not as expected.
The validity of this model is limited to how much the included interactions represent the underlying biology. Another limitation is the amount of gene expression signal captured by the knockdown and perturbation data. We carefully curated existing studies to have as many and as accurate as possible links between metastasis suppressors. This model is guaranteed to evolve with the rapidly accumulating biomedical literature. As there is no existing high-throughput method to measure the activity of diverse biological entities in parallel, we relied on a reverse engineering approach to infer the activity from gene expression changes in their downstream nodes. The reverse reasoning approach assumes that the changes in gene expression of a group of nodes reflect the activation or inhibition of their upstream proteins.
The regulatory links in the final model were derived from the literature and prioritized by our analysis. However, further validation of these interactions is necessary first to test their relevance and the reliability of this approach. We supported our observations by devising two notions of consistency, concordance, and coherence. The first of the two shows the general agreement between expectations and observations in the subgraphs and the response to drug treatment. On the other hand, coherence shows consistency in particular pathways and can be thought of as a measure of the relevance of specific links. We elected to use those two measures of agreement for ease of interpretation rather than more sophisticated notions of consistency. We carried out experimental validation and found that our model is most reliable in predicting the drug effect on a target node, but uncertainty increases as to through which pathway and how many upstream nodes are regulated.

Conclusions
In conclusion, we suggested that the metastasis suppressor PEBP1 is on the receiving end of two key regulatory mechanisms. One inhibitory pathway involves RELA (NFKB p65 subunit) and SNAI1, previously reported to interact with PEBP1. The other pathway involves the estrogen receptor (ESR1), which induces PEBP1 through the kinase NME1. Those two pathways were highly responsive to pharmacological perturbations and would be likely targets for intervention against metastatic breast cancer.
Author Contributions: M.A., conceptualization, methodology, software, formal analysis, visualization, and writing. T.H.L., conceptualization, methodology, data curation, data generation and writing. W.K., methodology, writing. D.R.K., conceptualization, writing, supervision, project administration, and funding acquisition. All authors have read and agreed to the published version of the manuscript.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript: