Logic-Based Modeling and Drug Repurposing for the Prediction of Novel Therapeutic Targets and Combination Regimens against E2F1-Driven Melanoma Progression

1. 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 unfavourable. Melanoma formation is driven by mutations in the BRAF and NRAS oncogenes [1]. However, these oncogenic aberrations are early events in melanoma genesis that alone do not seem to be sufficient to drive metastasis[2]. Over the past years, we [3,4,5] and others [6,7] 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,8], neoangiogenesis [9], extravasation [7], 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 E2F1-driven tumors and patient cohorts displaying cancer aggressiveness, we identified tumor-type specific receptor signatures associated with EMT, where the combined action of highly expressed E2F1, TGFBR1, and FGFR1 triggers the most invasive phenotype [11,14]. Several other protein-coding genes, miRNA genes, and lncRNA genes have been identified as constituents of E2F1-activated prometastatic 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,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 [19]. As a result, the prediction of potential therapeutic targets within this network requires the aid of systems biology-based and computational methods.

Uncovering major epigenetic features and the immune contexture of melanoma has catalyzed the development of anti-melanoma therapies within less than two decades. In 2004, no systemic therapies for melanoma had been shown to provide a survival benefit. Now, at least four regimens of targeted therapy and three for immunotherapy improve overall survival and disease-free survival, with each modality presenting distinct benefits and limitations. Particularly in 2011, vemurafenib became the first BRAF-targeted therapy 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 the 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 long-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 short-lived [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 anti-melanoma 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 single-agent 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 pre-exist (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 of 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 [26,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 [25,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 high throughput 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 validations [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 modeling for detecting protein signatures which play an important role in interconnecting many of the responsive genes that are typically not identified through gene-based differential expression analysis. The subsequent virtual screening which is a major contributor to computer-aided drug design (CADD) and drug repurposing concept, are increasingly popular techniques that improve the speed and efficiency of the drug discovery process [30,31], was applied to find FDA-approved drugs against protein signatures. This combined approach allows us to take advantage of existing safety profiles and established pharmacokinetic properties of approved drugs [32]. The selected protein-drug complexes were subjected to molecular dynamics simulation studies for investigating the stability of the complex formed. Our approach supports the role of AKT1 and MDM2 protein signatures as drivers of EMT in melanoma cancer and suggests 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.

2. Results

2.1. Establishment of a computational pipeline for the prediction of drug-targetable components of the E2F1-governed 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 E2F1-mediated EMT transition. The E2F1 map and the previously used workflow were 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 [34] 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) structure-based virtual screening and molecular dynamics simulation (MDS) studies to find repurposed drugs against protein signatures that elicit measurable biological responses, and (vi) to predict ADMET behaviors and pharmacokinetic parameters of candidate drugs. The workflow in Figure 1 was used for prioritizing therapeutic targets and also to screen potential candidate drugs that may be validated for advanced melanoma tumors.

2.2. 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 calculating of topological and non-topological properties of each node (Figure 1A) (additional file 1a). 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 (additional file 1b-c). 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 network (additional file 1d) to measure the EMT response (Figure 1B).

2.3. Boolean modeling 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 a potential drug target that can bring the phenotype to the lower possible level. In the boolean model, the state of a node is represented with two possible conditions 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 [33] of the publicly available dataset GSE46517 [34] from Gene Expression Omnibus (GEO).

To evaluate the input-output behavior, we divided the model into three layers: 1) the input layer, containing receptor molecules, 2) the regulatory layer, comprising nodes constituting a core-regulatory network, and 3) the output layer, including EMT phenotype (Figure 1B). The input layer was initialized with an 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.

2.4. In silico perturbation simulations using boolean modeling

The 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 additional file 2a-e). The EMT phenotype was regulated by nodes ZEB1, CDH1, VIM, and SNAI1 [10], and represented by 5 ordinal levels ranging from 0 (minimum) to 4 (maximum):

EMT = ZEB1 + NOT (CDH1) + VIM + SNAI1

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 a minimum level. We identified that for a single perturbation (in this case inhibition) of MDM2 or MIR25, EMT can be reduced to level 1 (Table 1). CDH1 is activated upon inhibition of MDM2 which inhibits EMT as well as inhibits CTNNB1 which subsequently inhibits SNAI1 [35,36] to further reduce EMT. A similar effect was observed upon inhibition of MIR25 [37,38]. On the other hand, a single perturbation (in this case activation) of AKT1 can increase the EMT to the highest level 4.

2.5. Assessment of protein signatures identified through boolean modeling

Our boolean model simulations suggested two key proteins AKT1 and MDM2 that upon inhibition can bring the EMT from level 3 to 1. Interestingly, AKT1 directly activates the VIM, a key marker for EMT. AKT1 also activates MDM2 which interacts with p53 to regulate the immune axis in metastatic melanoma. MDM2 also indirectly activates 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 the Kaplan-Meier curve (Figure 2) using the TCGA melanoma SKCM dataset. We found that higher expression of both AKT1 and MDM2 resulted in poor patient survival. These observations also confirm that the boolean model simulation was successful in predicting potential proteins that may be targeted for the treatment of metastatic melanoma.

2.6. Screening of FDA-approved drugs to block protein signatures

To identify drugs that are most likely to bind to AKT1 and MDM2 protein signatures, molecular docking was performed with the FDA-approved drug library from the ZINC database (additional file 3a-c). 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 ATP-competitive inhibitors [39,40] (Figure 1C). 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 1C).

From the achieved docked poses for each ligand-protein complex, we selected the pose with the lowest energy value (out of 9 conformations) and compared it with the crystallographic pose, RMSD ≤ 2.0 Å. Further, the best predicted binding mode and the corresponding binding affinity (in kcal/mol) are selected for each complex. In case of AKT1 Figure 3A (i-v), docking analyses revealed that the candidate drugs are packed against the residues LEU156, GLY157, PHE161, VAL164, ALA177, LYS179, GLU191, HIS194, GLU198, TYR229, ALA230, GLU234, ASP274, ASN279, MET281, ASP292, GLY294, LEU295, TYR437, PHE438, ASP439, and PHE442 and was stabilized by the hydrogen bonds, electrostatic, hydrophobic, and van der walls interactions. The best binding affinity of AKT1 is obtained with Cialis (-11.1 kcal/mol), followed by Paliperidone (-10.8 kcal/mol), Cobimetinib (-10.6 kcal/mol), Troglitazone (-10.5 kcal/mol), and Sertindole (-10.3 kcal/mol). The binding affinities and the number of interactions of these candidate drugs towards the ATP binding pocket of AKT1 are comparable to ATP competitive inhibitors [40]. Particularly, GLU234 of the protein backbone is necessary for the AKT1 biological activity. Secondly, the electrostatic interactions and hydrogen bonds to ASP292 in AKT1 are critical because this position is typically occupied by a divalent cation (Mg2+) bound to ATP [45]. Both of the residues (GLU234 and ASP292) of AKT1 were found to bind with previously known ATP-competitive kinase inhibitors [45,46] and also found with our candidate drugs. Other active site residues could be seen in binding were ALA230, GLU228, GLU278, ASN279, PHE161, and LYS179. Cialis/Tadalafil is a US FDA-approved drug for the treatment of pulmonary arterial hypertension [47]. Also, the drug shows biologic activity in human melanoma and pilot trial studies reported improved clinical outcomes in patients with head and neck cancer (HNSCC) and metastatic melanoma [48,49]. The next drug Paliperidone is an antipsychotic drug that reportedly inhibits glioblastoma growth in mice [50]. However, existing evidence is conflicting with its role as the drug increases prolactin levels, which might increase the risk of breast cancer [51]. Further in the list, Cobimetinib/GDC-0973 is FDA-approved for the treatment of patients with metastatic melanoma with a BRAF V600E or V600K mutation, in combination with vemurafenib [52]. Also, in combination with (i) chemotherapy and (ii) Niraparib, with or without atezolizumab are used as a treatment for patients with breast cancer and with advanced platinum-sensitive ovarian cancer (red) respectively [53]. The next drug hit is Troglitazone, which is a type II diabetic drug and, in clinical trials used in combination with lovastatin and their cotreatment was found to induce cell cycle arrest at the G0/G1 phase in Anaplastic thyroid cancer [54]. On the contrary, there have been concerns as it has detrimental side effects such as causing liver toxicity [55]. The last drug hit was Sertindole which exhibits antiproliferative activities in breast cancer with a potential application for the treatment of breast-to-brain metastases [56] and by inhibiting the STAT3 signaling pathway in human gastric cancer cells [57].

In the case of MDM2 Figure 3B (i-v), the molecular docking was performed into lining residues of this pocket containing amino acids (LEU54, LEU57, ILE61, MET62, TYR67, GLN72, VAL75, PHE86, PHE91, VAL93, HIS96, ILE99, TYE100, and ILE101). These residues form a hydrophobic cavity on the MDM2 protein structure and are potentially occupied by known inhibitors [58,59]. The best binding affinity of MDM2 is obtained with Finasteride (-10.7 kcal/mol), followed by Cobimetinib (-10.6 kcal/mol), Troglitazone (-10.4 kcal/mol), Loratadine (-9.8 kcal/mol), and Drospirenone (-9.7 kcal/mol). All top five candidates exhibit convincing binding mode into the druggable pocket of MDM2, and specifically hydrophobic interactions with key residues (VAL93 and LEU54) of MDM2 are obtained in all drug hits. Other important interface residues could be seen in binding were LEU57, ILE61, MET62, TYR67, GLN72, VAL75, PHE86, PHE91, HIS96, ILE99, TYE100, and ILE101. The binding affinity of Finasteride with MDM2 is in a similar range with Cobimetinib and Troglitazone; however, binding affinities of Loratadine and Drospirenone are comparatively lower. Finasteride/ Proscar is used for the treatment of alopecia and prostate cancer. In melanoma, the protective effect of finasteride on melanogenesis via downregulation of tyrosinase, TRP-1, MITF, and ACs expression has been demonstrated [60], further in vivo animal experiments are required for confirmation. Another study similar to this reported that Finasteride declines the risk of melanoma following prostate cancer [61]. The effects of Cobimetinib and Troglitazone were described above. The next drug hit is Loratadine which is routinely given to cancer patients in combination with other drugs and substantially improved survival in both breast cancer and cutaneous malignant melanoma [62]. In some instances, Loratadine has been linked to clinically apparent acute liver injury [63]. The use of Drospirenone in birth control pills and causing severe blood clots (thromboembolism) in women is conflicting [64]. Additionally, long-term use may increase the risk of breast cancer [65].

2.7. ADMET profile of top candidate repurposable drugs

The bioavailability radar plots in Figure 4(A) show that the hits are falling entirely within the physicochemical range on each axis and hence seem to fit into the bioavailability criteria. The pharmacokinetic profile in Figure 4(B) shows that all the top candidate drugs have high gastrointestinal (GI) adsorption. Except for Paliperidone and Troglitazone, all the top-hit drugs cross the blood-brain barrier (BBB) and can impair tumor development in brain metastasis. Three drug hits Cialis, Paliperidone, and Finasteride showed no violations for Lipinski, Ghose, Veber, Egan, and Muegge rules which suggests that these drugs are likely orally active drugs. Also, no alerts for PAINS and Brenk filters provide information that Cialis, Paliperidone, and Finasteride drugs didn’t contain potentially problematic fragments which are toxic or metabolically unstable. It is also essential to predict the interaction of drug hits with cytochromes P450 (CYP) and P-glycoprotein (P-gp) as these are key players in drug elimination through metabolic biotransformation [66]. Cialis, Paliperidone, and Finasteride all three are substrates of P-gp; and particularly, Cialis and Finasteride are not expected to cause unwanted adverse effects due to the lower clearance or accumulation of the drug metabolized by CYP450 isoforms (CYP1A2, CYP2C19, CYP2C9, CYP2D6, CYP3A4).

Overall, the results of molecular docking and ADMET profile suggested two potential strong-binding drug candidates (Cialis and Finasteride) of protein signatures (AKT1 and MDM2, respectively), shown to be orally bioavailable, non-toxic, and have good absorption and medicinal properties.

2.8. Molecular dynamics simulation (MDS) and docking validation

Since docking gives a static view of the binding interaction of compound hits into the active site of protein signatures, MDS gives a more clear idea about the physical movements of atoms and molecules with time by integration of Newton’s equation of motion [67]. Therefore, the two docked complexes (AKT1_ Cialis and MDM2_Finasteride), with the highest predictive binding affinity of -11.1 and -10.7 kcal/mol, respectively, along with safe therapeutic properties, were used for molecular dynamics study.

A 100 ns molecular dynamic simulation was performed to understand the molecular insights involved in the binding of Cialis in the active pocket of AKT1. The trajectory analysis (Figure 5A, (i) PL-RMSD) graph shows that the Ca atoms fluctuated in the range of 2.00–3.50 Å and stabilized after 80 ns of simulation with an RMSD value of 3.0 Å. Higher fluctuation (up to 3.50 Å) in RMSD was observed at 88 ns for the protein. The AKT1 backbone Cα experiences a gradual increase in deviation for the entirety of the simulation period. Also, Cialis in complex with AKT1 undergoes a gradual increase in deviation up to 40 ns and then decreases with minor deviations for a period of 40- 65 ns; and further maintained a deviation of ~2.20 Å for the remainder of the simulation period. The histogram (Figure 5A, (ii) P-%SSE) shows the secondary structure element (SSE) distribution for each trajectory frame over the course of the simulation. The percentages of the helix, strand, and overall secondary structure elements were found to be 25.06%, 14.63%, and 39.68% respectively. A total of 22 ligand contacts were formed with amino acids of protein (Figure 5A, (iii) PL-Contacts), from Leu156, Lys158, Gly159, Phe161 to Val164, Ala177, Lys179, Leu181, Met227, Tyr229, Ala230, Glu234, Glu278, Met281, Thr291, Asp292, Phe438, Asp439, Glu441, and Phe442. Figure 5A, (iv) Timeline graph represents that the ligand is stabilized by forming a majority of hydrophobic interactions with residues Leu156, Phe161, Val164, Ala177, Lys179, Leu181, Met227, met281, Phe438, and Phe442 with 1%–68% of simulation time. Hydrogen bonds are formed with residues Gly159, Lys179, Tyr229, and Asp292 throughout 1%–60% of the simulation trajectory. The RMSF value of the protein is coupled to the ligand in Figure 5A, (v) P-RMSF graph of simulated protein AKT1 for 100 ns, with X-axis measuring the average deviation of all protein residues (from 144-478) over time and Y-axis indicates RMSF (Å) values. The residues with higher peaks fluctuate the most during the simulation as determined by MD trajectories. Green-colored vertical bars indicate AKT1 residues that interact with Cialis and correspondingly, low RMSF values indicate the stability of the binding. Figure 5A, (vi) 2D-trajectory interaction diagram depicts that the hydrogen bond formed by the docking pose with Tyr229 is preserved in the MD trajectory pose with 59% of the total simulation time.

Further, again a 100 ns molecular dynamic simulation was performed to understand the molecular insights involved in the binding of Finasteride in the p53-MDM2 interaction surface. Figure 5B, (i) PL-RMSD graph shows that the Ca atoms fluctuated in the range of 1.00–2.25 Å and finally stabilized after 65 ns of simulation with an RMSD value of 1.5 Å. Higher fluctuation (up to 2.25 Å) in RMSD was observed at 15 ns. A slight divergence can be seen towards the end of the simulation around 92 ns. Since, the fluctuation lies under the permissible range of 1–3 Å, hence, can be considered as non-significant. The RMSD plot of Finasteride and MDM2 backbone were lying over each other. Hence, the formation of a stable complex can be inferred. Figure 5B, (ii) P-%SSE histogram shows the percentages of helix, strand, and the overall secondary structure elements to be 36.11%, 10.12%, and 46.24% respectively. Figure 5B, (iii) PL-Contacts histogram shows that a total of 20 ligand contacts were formed with amino acids of protein, from Gln18 to Ala21, Gln24, Leu54, Phe55, Ile61, Met62, Arg65, Tyr67, Asp68, Gln72 to Val75, Val93, His96, Ile99, and Tyr100. Figure 5B, (iv) Timeline representation suggested that the ligand is stabilized by forming Hydrogen bonds with Leu54, Gln72, His96, and Tyr100 with 1%-95% of simulation time. Hydrophobic interactions were formed with residues Ile19, Pro20, Leu54, Phe55, Ile61, Met62, Tyr67, Val75, Val93, and Ile99 over the course of 1%-36% simulation trajectory. Figure 5B, (v) P-RMSF graph of simulated protein MDM2 for 100 ns, with the X-axis measuring the average deviation of all protein residues (from 18-110) over time and Y-axis indicating fluctuation RMSF (Å) values. More fluctuations were seen at the C-terminal of the protein during simulation which indicates greater flexibility of the residues while other parts of the protein are more rigid and fluctuate less. Protein residues that interact with the Finasteride are marked with green colored vertical bars. Figure 5B, (vi) 2D-trajectory interaction diagram depicts that the hydrogen bonds formed by the docking pose with Tyr100 and Gln72 were preserved in the MD trajectory pose with 47% and 66% of total simulation time respectively. The RMSF value of the protein-coupled to the ligand.

3. 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 [68]. 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 in clinical trials [69,70]. 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 [71]; (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 integrated in 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 [72]. 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 a faster introduction to clinical trials and significantly reduce the costs for drug development.

Having in mind that metastasis is linked with the 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 in 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 important for network organization and molecular connections [73]. It is reported that cancer-associated proteins have large betweenness centrality as they control the communication between different components of a network [74].

Among non-topological properties, we have calculated the involvement of the motif constituents in the disease pathway, the gene prioritization score, and the average Log2 fold change for each motif based on the change in expression values of each node from non-invasive 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 neighbors 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 screening of FDA-approved drugs was employed and after binding affinity analysis top five candidate drugs selected for both the proteins AKT1 and MDM2, have evaluated for their safety profile. Following the candidate drugs with safe therapeutic properties and the least binding affinity with the signatures were subjected to MD simulations to check the conformational stability of the complex and dynamics/flexibility of the protein at 100 ns. The trajectory analysis confirms that the candidate drugs (Cialis and Finasteride) stabilized in the active pocket of protein signatures (AKT1 and MDM2) over the course of the simulation (at 80 ns and 65 ns respectively). Analysis of PL-contacts histogram and 2D-trajectory diagrams indicate that the hydrogen bonds formed by the residues TYR229 (in AKT1) and, TYR100 & GLN272 (in MDM2) are critically important residues and observed to play a predominant role in drug binding, thus contributing to the high stability of the complexes and could be further explored for in-vitro or in-vivo studies.

The PI3K-AKT signaling pathway is predominantly activated in melanomas, upon mutations of crucial melanoma drivers, in parallel with the MAPK (Ras-Raf-MEK-ERK) pathway. The BRAF/MEK inhibitors, which constitute the standard of care in the targeted therapeutics of melanoma, suppress the MAPK pathway. However, an active PI3K-AKT pathway offers alternative avenues for the evasion of MAPK-targeting regimens and disease progression. Indeed, members of the PI3K/AKT pathway are implicated in melanoma progression, metastasis, and acquired resistance to MAPK-targeting therapies. AKT1 activation frequently occurs through the silencing of PTEN, a tumor-suppressor gene encoding a phosphatase that acts on lipid and protein substrates. The major PTEN substrate is PIP3, which recruits AKT to the membrane and activates it by phosphorylation. By specifically dephosphorylating PIP3, PTEN suppresses the downstream signaling of AKT and, vice-versa its loss results in increased PIP3 levels and subsequent AKT1 activation. The combination of PTEN silencing with BRAFV600E expression in vivo leads to melanoma formation, while further ectopic expression of activated AKT1 in this genetic context co-operates with PTEN silencing to accelerate metastasis to the lungs and brain [75]. Progressing melanomas, in their vast majority, develop resistance to therapy due to the reactivation of MAPK signaling, commonly via alterations in BRAF, NRAS, and MEK1/2. However, a small proportion of resistant melanomas rely on the activation of the compensatory PI3K/AKT signaling cascade. PI3K/AKT-activating mutations in melanoma do not override proliferative arrest induced by BRAF/MEK inhibition but rather enable the survival of a small dormant subpopulation of MAPK-inhibited melanoma cells. The persistent exposure of this cell subpopulation to the constant selective pressure of BRAF/MEK inhibition eventually promotes the reorganization of signaling circuits and accelerates the evolution of tumor subclones that are highly resistant to targeted therapies [76]. Upregulation of the PI3K-AKT pathway is a critical event during the early and late evolution of resistance to MAPK pathway inhibition [77]. It has been proposed that AKTi combined with BRAFi-based therapy may benefit patients with tumors harbouring BRAF mutations along with PTEN deletions or AKT mutations [78]. In agreement with these studies, our analysis highlighted inhibition of AKT1 as an attractive strategy for preventing EMT-driven metastatic progression of melanomas with 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 the potential to overcome the suboptimal response rates for MAPK-targeting therapies [79]. The p53 gene represents a well-established player in carcinogenesis. Over 50% of all tumors carry p53 loss-of-function mutations, leading to the synthesis of functionally impaired protein products that are unable to transactivate genes that induce cell cycle arrest and apoptosis in response to oncogenic stress [79]. Tumors maintaining wild-type p53 exhibit other types of downstream p53 inactivation, such as hyperactivation of its endogenous post-repressor MDM2, a ubiquitin ligase that catalyzes p53 degradation and inactivation. Most melanomas retain a wild-type p53 and exhibit low p53 mutation rates, but instead show frequent inactivation of the CDKN2A, which eventually lead to MDM2 upregulation and subsequent p53 inhibition [80]. 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 wild-type tumors [81,82]. Shattuck-Brandt et al. [83] recently highlighted the therapeutic efficacy of MDM2 inhibition against p53WT melanomas with either a wild-type or a mutant BRAF background. They showed that an MDM2 antagonist (namely KRT-232) alone or in combination with BRAF and/or MEK inhibitors inhibits tumor growth in patient-derived xenografts (PDX) from 15 melanoma patients 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 [83]. On a similar note, our study revealed that a number of small-molecule MDM2 inhibitors show potential for preventing E2F1-driven metastatic progression.

Overall, our computational approach predicted AKT1 and MDM2 inhibitors as promising anti-melanoma drugs. If combined 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 the disease-free survival of melanoma patients. Future experiments are essential to confirm the metastasis-preventing potential of such combinations in melanoma animal models, in combination with standard-of-care melanoma regimens. Given that combinations of targeted therapeutics with immunotherapeutics could produce durable responses to clinical patients [84], another significant question that is worth further investigation is whether AKT1 and/or MDM2 inhibitors 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. 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 [85]. Similarly, 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 [86]. 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 [87] towards developing personalized therapies.

4. 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 the designing of next-generation targeted therapeutics. With the computational pipeline used in this study, we were successful in the identification of key protein signatures (AKT1 and MDM2) in melanoma from a core regulatory network that is based on a published E2F1 interaction map. In this work, a hybrid approach of logic-based modeling coupled with computer-aided drug design techniques was applied for the identification of drug candidates that can modulate the protein’s activity and could be possibly used for melanoma research. From virtual screening, the top candidate drugs based on the lowest binding affinity values against protein signatures were reported and evaluated for their safety profiles. MD simulations confirmed the stability of the two candidate drugs (Cialis and Finasteride) in complex with protein signatures over the course of 100 ns trajectory analysis. In conclusion, Cialis and Finasteride were predicted to be potent drugs to target AKT1 and MDM2 respectively; and may increase cell death in melanoma cancer cells, and this effect is mediated in the presence of E2F1. These findings would facilitate the development of effective inhibitors for clinical use in melanoma metastasis.

5. Methods

5.1. 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 the 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 [88]. More specifically, we calculated the average number of neighbors for each node in the network (degree) [89] and the density of connections among the neighbors of a node (betweenness centrality) [90] 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 [91].

5.2. Array data from aggressive melanoma cell lines

We used gene expression data from a previous study generated in SK-Mel103 and SK-Mel-147 cell lines (obtained from Dr. M. Soengas) with and without endogenous E2F1 depletion as described [3,85].

5.3. Motif prioritization

The regulatory motifs were prioritized using a ranking score for each motif considering key topological and non-topological properties with respect to the relevance of the melanoma phenotype. The motif ranking score is calculated using Eq. (1).

$${Rankingscore}_{ij}=\frac{W_{1j}}{2}\left(\frac{⟨ND⟩i}{\max ⟨ND⟩}+\frac{⟨BC⟩i}{\max ⟨BC⟩}\right)+W_{2j}\frac{⟨DP⟩i}{\max ⟨DP⟩}+W_{3j}\frac{⟨GP⟩i}{\max ⟨GP⟩}+W_{4j}\frac{⟨|FC|⟩i}{\max ⟨|FC|⟩}$$

(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 ⟨ND⟩ and betweenness centrality ⟨BC⟩) were divided into half for avoiding overemphasis on the topological properties and assigning equal weighting factors W_2j-W_4j to give equal importance to other properties (disease pathway association ⟨DP⟩, gene prioritization score ⟨GP⟩, Log2 fold change ⟨|FC|⟩) in motif prioritization. The equation generates a ranking score for each motif i (1…n) depending on the sets of values chosen for the weighting scenarios j (1 to13) shown in additional file 1a-c. Later, the top 10 motifs were selected from each of the weighting scenarios (13*10=130 motifs). Furthermore, a unique set of motifs were identified and processed for the construction of a melanoma-specific core regulatory network. The optimization of multi-objective function is discussed in detail [11].

5.4. Derivation of a core regulatory network

All the top-ranked motifs 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 are present in the E2F1 map. In addition, we added four EMT marker proteins (CDH1, VIM, ZEB1, and SNAI1) and direct connections with motif nodes (additional file 1d) in our regulatory core.

5.5. 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 [92]. 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 are 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 with other nodes in the regulatory layer to see the impact on the invasiveness (results in additional file 2a-e). 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.

5.6. Virtual screening of repurposable drugs

Virtual screening was performed as follows:

5.6.1. FDA-approved drug library preparation

The FDA-approved drug library was downloaded from the ZINC12 (http://zinc.dock-ing.org/zinc/) database. Since the library contains 2D structures, Open Babel 3.1.1 (https://pypi.org/project/openbabel/3.1.1/) was used to generate 3D energy-minimized structures to be utilized for docking studies.

5.6.2. Protein 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 [93]. Further, the coordinates of the active site residues were determined.

5.6.3. Binding affinity prediction using molecular docking

Virtual screening was carried out in PyRx v0.8 (AutoDock Vina-based) software [94]. The library compounds were first imported as SDF files in the open babel of PyRx and further energy minimization (using Universal Force Filed) 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 prepared FDA-approved drug library was subjected to docking against AKT1 and MDM2. To efficiently explore the docking conformational space, the search efficiency was set at 100%. For docking calculations, 9 con-formers were generated for each ligand-protein complex. The resulting ligand-docked poses were compared with the crystallographic poses based on ≤ 2.0 Å RMSD tolerance on the heavy atoms. The best predicted binding mode and the corresponding binding affinity (or binding free energy) in kcal/mol were selected. The more negative numerical values of binding free energy represent the better binding between a ligand and a protein signature. The docked complexes and graphical visualization were done in DS Visualizer [95].

5.6.4. Safety profile assessment of candidate drugs

The candidate drugs that bind to the protein signatures (in total, 1254 with AKT1 and 1257 with MDM2 respectively) were subjected to safety profile assessment based on ADMET risk, pharmacokinetics, drug-likeness, and medicinal chemistry friendliness prediction using the SwissADME tool (http://www.swissadme.ch/) [96]. The detailed properties/score values of the candidate drugs are provided (in the additional file 3a) along with the binding affinities (in kcal/mol) for both protein signatures separately (additional file 3b-c).

5.6.5. Molecular dynamics simulation (MDS)

Schrödinger LLC Desmond software was used to simulate the docked complexes for 100 nanoseconds [97]. For this, the complexes were pre-processed in Maestro Wizard for preparation and refinement steps. Missing side chains and loops were added using Prime, optimized with ProtAssign, and minimized using the OPLS_2005 force field. Further, the System Builder tool was implemented to build a system for simulation. The system was created in an orthorhombic box and the TIP3P water solvation model was used with default boundary conditions. The system was neutralized by adding counter ions and salt concentration set to 0.15M NaCl to mimic physiological conditions. The simulation parameters were defined as follows: simulation time 100 ns, recording interval trajectory 100 ps, ensemble class NPT, pressure bar 1.01325, and temperature 300 K. The stability of the simulation was verified by comparing the protein and ligand RMSD over time. The interactions between the protein and ligand were analyzed in the Desmond simulation interaction diagram tool.