Logic-Based Modeling and Virtual Drug Screening for the Prediction of Novel Therapeutic Targets and Combination Regimens Against E2F1-driven Melanoma Progression

Skin melanoma presents increasing prevalence and poor outcomes. Progression to aggressive stages is characterized by overexpression of the transcription factor E2F1 and activation of downstream pro-metastatic gene regulatory networks (GRNs). Appropriate therapeutic manipulation of the E2F1-governed GRNs holds potential to prevent metastasis, however these networks entail complex feedback and feedforward regulatory motifs among various regulatory layers, which challenge the characterization of drug targetablemake it difficult to identify druggable components. To this end, computational approaches such as mathematical modeling and virtual screening are important tools to unveil the dynamics of these signaling networks and comprehensively identify critical components that could be further explored as therapeutic targets. Herein, we integrated a well-established E2F1-mediated epithelial-mesenchymal transition (EMT) map with transcriptomics data from E2F1-expressing melanoma cells to reconstruct a core regulatory network underlying aggressive melanoma. Using logic-based in silico perturbation experiments of a core regulatory network, we identifiedy that simultaneous perturbation of AKT1 and MDM2 drastically reduces EMT in metastatic melanoma. Using the structures of the two protein signatures along with virtual screening of lead-like compound library available in ZINC12 database, we identified a number of lead compounds that efficiently inhibit AKT1 and MDM2 without eliciting toxicities. We propose that these compounds could be taken into account in the design of novel therapeutic strategies for the management of aggressive melanoma. were identified using virtual screening of lead-like compound library available in ZINC12 database. Subsequent high-throughput virtual screening of drug library using the structures of the two protein signatures predicted a number of lead compounds that efficiently inhibit AKT1 and MDM2 without eliciting toxicities. These can be experimentally evaluated and further considered as new anti-melanoma metastatic agents, in monotherapies or combination regimens.


Introduction
Cutaneous melanoma arises from melanocytes and represents the deadliest form of skin cancer, with increasing prevalence. Once it becomes metastatic, the prognosis is very unfavorable. Melanoma formation is driven by mutations in the BRAF and NRAS oncogenes [1]. However, these oncogenic aberrations are early events in melanomagenesis that alone do not seem to be sufficient to drive metastasis [2]. Over the past years, we [3][4][5] and others [6] have demonstrated that in addition to these driver events, melanoma progression is catalyzed by the abundant expression of E2F1, a member of the E2F transcription factor family. Although this transcription factor activates tumor-suppressive pathways at early oncogenesis, upon disease progression unbalanced E2F1 activity is rewired to deregulated cancer networks that underlie hallmarks of metastatic progression, such as resistance to apoptosis, chemoresistance [3,7], neoangiogenesis [8], extravasation [9], epithelial-mesenchymal transition (EMT) [10,11], metabolic reprogramming [12], and genomic instability [13].
By integrating logic-based network modeling and gene expression profiles of cancer cell lines from E2F1driven tumors and patient cohorts displaying cancer aggressiveness, we identified tumor-type specific receptor signatures associated to EMT, where the combined action of highly expressed E2F1, TGFBR1 and FGFR1 triggers the most invasive phenotype [11,14]. Several other protein coding genes (PCGs), miRNA genes and lncRNA genes have been identified as constituents of E2F1-activated prometastatic gene regulatory networks (GRNs) [11][12][13]15,16]. Within the E2F1-governed GRNs, non-linear feedback and feedforward regulatory motifs are formed among various regulatory network layers, entailing protein-coding and non-coding RNA genes [10][11][12][13]15,16]. Such regulatory motifs, which are commonly encountered in cancer networks [17,18], induce a whole range of dynamic behaviors, thereby evading the use of conventional data analysis methods approved by the Food and Drug Administration (FDA) for the treatment of melanoma [20]. Unfortunately, responses to BRAF inhibitor monotherapy, although impressive and rapid, were usually transient. In most cases, this was due to development of resistance via reactivation of the mitogen activated protein kinase (MAPK) pathway. Combined BRAF and MEK inhibition addresses this MAPK-mediated mechanism of resistance and constitutes the current standard-of-care for targeted melanoma therapy. Regimens of BRAF plus MEK inhibitors achieve longer-lasting disease control and are better tolerated than BRAF inhibitor monotherapy, but a major concern is that resistance, although delayed, is eventually developed [21,22]. Likewise, the treatment response of patients with mutant NRAS-positive metastatic melanoma to MEK inhibitors is transient and shortlived [23]. In the context of the paradigm-changing advances in cancer immunotherapy, several next-generation immune-based formulations, such as the checkpoint inhibitors ipilimumab, pembrolizumab, and nivolumab, have received FDA approval for the indication of metastatic melanoma and ensure durable responses. However, they are linked with immune-related toxicities and pose limitations for use in patients with either an overactive (autoimmune disease patients) or a suppressed (organ transplant recipients) immune system [24]. In summary, targeted therapy shows a rapid response time and generally less toxic off-target effects. However, resistance can be developed. While immunotherapy offers increased durability of benefit in all patients irrespective of tumor genotype, this is also frequently associated with immune-related side-effects, especially when combinations of immunotherapeutic drugs are used [21,22]. It is therefore essential to develop both, effective and safe strategies, that specifically interfere with the complex melanoma networks. This should enrich the antimelanoma drug arsenal with more personalized therapeutic options. The past decade's clinical experience taught us that, in general, combination therapies may be superior in terms of efficacy and/or safety than monotherapies. Combining anticancer drugs is currently seen as the approach most likely to overcome singleagent resistance, to produce sustained clinical remissions via multi-targeting effects on distinct mechanisms of action, and to reduce unwanted side-effects by usage of lower drug doses [25,26]. In fact, the need for combinatory therapies is an inevitable consequence of the evolving nature of tumors. Clonal evolution is particularly active when tumors are under selective pressures due to medical treatments, thereby promoting resistance to therapy. Resistant cell clones often preexist (although undetectable) at the start of treatment, supporting the idea that early administration of combinatorial treatments stands a higher chance of eradicating such clones when their number is very low, before acquired resistance is overtly diagnosed. Simultaneous targeting the driver oncogenic mutations along with the expected secondary resistance may provide a significant advantage in survival compared with administration at relapse. However, ab initio combination therapies are challenging in the clinical oncology setting because of the narrow therapeutic window between tumor cells and host, which overall limits the number of agents that can be simultaneously tested [27]. With recent advances in high-throughput screening methods, a systematic evaluation of combinations among large collections of chemical compounds in vitro has become feasible. This typically requires large-scale experiments, in which the combinatorial responses are tested in various doses on cancer cell lines or patient-derived cells, resulting in dose-response matrices that capture the measured combination effects for every concentration pair in a particular sample [26]. For example, systematic screening of pairwise combinations of 104 FDA-approved oncology drugs has been performed in the NCI-60 panel of human tumor cell lines to produce a comprehensive platform for defining the pairs with enhanced therapeutic efficacy [25]. However, even with modern highthroughput instruments, experimental screening of drug combinations can become a 'Herculean task', as the number of conceivable drug combinations increases rapidly with the number of drugs under consideration. In addition, the inherent heterogeneity of cancer cells further challenges the experimental efforts, as the combinations need to be tested in various cell contexts and genomic backgrounds. Hence, computational methods are often recruited to guide the discovery of effective combinations that can be prioritized for further pre-clinical and clinical validation [26,28].
Herein, aided by in silico workflows, we sought to predict efficient and safe compounds that either alone or in combination prevent melanoma progression by specifically targeting components of the prometastatic E2F1-governed GRNs in melanoma. Using a comprehensive regulatory and functional map of E2F1 in tumor progression and metastasis [11] which contains different types of regulatory factors, including genes, proteins, microRNAs, or complexes, we identified a core regulatory network in melanoma [28,29]. The core regulatory network was subjected to logic-based modelling for detecting protein signatures which play an important role in interconnecting many of the responsive genes that are typically not identified through genebased differential expression analysis. Subsequent virtual screening, an increasingly technique that improves the speed and efficiency of the drug discovery process [30,31], was applied to find compound hits against protein signatures eliciting measurable biological responses and to estimate their safety profiles. Our approach predicted that MDM2 plus AKT1 inhibitors could be a promising combination that is worthwhile to be further investigated as a novel anti-metastatic regimen in high E2F1-expressing melanoma patients.

Establishment of computational pipelines for prediction of drug-targetable components of the E2F1governed prometastatic GRN in melanoma and in silico screening of different inhibitors, alone or in combination
Previously, we have designed a comprehensive regulatory and functional map of E2F1 in tumor progression and metastasis [11] which contains different types of regulatory factors (n= 879) including genes, proteins, microRNAs, or complexes; and interactions (m=2278) based on information retrieved from published literature and databases. The map was modularized into three E2F1 regulatory compartments such as extra-/intracellular receptor signaling, post-translational modifications, regulators of E2F1 activity; and seven functional compartments including cell cycle, quiescence, DNA repair, metabolism, apoptosis, survival, and angiogenesis/invasion. Using a computational pipeline, we used the map to unravel a tumor type-specific regulatory core and to predict receptor protein signatures in bladder and breast cancer underlying E2F1mediated EMT transition. The E2F1 map and the previously used workflow is applied to identify a key functional module (core regulatory network) in melanoma. This core regulatory network is composed of regulatory motifs and critical molecular interactions that drive phenotype switching in melanoma. The core regulatory network for melanoma was subjected to our computational pipeline to detect protein signatures that play an important role in interconnecting many of the responsive genes that are typically not identified through gene-based differential expression analysis.
Our workflow includes (i) network-based analysis of topological parameters to characterize the pattern of factors in a networked system, (ii) mapping of the gene expression profiles from melanoma cell lines onto the E2F1 map, (iii) network reduction via a multi-objective function [32] to provide motif ranking by user-defined weights in an iterative manner, (iv) boolean modeling to analyze and predict the protein signatures linked to aggressiveness in melanoma, (v) virtual screening to find compound hits against protein signatures that elicit measurable biological responses, and (vi) to predict ADME behaviors, pharmacokinetic parameters, and druglikeness of compounds. The workflow in Figure 1 was used for prioritizing therapeutic targets and also to screen potential drug candidates that may be validated for the advanced melanoma tumors.

Identification of the metastatic melanoma-specific core regulatory network
We used our previously published network-based approach to construct a melanoma-specific regulatory core from the comprehensive map of E2F1 [11]. Here, we utilized the workflow and the E2F1 map to identify key network motifs and critical molecular interactions that drive a highly invasive melanoma cell phenotype. To do this, we have used the data extracted from the E2F1 map and identified important network motifs by calculation of topological and non-topological properties of each node (Supplementary file 1).
The motifs were prioritized based on a multi-objective optimization function where the function contains parameters accounting for each property. Weights are assigned to each property in terms of their importance in a user-defined manner and then, the motifs are ranked according to the value of the objective function. The top ten high-scored motifs were selected from each weighting scenario (Supplementary file 1). Finally, we merged all the top ranked motifs to obtain a melanoma-specific regulatory core. We expanded the regulatory core by adding receptor proteins which are the first neighbors of ranked motif nodes in the E2F1 map. We also added four well-known markers CDH1, VIM, ZEB1, and SNAI1 [10] in the core to measure the EMT response ( Figure   2a).

Boolean modelling of the melanoma-specific core regulatory network
We encoded the core-regulatory network into a Boolean model for stimulus-response and perturbation analyses.
Stimulus-response analysis was used to identify the effect of up/down expressed receptors on the EMT phenotype, and perturbation analysis predicted potential drug target that can bring the phenotype to lower possible level. In Boolean model, state a node is represented into one of two possible states i.e., 0 (OFF, inactive) or 1 (ON, active) [33]. The regulatory relationships between upstream nodes (i.e., sources) to downstream nodes (i.e., targets) are encoded into Boolean functions using logical gates 'NOT', 'OR' and 'AND'. Further, we calibrated the Boolean functions with fold-change (FC) expression data [11] of publicly available dataset GSE46517 [34] from Gene Expression Omnibus (GEO).
To evaluate the input-output behavior, we divided the model into three layers: 1) Input layer, containing receptor molecules, 2) regulatory layer, comprising nodes constituting a core-regulatory network, and 3) Output layer, including EMT phenotype ( Figure 2b). Input layer was initialized with FC expression profile i.e., a node with negative FC was represented by a state 0 and a node with positive FC was represented by a state 1.

In silico perturbation simulations using Boolean modeling
Model was simulated with initial values derived from the expression profile of input nodes and confirmed that the logical state of nodes in a regulatory layer represents the data (for details see Supplementary file 2). The EMT phenotype was regulated by nodes ZEB1, CDH1, VIM and SNAI1 [10], and represented by 5 ordinal level ranking from 0 (minimum) to 4 (maximum): For the initial condition, model simulations result in EMT of level 3, where ZEB1 and SNAI1 are active and CDH1 is inactive (see Table 1a). Further, we performed perturbation analysis of all nodes (except ZEB1, CDH1, VIM and SNAI1) in the regulatory layer of the model to bring EMT from level 3 to minimum level. We identified that for a single perturbation (in this case inhibition) of MDM2 or MIR25, EMT can be reduced to level 1 (see Table 1b). CDH1 is activated upon inhibition of MDM2 which inhibit EMT as well as inhibit CTNNB1 which subsequently inhibit SNAI1 [35,36] to further reduce EMT. Similar effect was observed upon inhibition of MIR25 [37,38]. On the other hand, single perturbation (in this case activation) of AKT1 can increase the EMT to the highest level 4.  indirectly activate the EMT by downregulating another hallmark protein CDH1. We investigated the expression profiles of AKT1 and MDM2 and their impact on melanoma patient survival using Kaplan-Meier curve ( Figure   3) using TCGA melanoma SKCM dataset (https://portal.gdc.cancer.gov/projects/TCGA-SKCM). We found that higher expression of both AKT1 and MDM2 resulted in poor patient survival. These observations also confirm that the Boolean model simulation were successful in predicting potential proteins that may be targeted for the treatment of metastatic melanoma.

Screening of small molecule inhibitors to block protein signatures
To identify compound drugs that are most likely to bind to AKT1 and MDM2 protein signatures, molecular docking was performed with the filtered library compounds (Supplementary file 3). The information about the active sites of proteins is retrieved from the literature and PDB database. More specifically, for AKT1 we performed screening against the kinase domain (150-408) which is previously selected to identify ATPcompetitive inhibitors [39,40] (Figure 4a). For MDM2, many recent studies indicate that its overexpression and subsequent deactivation of p53 result in failure of apoptosis and cancer cell survival [41][42][43]. We investigated the p53-Mdm2 interaction surface which is ~700 Å 2 . This druggable pocket of MDM2 where p53 binds provides a great opportunity for compound inhibitors to disrupt p53-MDM2 interaction [44] (Figure 4b). Three hits namely, ZINC000043178353 [45], ZINC000040429080 [46], and ZINC000043202934 [47] are reported as selective AKT1 inhibitors with IC50 values of 0.5, 0.03, and 1.0 nM respectively, and displayed potency against AKT1, AKT2, and AKT3 within cells [40,47]. However, we found two novel hits (ZINC000001491367: binding energy -12.6 kcal/mol; ZINC000003939645: binding energy -12.5 kcal/mol) that were not investigated as AKT1 inhibitors before. Their binding energies towards the ATP binding pocket of AKT1 is comparable to ATP competitive inhibitors [40]. These two novel hits are reportedly inhibitors of CDKs and their cyclin partners, particularly CDK7/cyclin H and CDK2/cyclin E which are often deregulated in cancer. Both hits showed a considerable activity when compared with Seliciclib, a drug in phase II clinical trial for the treatment of cancer [48]. Interestingly, the number of interactions that strongly bind these two compounds into the cavity was also reasonably high as compared to other hits ( Figure 5).
To determine the difference between the binding mode of novel hits to AKT1 from that of previously known inhibitors, a comparative analysis was performed. It was found that the common interacting residues in all inhibitors were GLU234 and ASP292. The amino acid residue GLU234 of the protein backbone is necessary for the AKT1 biological activity and this interaction was found with most of the previously known ATPcompetitive kinase inhibitors [49,50]. The second set of electrostatic interactions and hydrogen bonds to ASP292 in AKT1 is critical because this position is typically occupied by a divalent cation (Mg2+) bound to ATP [50]. Other common amino acid residues for all the five hits were PHE161 and LYS179; which indicates that the top screening hits ZINC000001491367 and ZINC000003939645 are reliable and promising for further evaluation.  with key residues and particularly, hydrophobic interaction with residues VAL93 and LEU54 can be seen in all five hits ( Figure 6).
The first compound hit ZINC000000537755 (Fluspirilene), exhibits a convincing binding mode into the MDM2 pocket. Fluspirilene is reported to show anti-proliferative activity at 10 µm in the NCI60 tumor cell line [53].
The second hit ZINC000169352550 is a compound containing morpholinone which is highly potent and

ADME/pharmacokinetic predictions and drug-likeness Bioavailability:
The bioavailability radar plots ( Figure 7) show a rapid appraisal of drug-likeness based on the physicochemical properties of the lead molecules. In the graphical output, the radar area (pink color) is in the optimal range for compound hits ZINC000001491367, ZINC000003939645, ZINC000043178353, ZINC000000537755, ZINC000095605306, and ZINC000084689539 giving information that the hits are falling entirely within the physicochemical range on each axis and could be considered drug-like. However, the compounds hits ZINC000040429080, ZINC000043202934, ZINC000000537755, and ZINC000096286451 are predicted to be not orally bioavailable because they are too polar and fraction Csp3 (in-saturation) is too high. Bioaccumulation: As shown in Table 2, the compound hits ZINC000001491367, ZINC000003939645, ZINC000095605306, ZINC000084689539, and ZINC000096286451 did not exhibit in silico inhibition of CYP2C9 and CYP2D6, members of the drug-metabolizing cytochrome P450 family of enzymes [49,56]. Though, other hits are interacting with CYP isoenzymes which could lead to bioaccumulation of compounds and toxicity.

Aqueous solubility and gastrointestinal absorption:
The aqueous solubility for all hits is estimated to be moderate except ZINC000000537755 which is predicted to be poorly soluble. Prediction of passive gastrointestinal absorption (GIA) was high for all the hits except ZINC000043202934 and it is based on the Intestinal Estimated permeation model [57]. It is observed that hits ZINC000001491367 and ZINC000003939645 are non-substrate to P-glycoprotein (multidrug resistance protein in the cell membrane) [58], suggesting that they are likely to have high intestinal absorption and bioavailability.

Discussion
A common approach to anticancer drug development has been based on a workflow, whereby molecules that are designed from scratch, to specifically interfere with a certain pathway, are anticipated to target and eradicate tumors in a highly selective manner, analogous to the "lock-and-key" specificity, hence maximizing efficacy and minimizing side effects [60]. Despite their promising results in the preclinical setting, the majority of innovative drugs are proven insufficient or suboptimal when administered in clinical patients, thereby leading to unacceptably low success rates of clinical trials [61,62]. The high failure rate of this approach is the consequence of several unpredictable parameters, mainly: (a) the individual genetic background of cancer patients, which limits the therapeutic benefits only to specific patient subpopulations and necessitates treatment personalization [22]; (b) the fact that cancer-related genes are highly interconnected and regulate each other through complex loops from different pathways [63][64][65]; (c) the inherent ability of tumors to adapt and evolve, which catalyzes acquisition of resistance to therapies, especially monotherapies [27]. To address these challenges, computational methodologies including, but not limited to, algorithms and machine learning tools, are now being increasingly recruited in many drug discovery programs. For example, computational approaches that 'dock' small molecules into the structures of macromolecular targets and 'score' their potential complementarity to binding sites are widely used in hit identification and lead optimization and are currently reforming the pharmacopeia landscape [66]. This approach allows for fast and comprehensive screening of the efficacy and safety profiles of a high number of leads, in the context of a particular cancer type. Prioritization of the top-resulting leads or combinations thereof could subsequently facilitate faster introduction to clinical trials and significantly reduce the costs for drug development.
Having in mind that metastasis is linked with activation of E2F1-governed GRNs, we applied a transcriptomics-aided bioinformatics workflow, followed by virtual drug screening to comprehensively characterize novel therapeutic targets in melanoma and predict their corresponding drug inhibitors. Due to the documented ability of targeted drugs to show superior safety and efficacy in combination schemes [22], we were particularly interested on drugs that can perturb these prometastatic GRNs when used simultaneously.
Using a well-established E2F1 map [11], we derived a set of three-node FBLs (n = 44) and used a ranking scheme that applies a weighted multi-objective function integrating topological and non-topological properties of each node. Topological properties such as node degree (number of edges connected to the node) are known for their importance in network organization and playing as central hubs in orchestrating molecular connections [67]. It is reported that cancer-associated proteins have large betweenness centrality as they control the communication between different components of a network [68]. Among non-topological properties, we have calculated the involvement of the motif constituents in the disease pathway, the gene prioritization score, and average Log2 fold change for each motif based on the change in expression values of each node from noninvasive to invasive phenotypes derived from in vitro experiments. Since the network was originally constructed around E2F1, the topological properties for some nodes are expected to be higher than other nodes. Therefore, to give equal importance to all nodes, we used different weighting scenarios in the multi-objective optimization function to avoid biases and ranked motifs accordingly. The top ranked motifs are merged to understand their combined effect on the regulation of EMT in melanoma. We further expanded the regulatory core network by adding receptor proteins the first neighbours of the ranked motif nodes and four marker proteins and their direct connections from the E2F1 map. Receptor proteins work as determinative factors and markers proteins are required to measure the EMT response. We developed a three-layered logic-based model of the regulatory core consisting of an input layer, a regulatory layer, and an output layer. We analyzed the regulatory core by using boolean logic for the input and regulatory layers, and multi-valued logic for the output layer which allows us to assess the combined effect of various network components on the EMT phenotype. Our model simulations identified two protein signatures AKT1 and MDM2 as potential drivers of EMT in melanoma. Further virtual drug screening with particular emphasis on the prediction of compounds with minimal toxicities revealed that AKT1 and MDM2 inhibitors, either alone or in combination with each other, can efficiently and safely suppress E2F1-driven invasion in melanoma.  [70]. Upregulation of the PI3K-AKT pathway is a critical event during the early and late evolution of resistance to MAPK pathway inhibition [71]. It has been proposed that AKTi combined with BRAFi-based therapy may benefit patients with tumors harboring BRAF mutations along with PTEN deletions or AKT mutations [72]. In agreement with these studies, our analysis highlighted inhibition of AKT1 as an attractive strategy for preventing EMT-driven metastatic progression of melanomas with a high-E2F1 content.
In addition to the PI3K/AKT pathway, the manipulation of the p53-controlled pathways is emerging as an alternative therapeutic option with a potential to overcome the suboptimal response rates for MAPKtargeting therapies [73]. The TP53 gene represents a well-established player in carcinogenesis [73]. Over 50% of all tumors carry p53 loss-of-function mutations, leading to synthesis of functionally impaired protein products which are unable to transactivate genes that induce cell cycle arrest and apoptosis in response to oncogenic stress [74]. Tumors maintaining wild-type p53 exhibit other types of downstream p53 inactivation, such as hyperactivation of its endogenous post-repressor mouse double minute 2 (MDM2), an ubiquitin ligase that catalyzes p53 degradation and inactivation [74]. Most melanomas retain a wild-type p53 and exhibit low TP53 mutation rates, but instead show frequent inactivation of the cyclin-dependent kinase inhibitor 2A (CDKN2A), which eventually lead to MDM2 upregulation and subsequent p53 inhibition [75]. As a result, small-molecule inhibitors that block the p53-MDM2 interaction have been pursued as a new cancer therapeutic strategy for restoring the tumor-suppressive function of p53 in TP53 wild-type tumors [76]. Shattuck-Brandt et al. [77] recently highlighted the therapeutic efficacy of MDM2 inhibition against TP53WT melanomas with either a wild-type or a mutant BRAF background. They showed that a MDM2 antagonist (namely  alone or in combination with BRAF and/or MEK inhibitors can inhibit tumor growth in patient-derived xenografts (PDX) from 15 patients with melanoma by suppressing p53 degradation. MDM2 inhibitor monotherapy was effective against BRAFV600WT tumors, while a combination of KRT-232 plus BRAF/MEK inhibitors exhibited a synergistic effect on BRAFV600E mutant PDXs [77]. In a similar note, our study revealed that a number of small-molecule MDM2 inhibitors show a potential for preventing E2F1-driven metastatic progression.
Overall, our approach predicted AKT1 and MDM2 inhibitors as promising anti-melanoma drugs. If combined with each other or with the standard-of-care regimens for melanoma, such as BRAF and MEK inhibitors, these substances could offer appealing alternative therapeutic strategies, potentially overcoming therapeutic resistance and improving disease-free survival of melanoma patients. Future experiments are essential to confirm the metastasis-preventing potential of such combinations in melanoma animal models.
Given that combinations of targeted therapeutics and immunotherapeutics hold a potential to produce durable responses to clinical patients [78], another significant question that is worth further investigation is whether AKT1 and/or MDM2 inhibitors can effectively synergize with checkpoint inhibitors to improve patient survival.
Melanoma is a highly heterogeneous and dynamically evolving cancer type. The increased intratumoral cell diversity can accelerate somatic evolution, because a tumor consisting of a genetically heterogeneous cell population has more possibilities to respond to microenvironmental changes, to evolve, and to spread [79]. It is noteworthy that melanomas can, in several cases, evolve through unorthodox pathways, via the same genes that can successfully inhibit melanoma proliferation. For example, while tyrosinase inhibition is seen as an approach to successfully target proliferative melanoma cells, its loss can trigger EMT-mediated melanoma progression [80]. In a similar manner, stabilization of the tumor-suppressor p53 can drive therapeutic resistance and it was recently suggested that inhibiting rather than activating wild-type p53 may sensitize previously resistant metastatic melanoma cells to therapy [81]. Advanced computational methods, such as artificial intelligence and machine learning are anticipated to shed more light on the roles of target genes and unveil the spatiotemporal complexity of networks underlying metastasis [82] towards developing personalized therapies.

Network analysis and motif identification
The Cytoscape version of the E2F1 map was downloaded from https://sourceforge.net/projects/e2f1map/files and converted into a format suitable for Cytoscape plugin NetDS v3.0 [17]. The purpose of this was to identify important nodes and network motifs in the network. The loop length was set to three nodes and feedback motifs (n = 444) were retrieved. We then used the Cytoscape plugin NetworkAnalyzer to evaluate the topological properties of nodes [83]. More specifically, we calculated the average number of neighbours for each node in the network (degree) [84] and the density of connections among the neighbours of a node (betweenness centrality) [85] to understand the overall organization of the network. Among non-topological properties, we calculated the number of nodes in a motif involved in the KEGG melanoma pathway (KEGG: 05218), and a prioritization score for each gene from the web resource DISEASES [86] shown in Supplementary file 1. Discovery Rate). Only targets displaying a minimum twofold induction or reduction (P < .05) by E2F1 knockdown were included for clustering [3].

Motif prioritization
The regulatory motifs were prioritized using a ranking score for each motifs considering key topological and non-topological properties with respect to the relevance for the melanoma phenotype. The motif ranking score is calculated using Eq. (1).
The equation uses a multi-objective function which is normalized to the maximum property value under consideration. We used a ranking scheme that is previously developed [11] by assigning different weights to various topological and non-topological parameters. In particular, the weights to two topological parameters (node degree ⟨ ⟩ and betweenness centrality ⟨ ⟩) was divided to half for avoiding over-emphasizes topological properties and assigning equal weighting factors 2 − 4 to give equal importance to other properties (disease pathway association ⟨ ⟩, gene prioritization score ⟨ ⟩, Log2 fold change ⟨| |⟩) in motif prioritization. The equation generates a ranking score for each motif (1…n) depending on the sets of values chosen for the weighting scenarios (1 to13) shown in Supplementary file 1. Later, top 10 motifs were selected from each of the weighting scenario (13•10=130 motifs). Furthermore, unique set of motifs were identified and processed for the construction of melanoma-specific core regulatory network. The optimization of multiobjective function is discussed in detail [11].
All the top ranked motif identified in the previous steps were merged to create a regulatory core. Additionally, we also considered receptor proteins as critical factors determining the EMT phenotype and directly interacting/ regulating nodes present in the top-ranked motifs. In total, we found and included ten receptor proteins (AR,   ESR1, FGFR1, FLT4, NR2F2, NR4A1, TGFBR1, TGFBR2, THRA, and THRB) into the regulatory core. These receptor proteins are the first neighbors of ranked motif nodes and present in the E2F1 map. In addition, we added four EMT marker proteins (CDH1, VIM, ZEB1, and SNAI1) and direct connections with motif nodes (Supplementary file 1) in our regulatory core.

Logic-based modeling to derive protein signatures
To identify protein signatures in the regulatory core, the network is translated into a logic-based model and in silico perturbation experiments were performed in the software tool CellNetAnalyzer [87]. For this, we derived Boolean rules for the input (receptor proteins) layer and propagation of signals from the input layer to the output layer through the nodes present in the regulatory layer. The network is simulated to determine the impact of the input layer vectors on the EMT phenotype (output layer). We performed single and double perturbation experiments iteratively for the initial conditions that is determined through the additional publicly available gene expression dataset (GSE46517) from Gene Expression Omnibus (GEO). The perturbation experiments were performed by changing the Boolean state of each node alone and in combination to other nodes in the regulatory layer to see the impact on the invasiveness. Those node(s) which upon inhibition change the EMT to minimum level or upon activation to maximum level are further evaluated as effective protein signatures associated with EMT transition in melanoma (Supplementary file 2).

Virtual screening of drugs
For virtual drug screening was performed as follows:

(ii) Structure preparation
The crystal structure of protein signatures AKT1 (PDB: 3OCB) and MDM2 (PDB: 3JZK) were downloaded from the RCSB Protein Data Bank (https://www.rcsb.org/). Proteins were pre-processed by removal of heteroatoms, adding polar hydrogens, and gasteiger charges using the AutoDock Vina [88]. Further, the coordinates of the active site residues were determined.

(iii) Molecular docking study
Virtual screening was carried out in PyRx v 0.8 (AutoDock Vina-based) screening tool [89]. The library compounds were first imported as SDF files in open babel of PyRx and further energy minimization of all the library compounds was performed followed by conversion into PDBQT format files. Later, a gird box was designed to cover the binding site residues within the protein signatures and then the library compounds were subjected to docking against AKT1 and MDM2. At each step, the energy of interaction of compound and protein was evaluated using binding energy (Kcal/mol) value (Supplementary file 3). The docked complexes and graphical visualization were done in DS Visualizer [90].

(iv) ADMET risk and pharmacokinetic prediction
The docked compounds were sorted based on binding energy (kcal/mol) and further filtered by computing a pool of ADMET risk, pharmacokinetics, drug-likeness, and medicinal chemistry friendliness prediction using the SwissADME server (http://www.swissadme.ch/) [57] (Supplementary file 3).

Conclusions
Cancer is a disease where multiple pathways are dysregulated, and its development and progression involve both independent and overlapping molecular targets. Advanced computational methods can unravel the properties of cancer-related proteins and their interactions in the molecular networks and enable designing of next generation targeted therapeutics. With the computational pipeline used in this study, we were successful in the identification of key protein signatures that derive melanoma metastasis phenotypes using in silico perturbation experiments. Using the virtual screening of lead compounds library, we identified key compounds that bind to AKT1 and MDM2 and suppress their metastatic activity. Some of the top hits were already investigated as potential inhibitors of the identified protein signatures. However, we also find novel hits, some of them were investigated for other caner types and can be further investigated for their potential to check melanoma metastasis. Among the top hits-ZINC000001491367, ZINC000003939645 are predicted as potential inhibitors for AKT1 and ZINC000095605306, ZINC000084689539 for MDM2 inhibition. These compound hits would facilitate the discovery and development of effective inhibitors for clinical use in melanoma metastasis.

Supplementary Materials:
Supplementary file 1: Four sheets including topological and non-topological properties of three-node FBLs and prioritization score; the top ranked FBLs identified in melanoma; weighting scenarios for motif prioritization; and regulatory core interactions in melanoma.