Boolean dynamic modeling of cancer signaling networks: Prognosis, progression, and therapeutics

Cancer is a multifactorial disease. Aberrant functioning of the underlying complex signaling network that orchestrates cellular response to external or internal cues governs incidence, progression, and recurrence of cancer. Detailed understanding of cancer’s etiology can offer useful insights into arriving at novel therapeutic and disease management strategies. Such an understanding for most cancers is currently limited due to unavailability of a predictive large-scale, integrated signaling model accounting for all tumor orchestrating factors. We suggest that the potential of Boolean dynamic (BD) modeling approaches, though qualitative, can be harnessed for developing holistic models capturing multi-scale, multi-cellular signaling processes involved in cancer incidence and progression. We believe that constraining such an integrated BD model with variety of omics data at different scales from laboratory and clinical settings could offer deeper insights into causal mechanisms governing the disease leading to better prognosis. We review the recent literature employing different BD modeling strategies to model variety of cancer signaling programs leading to identification of cancer-specific prognostic markers such as SMAD proteins, which may also serve as early predictors of tumor cells hijacking the epithelial-mesenchymal plasticity program. In silico simulations of BD models of different cancer signaling networks combined with attractor landscape analysis and validated with experimental data predicted the nature of short-and long-term response of standard targeted therapeutic agents such as Nutlin-3, a small molecule inhibitor for p53-MDM2 interaction. BD simulations also offered a mechanistic view of emerging resistance to drugs such as Trastuzumab for

toward better decision-making in the clinical settings, and thereby help identify novel therapeutic strategies for improved cancer treatment at personalised levels.

K E Y W O R D S
Boolean dynamic modeling, cancer signaling, drug resistance, drug-response, prognostic markers, tumor progression INTRODUCTION Cancer is a multifactorial disease [1][2][3]. Cancer can be viewed as a result of disturbances of a delicate balance between cellular states, specifically proliferation and celldeath [4]. Identifying these disturbances and understanding their deleterious effects can help in better cancer prognosis, arresting tumor progression and improving the therapeutics needed for combating the disease. Systems biology-based Boolean dynamic (BD) modeling is an emerging approach for unraveling these disturbances which is paramount to developing strategies for restoring the delicate balance.
Worldwide, cancer incidence is 19.3 million per year, and cancer deaths stands at ∼9.9 million in 2020 [5]. Cancer incidence is expected to increase by ∼55% to 30.2 million per year by 2040 [6,7]. Current cancer care costs US$ ∼1. 16 trillion per year [8]. Large incidence rates coupled with rising associated costs makes the disease a public health priority and calls for arriving at effective strategies for preventing, managing, and combating it [7].
Cellular states are dynamically orchestrated by the underlying molecular network that responds to various external and internal stimuli [9]. Disturbances leading to dysregulation of a normal cellular state could be attributed to mutations affecting wiring of the network [2,9,10]. These mutations could be in proto-oncogenes and tumor supressing genes present in normal healthy cells. Influence of these disturbances could be far reaching, such as triggering metabolic reprogramming [10][11][12][13][14][15], evading apoptosis, a form of regulated cell-death [4,[16][17][18][19]. These disturbances are caused by oncogenic or non-oncogenic factors which alter the information flow through the network [20]. Oncogenic factors include overactive forms of certain protooncogenes whose gain-in-function mutation drives cancer initiation. Non-oncogenic factors refer to loss-of-function of tumor suppressor genes leading to cancer development [21]. Identifying sections of the molecular network that are dysfunctional, the normal processes that are disturbed and the extent of ensuing damage can offer novel insights to diagnose, combat, and prevent these diseases [22].
Therapies for treating cancer capitalize on the specific molecular aspects in cancer cells that delineate them from normal cells. Conventional therapies using first or second line of treatment take advantage of the genetic instability reflecting deficiencies in DNA repair in the cancer cells [23]. Alkylating agents such as chlorambucil, cisplatin are used as first line chemotherapeutic drugs for inducing DNA damage and thereby arresting cancer cell proliferation. A second line of treatment includes use of drugs such as methotrexate or of ionizing radiation. Traditional cytotoxic drugs and radiation therapy are typically weakly selective [21,23]. As a direct consequence, these therapies may affect other normal cells and thereby causing severe undesired side-effects [24][25][26]. This led to the advent of targeted therapy which involves use of drugs targeting a cancer-specific malfunction such as those caused by a mutation. Targeted therapies for various cancers offer improved remission and overall survival as compared to the conventional approaches [24][25][26]. However, they are often hampered by natural or acquired resistance to drugs. Knowledge of the mechanisms governing a cancer can offer insights toward addressing this challenge. While identification of precise causal mechanisms governing a particular cancer is often difficult, even partial knowledge of these can help finding effective therapy that may circumvent emergence of drug resistance [27]. A systematic assessment of a system-level model of oncological signaling which orchestrates cell fate can help discern such causal mechanisms. Further, such a model can serve as an in silico platform for rigorously testing systemlevel effects of a certain drug and also provide pointers for novel therapeutic approaches and other clinical decisions [28].
Systems biology driven cancer models permit integration of oncological signaling in a tumor cell and its microenvironment, and regulation at protein-protein and genetic levels. Among various systems-level modeling approaches [29], logical models, such as those based on Boolean networks, enable unraveling qualitative principles that elucidate mechanisms governing various behaviors elicited by cells [30][31][32][33][34] and particularly that of oncological outcomes [35]. Since its first use [36], BD models are routinely used in characterizing several naturally observed attributes of biological systems such as underlying complexity [37][38][39], self-organizing principles [40], redundancy, [41] and other nonlinear aspects [42,43]. BD approach offers a significant promise for understanding cancer's etiology from a signaling dynamics perspective and help identifying novel targeted therapies [35,44].
The present review is concerned primarily with BD modeling of cancer signaling networks with an emphasis on identifying principles governing various cancer-related processes. In the next section, we present a primer on BD modeling, which can be skipped by readers who are either familiar with it or are primarily interested in the BD modeling-based predictions. In the subsequent three sections, we review various BD modeling-based contributions leading to identification of diagnostic markers, deciphering mechanisms governing cancer progression, and understanding drug responses elicited by tumors.

BD MODELING OF BIOLOGICAL NETWORKS
Building a BD model of a normal or diseased cell necessitates availability of an underlying molecular network governing the signal transduction process. A biological network consists of entities such as proteins, genes, small molecules, and the causal relationships as interactions between them. In this section, we first present a general formalism for BD modeling. Subsequently, we describe systematic approaches for constructing cancer-specific annotated signaling network.

A primer on BD modeling
In this primer, we describe the BD modeling approach, its implementation and associated strategies using a cancer progression signaling network as a motivating example. These are applicable to signaling networks such as proteinprotein networks (PPNs) [45], gene-transcriptional networks (GTNs) [46], multi-cellular networks, [47,48] and a combination thereof [49,50]. PPNs consisting only of interacting proteins are typically employed in cancer modeling when the objective is to understand how signal flow from a receptor is disturbed due to a malfunctioning entity in the network. Moreover, it is considered for identifying potential diagnostic and therapeutic protein markers. On the other hand, GTNs typically consist of transcription factors as entities. An interaction in a GTN represents regulation of its synthesis via the transcription and translation machinery by other transcription factors. GTNs are considered in cancer modeling for a variety of purposes such as understanding reprogramming triggered by oncogenic genetic abnormalities [51,52]. Multi-1 , 2 , … , 6 represent the Boolean values corresponding to the six nodes in (A) cellular networks consist of multiple cells themselves as nodes or as compartments with sub-networks in each of them.
For demonstrating BD modeling of biological networks, we consider a six-node network consisting of signed, directed interactions based on TGFβ signaling in a cancer cell [53] as a motivating example ( Figure 1A). TGFβ signaling is implicated in cancer progression [54]. TGFβ, a secreted cytokine protein, upon binding to its receptor triggers signaling via MAPK cascade that culminates in activating the synthesis of Snail1 protein. Snail1 transcription factor activates Zeb1 and represses miRNA200. Moreover, Snail1 represses the transcription of miRNA34a and in turn miRNA34a inhibits Snail1 protein by modulating its translation. Thus, Snail1 and miRNA34a are locked in a double negative feedback loop. A similar double negative feedback loop exists between Zeb1 and miRNA200. On the other hand, Ovol2 and Zeb1 interlock each other via a transcriptional double negative feedback loop. Further, Zeb1 and Ovol2 transcription factors repress miRNA34a and the synthesis of TGFβ, respectively. For the sake of brevity, we represent TGFβ, miRNA200, Snail1, Zeb1, Ovol2, and miRNA34a, respectively as n 1 , n 2 , . . . , n 6 .
BD network model : A node (n i ) in the network is assumed a Boolean variable taking a logical value FALSE or TRUE captured by { = 0 1}, respectively, representing inactive, that is, not-expressed (OFF) or active, that is, expressed (ON) forms of the entity. In case of proteins, expressed and not-expressed reflect high and low concentrations, respectively. Node n i is associated with a Boolean function f i consisting of logical operations AND, OR, or NOT governed m-logical-rules corresponding to m binary inputs to it. f i determines the TA B L E 1 States reached in the first few timesteps starting from initial state ( s 0 = 101001) using synchronous (S), deterministic asynchronous (DA), general asynchronous (GA), and random order asynchronous (ROA) update schemes. Nodes updated at a time t using the Boolean functions in Figure 1B are in red with the value of all others retained as is. While FP reached for S, and DA methods are specified by the state at t = ∞, that for GA and ROA methods depend respectively on randomly chosen update order and the randomly chosen permutations at every . For the case of DA, timescale is specified in Figure 1B. Permutations used for ROA for = 1, 2, 3, respectively, are 431562, 564132, and 164325  Figure 1B). For example, n 2 and n 5 inhibit n 1 resulting in 1 ( 2 , 5 ) = ( 2 ) ( 5 ), where an OR logic is assumed for the two inhibitory inputs (Figure 1). The dynamics to the Boolean model is introduced by updating using evaluated at the prior values of those nodes contributing to the logical-rules in it [36]. The value taken by depends on the update scheme adopted, which we discuss next.
A synchronous (S) update scheme is one where ∀ is updated simultaneously, that is, at time + 1, +1 = ( 1 , 2 , .., 6 ) [46,[55][56][57]. Starting from a random initial state 0 = { 0 1 , 0 2 , 0 3 , 0 4 , 0 5 , 0 6 } = {1, 0, 1, 0, 0, 1}, one S update results in network's state ( 1 ) being 110111, with each digit representing the Boolean value (Table 1) taken by the corresponding node [43,58]. For the sake of brevity, henceforth, a network's state will be presented as six binary digits randomly assigned to six nodes in Figure 1A. Upon repeated updating, after visiting several transient states, the network reaches a Boolean steady state of 010011 (Table 1). Such a steady state is referred to as an absorbing-state or fixed-point (FP) attractor or FP. Set of initial states culminating in FP forms its basin of attractor. Note that value taken by nodes at FP indicates the active or inactive entities in the cell under homeostasis and thereby reflecting its phenotype. Besides FPs, a network can also reach cyclic attractors, in which transients hover around periodic sequence of states [59]. S update method is known to lead to several spurious cyclic attractors [59].
In biological systems, all nodes are not necessarily activated simultaneously. Asynchronous update scheme via deterministic asynchronous (DA), general asynchronous (GA) and random-order asynchronous (ROA) methods permits staggered updating of nodes [43,60]. Nodes to be updated are either randomly chosen or specified by a priori known discrete activation timescale for every entity [43,60]. For DA method, at every , only nodes with ≤ are updated [48]. Table 1 shows 1 , 2 , 3 , and ∞ (FP) obtained by starting from 0 . FP reached from an initial state is identical and unique for both S and DA methods. On the other hand, GA and ROA methods are applicable when is unavailable and further, permit capturing stochasticity, a phenomenon naturally present in biological systems [61,62]. In GA method, at every , a randomly chosen node is updated leading to reaching an update-order dependent FP from an initial state [45]. Same 0 leads to different dynamic states for S, DA, and GA methods (Table 1). In ROA method, at every , 1. updating is further resolved into as many intervals as that of the number of nodes, and 2. all nodes are updated in a certain sequence (permutations), chosen uniformly randomly [50]. GA approach captures the stochasticity due to the cascading effect of updating randomly chosen individual nodes over several timesteps. On the other hand, ROA allows incorporating the stochastic behavior originating from the interdependency of the nodes at all timesteps. Table 1 shows states reached in three successive instances, each employing different permutations, with 0 = 101001. All three update methods can lead to cyclic attractors. In fact, the spurious cyclic attractors predicted by the S method can be circumvented by using asynchronous update schemes [59]. However, GA and ROA can in addition lead to complex attractors wherein the transients meander over a set of states in no specific order.
An extension of a BD model is the multivalued discrete dynamic network model wherein can take multiple integers [31]. This allows capturing varying levels of each of the nodes. Further, every rule in the function governing the dynamics of a node is scaled with a weight factor that F I G U R E 2 (A) Complete state transition graph (STG) showing states (dots) and one-step transitions (arrows) between them obtained using synchronous update method. Self-arrow to a state indicates FP. Red dots show the transients obtained to reach the FP 010011 (green) with initial state s 0 being 101001. Orange dot represents FP 101100. (B) Complete STG obtained using ROA update method. (C) Partial STG, a section of the complete STG in (B), consisting of seven states (dots highlighted in [B]) and all one-step transitions between them. Initial state is 0 = 101001. Probability of reaching FPs from 0 is shown within the self-arrow. Numbers (in red) placed next to the one-step transition arrow captures the associated transition probability calculated by incorporating the corresponding number of permutations in Table 2 into Equation (1) captures the multivalued nature of . Simulations of the discrete dynamics of the model using the multivalued and the modified functions are performed using synchronous and asynchronous update schemes described above. Besides the ensuing computational complexities, deciding the range for the multivalued variables is often a challenge. A comprehensive list of toolboxes for performing discrete dynamic model simulations is available at The CoLoMoTo consortium [63].
Long-term behaviour of a BD model : Analyzing the transients, due to one-step transitions, reached while traversing from 0 to desired FP can offer insights into network's inclination to reach a phenotype. A collection of such one-step transitions achieved by starting from all possible 64 ( = 2 6 ) initial states is the six-node network's statetransition-graph (STG) which has states as nodes connected by one-step transitions as directed interactions and can be constructed for all four update methods [59]. For the case of synchronous update, Figure 2A shows the complete STG. The states reached while traversing from 0 = 101001 to reach the FP 010011 are captured in Figure 2A (red dots). STG also shows that the FP 101100 (orange dot) can be reached via synchronous update only when the initial state is itself. Thus, all states other than FP 101100 belong to the basin of attraction of FP 010011. Note that the path taken to reach an FP from an initial state is unique when S update method is employed, as corroborated by the STG for the six-node network ( Figure 2A). On the contrary, the path traversed to reach an FP starting from an initial state is non-unique for the case of GA and ROA methods. Thus, stochasticity incorporated in GA and ROA methods permits expansion of the basin of attractors for reaching different FPs, a feature exhibited by real networks [45]. In Figure 2B, we show the complete STG for the six-node network when ROA method is employed. STG captures the permutation-dependent path to reach an attractor from any state that the network can take. For example, starting from 0 = 101001, one-step transitions caused by permutations 451236 and 416325 lead to states 100101 and FP 101100, respectively. Thus, depending upon the path dictated by the permutations chosen for the intermediate timesteps, ROA method takes 0 = 101001 to either 010011 or 101100 FP, as shown in the partial STG in Figure 2C. Thus, multiple paths could lead to same FP. Note that multiple permutations could cause same one-step transition. For ROA method, since every one-step transition will have an underlying permutation, STG will have utmost 46, 080 ( = 64 × 6! = 64 × 720) transitions. For a large network, finding complete STG is computationally tedious. However, experience suggests that a much smaller STG could aptly capture the required characteristics, particularly that of the driver entities that govern the network's dynamics and of the attractors.
Ability of a network to settle into an FP is quantified by the probability which can be estimated from the STG. This probability can help unravel network's features such as driver entities. The first step in this quantification process is finding the 1-step state transition matrix . An element of this matrix specifies the probability of achieving a one-step transition from to and is given by [59,64] = ∑ where, is the number of permutations causing the transition. Number of permutations causing the one-step transitions in the partial STG is in Table 2, and the TA B L E 2 Number of permutations causing the one-step transitions in the partial STG in Figure 2C  corresponding transition probabilities estimated using Equation (1) are presented in Figure 2C. An iterative solution of gives the probability of reaching different FPs starting from 0 . Setting the initial vector where is the number of states in the STG, gives the average or steady-state probability of absorbing into different FPs of the STG. Steady-state probabilities of reaching the two attractors 010011 and 101100 estimated using Equation (2) are 0.9 and 0.1, respectively ( Figure 2C).

Biological network construction
Often network of interest may be unavailable. As a result, constructing a reliable, domain-specific signaling wiring diagram of a normal healthy cell is the logical starting point for building a BD model. Note that domain-specific refers to specific tissue of interest and the type of the wiring diagram considered. The first step toward constructing a network is curation of information about the nodes and interactions between them. A comprehensive review of nature of different entities and interactions is in Papin et al. [65]. Two broad approaches for biological network curation and construction are data-driven objective and knowledge-driven objective [66,67]. Data-driven objective employs available experimental data, for example, in-house generated microarray/proteomics data, those in GEO repository [68], ExpressionAtlas [69], GeNet [70,71].
Entities and interactions thus distilled can be visualized as a wiring diagram using software such as CellDesigner [83], Cytoscape [84], which stores it in a portable, shareable, machine-readable format [85]. Collected information is further structured and refined to include cancer-typespecific annotations such as tissue, nature of activation, associated logic. Annotations form the basis for arriving at a logical function governing activation of a node in BD modeling framework, detailed in section 2.1. After construction, a series of network reduction by identifying key structures [86] or removing certain variables such as frozen nodes [39,43,[87][88][89] curtails computations. Systematic step-wise network construction methodology and compendium of resources are available in Türei et al. [67] 3

MARKERS GOVERNING PHENOTYPE SWITCHING
Early diagnosis of tumor formation is extremely useful in employing effective treatment regimens [90]. Cancer being a multifactorial dynamic signaling problem involving multitude of interacting entities, aberrant nodes serves as markers for (early) diagnosis, thereby facilitating better disease management. For cancers exhibiting strong variability in clinical response to therapies, markers distinguishing disease severity and reflecting therapeutic outcome help stratifying patients for an appropriate treatment [90][91][92][93]. In this section, we review BD approaches to model switching of inflammation or apoptosis to proliferation phenotype and thereby identify diagnostic markers. Networks considered in this section span from only-apoptosis TA B L E 3 Cancer diagnostic markers predicted by Boolean dynamic modeling of various networks. Construction objective adopted and the network statistics are specified. S and GA refer respectively to synchronous and general asynchronous update methods described in section 2. to pan-cancer to a specific cancer-tissue. These networks could consist of entities within a cell or additionally those accounting for multicellular signaling events in the tumor microenvironment. For every case considered, network details, the Boolean modeling update scheme adopted and key findings are presented in Table 3. Malfunctioning and overexpression of an entity due to a certain mutation is captured by node deletion or node over-expression (ON), respectively.

Inflammation to proliferation phenotype switching
Cancer cells secrete inflammatory cytokines and chemokines and thereby induce inflammation in the neighbouring cells [2,94]. Inflammatory cells in turn promote cancer progression by making them to commit to proliferation phenotype [95,96]. This phenotype switching is an important step in cancer progression. Entities involved in this switching process can thus serve as diagnostic markers [96].
Analysis of FP, cyclic, and complex attractors of a colitisassociated colon cancer network showed that transient dendritic cell activation, an essential step for a successful immune response, is needed for a release of inflammatory cytokines and chemokines [49]. Entities considered in the network include immune cells in tumor microenvironment such as dendritic cells, secreted cytokines and chemokines, and intracellular signaling nodes. This analysis further showed that simultaneous sustained dendritic cell activation and p53 protein inactivation may govern the inflammatory to proliferation phenotype switching [49]. BD attractor analysis of an (pan-cancer) apoptosis network corroborated that inactivation of p53 results in breaking key feedback loops which serve as an alternative for irreversible apoptotic cell fate when those with normal caspase-3 are dysfunctional [97].
Recently, a novel approach of finding the Hierarchical Partitioning for the Phenotype, implemented on the colitisassociated colon cancer network [49], offers promise in finding the global attractors corresponding to a specific phenotype [98]. This approach involved pre-evaluating the Boolean functions with the phenotype-specific value corresponding to the external nodes. This results in a simplified sub-network with interactions governed by fullyor semi-updated rules that uniquely specify the attractors corresponding to the phenotype. In order to identify the markers that may be controlling the proliferation phenotype, the model was analysed by assuming sustained activation of adenomatous polyposis coli protein representing premalignant epithelial cells and dendritic cells.
This analysis suggested that SOCS, JAK, and STAT3 may be controlling the proliferation global attractor. Moreover, the analysis also revealed that during sustained dendritic cell activation, the model does not predict apoptotic phenotype.

Proliferation phenotype by evading apoptosis
One of the hallmarks of cancer is evading apoptotic signaling [2]. In normal cells, TNF superfamily cytokines TNFα and FASL [99] maintain a balance between the apoptotic, necrotic, and proliferative phenotypes [100]. Note that the precise role of TNF superfamily cytokines in driving different cell-types into these different phenotypes is still under investigation [101]. Under various pathological conditions, cells secrete large quantities of TNFα leading to an elevated leads of the cytokine in the tissue microenvironment. Steady-state probability of a TNFα-FASL signaling network for reaching different phenotypes could delineate the roles played by various cancer-associated proand anti-apoptotic genes [102]. Even when TNFα is active, deletion of APAF1, BAX, FADD, or Caspase 8, or overexpression of BCL2 or NFκB is sufficient to switch the steady-state probability significantly toward anti-apoptotic phenotypes [48,102]. FP attractor analysis of a BD model of neuronal cell network activated simultaneously with FASL and neuronal growth factor showed that overexpressed HSP70 shock protein dominates over HSP27, HSP40, and HSP90 in promoting proliferation of cells, which would otherwise commit to apoptosis [103]. Systematic experimentation on Neuro2a cell line substantiated this prediction. Further, the attractor analysis identified BCL2, IAP, cFLIP, and NFκB as key players of pro-survival cell fate.
An integrated PPN of signaling pathways was analyzed systematically for assessing the nature of stability of attractors corresponding to apoptosis, proliferation, and quiescent phenotypes [104]. Analysis of FP and complex attractors of the network under 32 distinct environmental conditions predicts that 1. under normoxia conditions, deletion or overexpression of 10 proteins such as EGFR, NFκB, Ras can transform quiescent cells into proliferative ones, and 2. under hypoxia conditions, the repertoire of deleted or overexpressed proteins promoting proliferation increases by 7 over and above those under normoxia conditions. Further, deletion of APC, p53, Smad4, pTEN and overexpression of RAS, Tcf, Akt in the network to mimic colorectal carcinoma [105] led to switching of 97% apoptosis/quiescence attractors under normal conditions to 50% taking proliferation phenotype. BD model of an EGFR+ breast cancer network predicted that proliferation phenotype requires malfunctioning of PTEN or GSK3 β [106]. Moreover, structural and topological analysis of the attractors of a T-cell large granular lymphocyte leukaemia network led to identification of entities such as RAS, PLCG1, IAP, SMAD as key for apoptosis to proliferation phenotype switching [45] with the role of SMAD substantiated experimentally.
MAPK pathways are known to be sentinels of cancer. TGFβ signaling exhibits dual role in cancer cells by context-specifically regulating apoptosis and proliferation phenotypes. While it acts as a tumor-supressor via p38-MAPK pathways in normal or early cancer cells by arresting cell-cycle or triggering apoptosis, its tumor-promoting role in the late-stages leading to cells choosing proliferation attractor is attributed to the dysfunction of SMAD2 and SMAD3 proteins [107]. Ability of signaling through activated MAPK pathways to facilitate a balance [108] between the proliferative and non-proliferative phenotypes makes the FGFR3-mutated bladder carcinomas less aggressive [109].
In summary, BD modeling approaches implemented on different cancer-tissue-specific networks led to identification of markers that drive inflammation to proliferation phenotype switching and lead to acquiring proliferation cell-fate by evading apoptosis. Tracking these markers can help early stratification of patients and thereby facilitating improved management of the disease.

CANCER PROGRESSION: EPITHELIAL-MESENCHYMAL PLASTICITY
Cancer progression involves tumor growth, evolution, invasion, and spreading. In order to achieve these, cancer cells capitalize on the normal cell-biological program called epithelial-mesenchymal plasticity (EMP) [110,111]. EMP facilitates adaptation of cancer cells to a new environment [112]. Moreover, it may even induce drug resistance in them [113]. EMP consists of multiple possible phenotypic states with complete epithelial-to-mesenchymal transition (EMT) and complete mesenchymal-to-epithelial transition (MET) at the two ends of the spectrum [53,[114][115][116][117]. A cancer cell hijacks the EMT process to leave the primary tumor site for invading surrounding tissues and thereafter, migrate, and enable distant metastasis [112,118]. Significant fraction of cancer-related mortality is associated with metastasis [119]. An outstanding question is, what are the factors that dictate distal-organ invading abilities of a cancer cell? Mechanistic insights into EMT can offer clues to arresting or controlling this cell-biological program hijacked by the cancer cell. An emerging strategy to unravel topological features and causal mechanisms TA B L E 4 Cancer progression markers predicted by Boolean dynamic modeling of various networks. Construction objective adopted and the network statistics are specified. S, GA, and ROA refer respectively to synchronous, general asynchronous, and random-order asynchronous update methods described in section 2. governing tumor invasion, EMT, and migration is to systematically analyze GTN constructed using large geneexpression data from different cancer patient cohorts. Table 4 contains the network details, the BD modeling update scheme adopted, and key findings for the cases reviewed in this section. BD model of bladder and breast cancer juxtaposed with independent experimental studies led to prediction of E2F1, TGFBR2, and EGFR being important molecular signatures for regulation of tumor invasion and thereby EMT [86]. TGFβ induces a trans-differentiation program enabling the epithelial phenotype characterized by adhesion properties such as tight junctions to create mesenchymal derivatives such as better mobility due to loss of E-cadherin and overexpression of vimentin [54]. This induction of EMT occurs in cooperation with the activation of SONIC Hedgehog and Wnt-SHH pathways [120]. Activation of Wnt signaling, implicated in EMT triggered by different cancer [121,122], may even induce drug-resistance [122,123]. A systematic combinatorial node perturbations on a BD model of EMT in Hepatocellular carcinoma complemented with siRNA-based experimentation showed that inhibition of SMAD along with RAS, NOTCH, or SOS/GRB2, which hamper the key feedback loops, is effective in suppressing TGFβ-induced EMT [124]. However, inhibition of SMAD alone results in a hybrid EMT state where cells exhibit both epithelial and mesenchymal traits.
Local invasion by and eventual late-stage migration of cancer cells must have myriad signatures in the early stages of metastasis. In the case of gut cancer, simultaneous Notch overexpression and p53 deletion in a TGFβ-induced EMT network drives settling of cells exhibiting invasion or EMT or other phenotypic attractors into a metastatic one [125]. This suggests that the synergistic effect could have favored migration, as has been corroborated by experiments in mouse gut [126].
Phenotypic variability expressed by small-cell lung cancer, which is aggressive in terms of its ability to relapse post-first line of treatment, constitutes presence of both neuroendocrine/epithelial (NE) and nonneuroendocrine/mesenchymal (ML) behaviours [127]. Attractor analysis of the underlying transcription network, constructed using expression data from clinical samples [82,128], revealed specific coexistence of NE and ML sub-populations. The coexistence of such epithelial and mesenchymal characteristics has been validated in vitro on several small-cell lung cancer cell lines [129]. Further, the Boolean value taken by the identified attractors suggested that while INSM1, POU3F2, SOX2, SOX11,FOXA2, OVOL2 genes may be responsible for the NE behavior, MYC, NFKB1, SMAD3 genes are active in the ML phenotype.
In vitro single-cell study further revealed that after cells were treated with first line (cytotoxic) drugs, the two sub-populations transited into hybrid EMT phenotype [129].
Cells are constantly exposed to internal and external noise, sources for which are plenty [61,62]. Even though BD models predict the attractors corresponding to different cellular states, it is well-known that EMT may involve multiple intermediates [130][131][132][133]. Given the inherent stochastic environment cancer cells may experience, gaining confidence on the causal mechanism governing these different attractors requires stability analysis of these multiple states [134]. Relative stability of an attractor is quantified by comparing effort required for the network to switch from one FP to another [112]. Recently, dynamic stability analysis of attractors of a generic GTN underlying EMT showed that attractor corresponding to the epithelial phenotype is more stable than the others [134]. The dynamic stability was estimated based on the one-degree neighborhood of cell states [134].
New emerging evidences suggest that blocking EMP in both directions, that is, EMT and MET, and holding the cancer cells in its current relative epithelial-mesenchymal state may help reduce its ability to adapt to fluctuating environments and thereby, decrease its metastatic potential [135]. Recently, a combined topological and BD analysis of the EMT network identified key positive feedback loops that may govern E/M phenotypic multistability, the extent of its influence, and the deletion of interactions controlling plasticity [136]. In order to account for heterogeneity, dynamic variability in the Boolean value taken by an entity was introduced by linking the decision for updating a node at a certain timestep to the relative number of instantaneous activating or inhibitory inputs at that instant [137].
In summary, various BD modeling approaches have been employed to model cancer progression via the EMT for different cancer types. Molecular signatures governing EMT were identified using systematic perturbation analysis. Dynamic stability analysis of the attractors revealed key underlying topological features modulating the phenotypic multistability involved in the EMT. These findings can be utilized to better control the cancer cell's ability to hijack the EMP.

THERAPEUTIC OUTLOOK: DRUG RESPONSE AND RESISTANCE
The plethora of drugs is available for different cancers [138], several of which are in various stages of clinical trials [139]. Various first and second line of therapeutic agentsfor example, Gefitinib for the non-small cell lung-cancer [140] approved for targeted therapies are typically based on the mutation that drives tumor formation and do offer some promise like improved remission compared to traditional chemotherapy [24,25,26,141]. However, the overall efficacies of these approaches are also fraught with varying susceptibility or side-effects or intrinsic resistance to drugs or induced drug-resistance, just as in the case of the conventional therapies. Thus, finding cancer-specific targeted therapies that circumvent these challenges continues to remain an open quest. BD modeling-based strategies for attractor switching, particularly to apoptosis, identified in silico are valuable in finding potential candidates for design of better interventional strategies that can be tested in vitro, in vivo, and in clinical settings [44]. Details of the network and the BD modeling scheme adopted and key findings for the cases reviewed in this section are summarized in Table 5.
Use of Trastuzumab or Herceptin, Pertuzumab and small-molecular inhibitor erlotinib drugs are well-known line of treatment for HER2-positive breast cancer specifically to block Erb2/Erb1 receptor induced deregulation of downstream signaling leading to tumor progression. BD modeling of a cancer cell line-specific signaling network constructed by constraining time-course proteomic data from three HER2-positive breast cancer cell-lines exposed to these drugs led to identification of specific new interactions that may strongly influence MAPK and PI3K activation patterns [142]. Specifically, in HCC1974, which is Trastuzumab resistant, the identified two new edges PDK1 → ErbB-2 and p70S6K→Akt introduced feedback loops resulting in suppression of apoptosis phenotype via PI3K pathway and thereby causing stabilization of Trastuzumab drug-resistance by promoting oncogenic effects. Observations in the clinical breast cancer samples substantiate these predictions [143]. Use of inhibitors such as celecoxib [144] against PDK1 could help overcome such a drugresistance. Herceptin-based treatment targets the extracellular domain IV of HER2 which in turn down-regulates activation of MAPK pathways, viz., specifically ERK1/2, p38, and JNK1/2, and thereby regulating signal flow toward mitogenic and survival phenotypes [145]. However, clinical data revealed that ∼70% of the initial responders to Herceptin subsequently experienced progression to metastasis [146] suggesting the possibility of drug-resistance [147,148]. BD modeling of the underlying network showed that dynamic variations in the levels of several dual-specificity phosphatases (DUSPs), which too regulate MAPK pathways, could provide insights into targeting certain DUSP to overcome resistance against Herceptin drug treatment [149]. In fact, gene expression analysis on clinical samples from a HER2-positive cohort receiving Herceptin treatment shows that higher expression of DUSP4 may be correlated with poor patient survival [150] making it a potential target. Therefore, inhibiting DUSP4 expression could help overcome Herceptin drug-resistance. Probabilistic Boolean modeling of a breast cancer signaling network corroborated with siRNA gene silencing experiments on MDA-MB-435 cell line suggested that Mcl1 is a good drug target for influencing cancer cell growth [151].
Systematic attractor landscape analysis of 45 distinct cancer cell-specific p53 networks perturbed with smallmolecular inhibitors predicted the phenotype switching affected [152]. Small-molecule inhibitors used include Nutlin-3 [153] for p53-MDM2 interaction, GSK2830371 [154] for Wip1, MK-2206 [155] for Akt, CDK2 inhibitors [156] for Cyclin E, and Navitoclax [157,158] for BCL2 family proteins along with or without DNA damaging drug TA B L E 5 Novel therapeutic targets predicted by Boolean dynamic modeling of various networks based on cancer drug response and resistance. Construction objective adopted and the network statistics are specified. While S and GA refer respectively to synchronous and general asynchronous update methods described in section 2.1, PB refers to probabilistic Boolean modeling [166,167]  etoposide. Networks used were constructed using genomics data from 83 human cancer cell-lines. Specifically, this study showed that single inhibition of Akt and combinations such as inhibition of p53-MDM2 and Wip1 with or without DNA-damaging drug may offer high efficacy. While Akt is a promising therapeutic target for overcoming cell-type-specific drug-resistance, given that it has multiple functional sites and is involved in several feedbacks [159], triple inhibition instead of single may be necessary [160]. A therapeutic drug inhibiting a target node can cause cells to offer adaptive resistance, wherein the inhibitor, such as that against RAS, may induce dynamic reprogramming of the signaling process to adapt to new treatment conditions [161][162][163]. In such cases, it is imperative to identify the aspects of the signaling network responsible for such dynamic reprogramming and arrest the same. For instance, employing novel probabilistic BD modeling on a colorectal cancer network showed that blocking a feedback regulation involving Src may overcome the adaptive resistance during targeted therapy involving inhibition of MAPK pathway [164,165]. Probabilistic BD modeling permits including variability in selecting the Boolean logic employed in the Boolean function of a node based on the weights of the interactions input it [166,167].
Emerging evidences suggest that interaction of cancer tissues with microenvironment consisting of inflammation-related cells such as macrophages may influence drug responses [168]. In vitro studies on MDA-MB-231 cell lines suggest that the efficacy of etanercept, a TNFα inhibitor may be compromised due to activation of other pathways triggered by macrophages-secreted cytokines and chemokines [169]. BD model of the TNFα signaling network integrated with JAK-STAT and PI3K-AKT pathways triggered by interleukins and other growth factors suggests that cells exposed to etanercept may continue to express survival phenotype due to NFκB activation triggered by these extracellular stimuli [170].
Systematic in silico analysis of perturbed BD model of colitis-associated colon cancer along with experimental validations revealed that concurrent activation of ceramide and inhibition of PI3K/AKT pathway could lead to an effective anti-cancer response in tumor cells [49]. Availability of several inhibitors [171] of PI3K/AKT pathway that are under clinical evaluation and of C2-ceramide that increases endogenous levels of ceramide [172] makes this combinatorial targeting attractive. Attractor analysis of a BD model of gastric adenocarcinoma cells suggested that besides the traditional combinatorial inhibitors such as that for PI3K and MEK, and for MEK and AKT currently investigated in clinical trials for different cancers, simultaneous targeting of TAK1 and PI3K, and of TAK1 and AKT could favor apoptotic phenotype [173].
Cancer tissue-specific BD models led to identification of the causes of poor efficacy of a treatment regime and the mechanism such as dynamic reprogramming underlying acquired drug-resistance during a certain therapy. BD analysis further led to identification of potential combinatorial treatment targeting specific nodes or interactions that could circumvent emergence of drugresistance.

PERSPECTIVES
Ability to predict and understand causal mechanisms governing cancer incidence and progression is paramount to preventing, managing, and combating the disease. Cancer is widely recognized as a result of disturbance in the signaling problem due to aberrant cues, source for which could be external or internal to cells [2]. BD modeling approaches are used extensively to predict prognostic markers, to decipher mechanisms governing cancer progression and to identify novel therapeutic targets for improving disease management strategies. While these predictions offer novel insights, most of the BD models in literature (Tables 3-5) primarily consider intracellular signal flow, which accounts for only a small fraction of the tissue-level cancer signaling process. It is well recognized that tumor formation and cancer progression are orchestrated by complex dynamic communication between various inter-and intra-cellular machinery [2,174,175]. We present in Figure 3 a block diagram capturing dynamic communication that may be present between different cancer-related signaling machinery. The functioning of each of these machineries is in turn governed by interacting entities such as protein, mRNA, small molecules at multiple levels.
Unraveling etiology of tumor formation and its eventual dynamic responses require a cancer-specific model that incorporates the multi-cellular signaling process elucidated in Figure 3. Such a model will juxtapose signal flow within a cancerous cell and that in other cells in the tumor microenvironment, with detailed communication between them. This poses challenge at two-levels, viz., constructing a reliable BD model of the comprehensive network and performing massive computations of the same.
A step towards addressing the first challenge could be constructing and analyzing cancer-specific BD model of each of the signaling machinery (in Figure 3) under isolated conditions. Subsequently validate the predictions using experimentally measured threshold levels of observable nodes and exhibited phenotypes. Key bottleneck here could be lack of knowledge of 1. entities and interactions participating in the considered machinery and 2. which sub-machineries may further be involved. Should a network under normal conditions be available, the first bottleneck could be circumvented by identifying the perturbations in the nodes or interactions or associated F I G U R E 3 Block diagram of different signaling machinery (boxes) influencing cancer incidence, orchestration and progression along with a few important entities in it. Signaling machinery in the lower (red) and upper (cyan) sections, respectively, correspond to the cases wherein signal transduction occurs within a tumor cell and in its microenvironment. Communication between signaling machinery within tumor cells and between those in its microenvironment are captured respectively by golden and green arrows. Blue arrows represent paracrine communication between tumor cells and its microenvironment, and vice-versa is captured by black ones Boolean logic or a combination thereof that may predict a certain observed phenotype. However, addressing the second bottleneck is harder. For example, consider the case of cancer progression machinery such as "Invasion and Metastasis" (Figure 3). While involvement of several intracellular entities such as Snail in this machinery is well-known [112], emerging evidences suggest that mechanotransduction -signal flow triggered by mechanical cues from local microenvironment -can modulate the signaling process [176][177][178][179][180][181]. Further, recent studies show that pro-survival autophagy machinery induced by shear [182,183] may govern the key short timescale processes such as intravasation and extravasation during metastasis, seeding of which could be by polyclonal cell clusters [184]. Detailed mechanisms that regulate mechanotransduction are as yet unavailable and has attracted attention only in the recent years [177,181], which calls for improving knowledgebase on these aspects. Moreover, it is unclear how mechanical cue dynamics involving diffusion of species in physical dimensions can be incorporated in BD modeling framework, which requires variables to be strictly discrete. Combining agent-based modeling approach which permits capturing such tumor microenvironment dynamics along with the BD model of intracellular signaling could address this limitation of discrete dynamic modeling approaches [185].
Well-benchmarked individual machinery models could then be used in a plug-and-play mode to create a comprehensive, integrated BD model. Such a model should incorporate intracellular signal transduction and other processes occurring in the tumor microenvironment, and more importantly communication between them (Figure 3). Use of ROA update method further permits incorporation of heterogeneity in both inter-and intracellular interactions [186,187]. This can enable incorporating ensemble-level features such as seeding metastasis by tumor cell clusters [184], cell-to-cell abundance variability. Heterogeneity is capitalized by cancer cells for various purposes such as relapse [127,188]. Such integrated models can be used for 1. precisely predicting the dynamic orchestration of attractors representing tumor formation, growth and progression and the associated critical nodes or interactions within and outside a cancer cell, 2. deciphering dynamic reprogramming of normal long-and short-range biological processes such as metabolism within tumor cell [10][11][12][13][14][15], EMP [110,111] and 3. identifying, from the critical orchestrators and those regulating them, with significant confidence, potential candidates that enable phenotype switching for targeted therapies.
Emerging novel therapeutic strategies range from using microorganisms such as bacteria [189,190], nanoparticles [191] as carriers for delivery of drugs deep inside tumors to using genome editing technology for different cancers [192] to reprogramming metabolic pathways [193]. Since delivery agents or metabolic reprogramming could exhibit deleterious influence or undesired outcomes on the host normal or cancer cells, an integrated model can help unravel their short-and long-range (side-) effects, presence of inherent resistance to drugs or even drug-resistance induced by them. These are paramount to effective drug administration and disease management. Further, a comprehensive BD model trained with multiscale measurements [194] such as proteomics, genomics, RNAseq at scales ranging from single-cell to whole-body levels for identifying attractor phenotypic patterns can reflect the physiological state of an individual [195,196]. An omics data-tuned patient-specific model can provide insights into designing personalized therapy and novel disease management strategies [35,[196][197][198][199]. While a BD model consisting of large number of interacting nodes along with associated Boolean logic can be assembled, capitalizing it for useful purposes requires finding the one-step state transition matrix ( ). Estimation of is a computationally tedious task, especially when random order asynchronous update method is employed. Reasons for this are three-fold: 1. finding all possible FP, cyclic, and complex attractors, 2. computing large number of states needed for constructing reasonably sized STG, 3. computing frequency of permutations corresponding to every pair of states in the STG. While efficient algorithms are becoming available for finding primary attractors in large-scale Boolean networks [98], development of an effective tool for estimating is needed. We posit that this tool will strongly hinge on the ability to identify the smallest size of the STG that preserves all the necessary characteristics. We further believe use of AI-based methods could address these challenges effectively. While most BD models of cancer network consider use of FP attractors, a question arises as to what are the practical implications of cyclic attractors, if any found and how its features can be harnessed for useful purposes.
BD modeling predictions are based on qualitative relationships between the entities in the network, and thus information pertaining to subtle underlying dynamics is sacrificed. For instance, BD models fail to predict precise bifurcation or switch-like behavior dictated by the actual kinetics and associated parameters governing signal flow. This is a well-recognized severe limitation of discrete dynamic models. Using multi-valued logical variables, which increases computational complexities pertaining to evaluation of the Boolean functions, or fuzzy logic instead of Boolean logic [33,200] could address this limitation to some extent [31,48]. However, it can at best capture these subtle behaviors as state-change patterns only [201]. As discussed earlier, BD models currently cannot account for certain continuous processes such as diffusion in physical dimensions. Diffusion of various species is an important contributor in communication between a tumor cell and its microenvironment. Capturing these factors is necessary from a holistic perspective. Development of new tools that can account for these cancer-related factors while efficiently estimating is needed to build a large predictive cancer-specific BD model. A recently proposed novel approach to combine BD modeling and agentbased framework via PhysiBoSS tool [185] accounting for diffusion of species in tumor microenvironment is a small, promising step in this direction.
writing-review and editing: Shubhank Sherekar and Ganesh A. Viswanathan.