Logic-based modeling and drug repurposing for the prediction of novel therapeutic targets and combination regimens against E2F1-driven melanoma progression

Melanoma presents increasing prevalence and poor outcomes. Progression to aggressive stages is characterized by overexpression of the transcription factor E2F1 and activation of downstream prometastatic gene regulatory networks (GRNs). Appropriate therapeutic manipulation of the E2F1-governed GRNs holds the potential to prevent metastasis however, these networks entail complex feedback and feedforward regulatory motifs among various regulatory layers, which make 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 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 identified that simultaneous perturbation of Protein kinase B (AKT1) and oncoprotein murine double minute 2 (MDM2) drastically reduces EMT in melanoma. Using the structures of the two protein signatures, virtual screening strategies were performed with the FDA-approved drug library. Furthermore, by combining drug repurposing and computer-aided drug design techniques, followed by molecular dynamics simulation analysis, we identified two potent drugs (Tadalafil and Finasteride) that can efficiently inhibit AKT1 and MDM2 proteins. We propose that these two drugs could be considered for the development of therapeutic strategies for the management of aggressive melanoma. Graphical abstract Supplementary Information The online version contains supplementary material available at 10.1186/s13065-023-01082-2.

The E2F1 network has 1015 nodes and 4174 edges.The average number of neighbors for each node in the network is 7.894 indicates that the network is well-connected (isolated nodes are 0).The average clustering coefficient of the network is 0.226 and a large network diameter of 8 shows the modular organization of nodes in the network [1,2].Further, a large average characteristic path length and a small average clustering coefficient suggest that the network has a small-world architecture and the small-world property reveals that signals can propagate very fast through the whole network [2,3,4].
The topological properties (i) node degree and (ii) betweenness centrality of the network were mapped for visualization.The node degree distribution (k) gives the number of neighbors for each node (n) i.e. the number of edges linked to a node in the network (Supp Fig 1).A power law (red line) y = ax b was fitted where y indicates the number of nodes that share a particular degree at x (a = 323.74and b = -1.274,R-squared = 0.842, correlation = 0.850).Genes with high degrees are in the right-down region of the graph.This leads us to the conclusion that the network has a scale-free topology, which is compatible with the network's sparse presence of a few high-degree nodes, called hubs, such as E2F1 (degree = 421).Such networks that contain a few hubs are generally heterogeneous in terms of node degree and are seen as being robust to single random perturbations [2,5].

Supp Fig 1: Node degree distribution of the network.
In Supp Fig 2, the betweenness centrality (Cb) measures the density of connections among the neighbors of a node (n) i.e. the amount of control that a particular node exerts over the interactions of the other nodes in the network [6].A power law (red line) y = ax b was fitted where y indicates the betweenness value for each node vs the number of neighbor nodes at x (a = 0.0 and b = 2.108, R-squared = 0.477, correlation = 0.953).Genes with high betweenness values are in the right=up region of the graph.From this, we conclude that E2F1, followed by E2F2 and E2F3, has very high node degrees and high betweenness values.This is understandable because the network was built by focusing on interactions around the E2F family.Other nodes with high node degrees are p53 [7] and MYC [8], are known for their role in melanoma tumorigenesis.Particularly, E2F1 has the greatest betweenness value of 0.423, indicating that it is crucial for the network's signal flow.

Robustness of E2F1 network:
We calculated the robustness of the E2F1 network using MORO, Cytoscape app [9] and compared it with the two existing large-scale signaling networks-the canonical cell signaling network (http://stke.sciencemag.org)and the human signal transduction network (http://www.bri.nrc.ca/wang).Directed networks are loaded for the analysis, the parameters are chosen (Update-Rule Scheme CONJ-DISJ and No. of initial states:10), and batch mode simulation is performed to calculate the robustness of the networks.Supp Table 2 displays the analysis results.E2F1 network's sRobustness (robustness against initial-state perturbation) and rRobustness (robustness against update-rule perturbation) values were found to be fairly close to the trend of the other two stable networks.It has been also suggested that a modular organization of cancer signaling networks is associated with patient survivability which suggests a relationship between modularity and network robustness [10].High modularity values (>0.5) indicate that the networks have dense connectivity between the nodes.

Regulatory network motifs:
The E2F1 network is of recurring structural patterns called network motifs which are used by a cell to process information and to govern dynamic response to external or internal fluctuations [11][12][13].Such network motifs include feedback and feedforward loops (FBLs and FFLs) that are present among various regulatory network layers and are commonly encountered in cancer networks [14,15].FBLs are characterized by the direct or indirect inhibition or activation of a node by its own target, and depending on the parity of the negative linkages in the loop, they can be positive or negative (Supp Fig 3) [16].However, FFLs are composed of an input node that controls an intermediate node, while the two also control an output node.Depending on the parity of the loop's negative linkages, the FFLs can either be coherent or incoherent [16].Both genes that code for proteins and those that do not are accompanied by such regulatory motifs [17].

Supp Fig 3: Examples of A(i) positive feedback loop and B(ii, iii) negative feedback loops from melanoma core network.
Such gene sets are highly interconnected, involved in specific phenotypes, and regulate one another through regulatory loops (motifs) from various pathways and the analysis of these can provide insights into the structure and dynamics of the network followed by the identification of therapeutic targets [18].These targets can be challenging to identify without the assistance of computational analysis due to the large and complex quantitative data, which requires computational algorithms for the analysis [19]; as well as the need to take into account a number of factors such as target-related safety issues, druggability, assayability, disease mechanisms, reduction in animal testing, time and cost considerations [20].Vera and coauthors, outlined the necessity of computational strategies in melanoma diagnostics and therapy [19].