Boolean modeling of mechanosensitive epithelial to mesenchymal transition and its reversal

Summary The significance of biophysical modulators of the epithelial to mesenchymal transition (EMT) is demonstrated by experiments that document full EMT on stiff, nano-patterned substrates in the absence of biochemical induction. Yet, current models focus on biochemical triggers of EMT without addressing its mechanosensitive nature. Here, we built a Boolean model of EMT triggered by mechanosensing – mitogen crosstalk. Our model reproduces epithelial, hybrid E/M and mesenchymal phenotypes, the role of autocrine TGFβ signaling in maintaining mesenchymal cells in the absence of external drivers, inhibition of proliferation by TGFβ, and its apoptotic effects on soft ECM. We offer testable predictions on the density-dependence of partial EMT, its molecular drivers, and the conflict between mitosis and hybrid E/M stability. Our model opens the door to modeling the effects of the biomechanical environment on cancer cell stemness linked to the hybrid E/M state, as well as the mutually inhibitory crosstalk between EMT and senescence.


INTRODUCTION
Carcinomas emerge from healthy epithelial tissue because of oncogenic mutations, their growth enabled by the loss of contact inhibition of proliferation. 1,2 The weakening of cell junctions is critical for tissue invasion and migration leading to metastasis. This transformation from epithelial to an invasive, migratory mesenchymal phenotype is called the Epithelial to Mesenchymal Transition (EMT). 3,4 A rate-limiting step in the development of metastatic cancer, EMT remodels the cytoskeleton to enhance migration, disrupts the extracellular matrix (ECM) with proteolytic enzymes to promote tissue invasion, and enhances resistance to apoptosis as well as senescence. 5,6 The reverse mesenchymal-epithelial transition (MET) often aids the settling of these migratory cells into distant metastases. 7 As metastatic cancers account for 90% of cancer deaths worldwide, 8 development of effective therapeutic approaches requires predictive models of the complex interaction between contact-mediated signaling, engagement of the EMT transcriptional program, and its reversal during MET.
A closer look at cancer cells within solid tumors revealed a striking cell-to-cell heterogeneity in the expression of epithelial and mesenchymal markers. [9][10][11][12] The initial picture of an all-or-nothing transition from an epithelial to a mesenchymal state gave way to the insight that at least one but perhaps several stable cell phenotypes between these two extremes exist in vitro and in vivo. 7 These hybrid E/M phenotypes demonstrate characteristics of both epithelial and mesenchymal cells, including an ability to maintain adherens junctions while polarizing horizontally to migrate in a mesenchymal fashion. 13 As a result, hybrid E/M cells migrate and invade as clusters, 14 which gives them an advantage in seeding colonies, 15,15-17 traversing narrow capillary vessels by forming single-cell trains. 18 Indeed, hybrid E/M cells were shown to have up to 50-fold greater metastatic potential than single circulating cancer cells. 15 In healthy epithelial tissues such as blood vessels, lung, or liver epithelia, balancing birth and death in response to injury is complemented by partial EMT and migration at monolayer wounds. 3 Accompanied by proliferation, this helps wounded epithelial sheets migrate as collectives, close wounds, and use MET to reestablish their junctional barriers. Yet, the precise cellular contexts that generate hybrid E/M rather than mesenchymal cells are not well understood. 7 EMT is triggered by a series of transforming signals such as TGF-b, Wnt, Notch, Hedgehog and IL6, aided by growth factors, hypoxia and stiff/fibrous ECM. [19][20][21] These influences converge on a core transcriptional circuit rich in positive feedback that can robustly maintain both epithelial and mesenchymal states. 22 The Snai1 (Snail) As Boolean models characterize regulatory signals as expressed/active or absent/inactive, 43 they rarely deal with processes that depend on the uneven distribution of regulatory molecules in different parts of a single cell. Yet, several key aspects of EMT controlled by the cell's biophysical environment involve such asymmetry. The transition from apical/basal polarity in an intact monolayer to horizontal polarity and directed migration at lower density requires a different set of regulatory molecules to be ON at the cell's leading versus trailing edge. 44 In a recent study modeling the crosstalk between epithelial contact inhibition of proliferation, migration, cell cycle progression, and anoikis we focused on the mechanosensitive regulation of cell behaviors by ECM stiffness and a cell's neighbors. 45 Rather than building a spaceaware 3D model, we captured key features of the biology in a Boolean framework by modeling the logic of pathways that push cells past discrete thresholds: anoikis as cells loose ECM anchorage, establishment of apical/basal polarity and contact inhibition in dense monolayers, and the coupling between contact inhibition of proliferation and migration. Our model predicted a metastable migratory state in non-dividing cells at low density but stopped short of integrating these mechanosensing mechanisms with the regulation of EMT-driving transcription factors.
The importance of the biophysical microenvironment in controlling EMT is demonstrated by a series of experiments published EMT model cannot yet reproduce. These include the stiffness-dependent effects of TGF-b, leading to apoptosis on soft ECM versus EMT on stiff substrates, 46 the demonstration of full EMT in the absence of any biochemical transforming signal on stiff nano-patterned ECMs, 47 mechanical stress-mediated EMT transcription factor induction, [48][49][50] or the effect of cytoskeletal and nuclear morphology on the transcriptional regulation of EMT. 51 Yet, we lack predictive mechanistic model that can reproduce and predict the effects of a single cell's biomechanical environment on its behavior in different signaling contexts, such as varying growth stimuli or transforming signals. Moreover, current computational models only address the process of MET in response to either the withdrawal of input signals that activate Snai1, 26,27,30 or to drugs/mutations. 40,41 To address these needs, we have built on our published mechanosensitive model of epithelial contact inhibition, proliferation, and apoptosis 45 to model the process of EMT triggered by biomechanical and ll OPEN ACCESS Figure 1. The transcriptional circuit driving EMT is a tree-way switch influenced by cell junctions, migration, and growth signaling (A) Regulatory switch controlling EMT (yellow) with relevant inputs from regulators of cell cycle (purple, lavender), adhesion (light & dark blue), growth (greens), migration (pink) and apoptosis (dark red; links from the EMT switch to the rest of the network not shown). (B) Three stable attractors of the isolated EMT regulatory switch. Orange/blue node borders: ON/OFF. (C) Modular network representation of our extended Boolean model. Gray: inputs representing environmental factors; dark blue: Adhesion signals; red: TGFb signaling; green: Growth Signaling (lime green: basal AKT & MAPK, bright green: PI3K/AKT oscillations, mustard: mTORC1, dark green: NF-kB, GSK3, FoxO1); dark red:Apoptotic Switch; light blue: Contact Inhibition; pink/light orange: Migration; light brown: Origin of Replication Licensing; lilac: Restriction Switch; purple: Phase Switch; dark orange: cell cycle processes; yellow: EMT switch; /: activation; x: inhibition. See also Data S1/File S1. ll OPEN ACCESS 4 iScience 26, 106321, April 21, 2023 iScience Article The regulatory environment outside the switch can bias or dictate the state(s) it can maintain in different environments. To model this we linked the EMT switch into our contact inhibition model, 45 where it is influenced primarily by (1) mechanosensing and junction-mediated contact inhibition via b-catenin localization, (2) migration via Pak1, and (3) growth signaling via Akt1-mediated NF-kB activation and GSK3-b suppression ( Figure 1A and 1C).
Growth factors cooperate with matrix stiffness and loss of cell-cell contacts to drive reversible mechanosensitive EMT The key difference between the EMT model on Figure 1C and previous models is that it does not include the traditional exogenous drivers of EMT such as TGF-b, Wnt, Notch, Sonic Hedgehog, or IL-6. 39, 71 Instead, we focus on signals that can trigger EMT in the absence of the above signals -namely the cell's biomechanical environment. Evidence from different cell lines indicates that epithelial cells undergo an epithelial to mesenchymal transition when plated on stiff substrates. [47][48][49][50] To model this, we started by scanning all stable phenotypes of our model on a soft ECM at different cell densities and growth factor exposures (sampling algorithm that maps the model's stable phenotypes/attractors detailed in STAR Methods). Because of an inability to form stress fibers and activate YAP on soft ECM, the only stable model phenotype is epithelial -even at very low densities that preclude adherens junctions ( Figure 2A, initial state).
Plating our model epithelial cells on a stiff matrix triggered stress fiber formation, YAP activation and Rac1/ Pak1-driven migration ( Figure 2A, transition in Migration module at the start of Stiff ECM exposure). In addition, YAP led to the induction of high PI3K 110 (p110) protein expression. 72 In growth-promoting environments (serum, EGF, HGF), high expression of p110 increases the potency of PI3K/AKT signaling, which activates NF-kB and turns off GSK3b. 73 NF-kB, in turn, lead to Snai1 transcription, stabilized by low GSK3b activity. 74,75 Snai1 nuclear localization is subsequently aided by Pak1-mediated phosphorylation, a key early event leading to EMT. 76 Thus, YAP-mediated sensing of a stiff ECM works cooperatively with growth signaling to promote Snai1 activation and subsequent miR-34 repression by Snai1. Without miR-34 to lower b-catenin and neighbors to trap it at adherens junctions, nuclear b-catenin levels rise sufficiently to support high Zeb1 activation, miR-200 repression, and full EMT (Figure 2A, final state). This effect is followed by cell cycle entry, which does not require EMT but is also supported by growth signals and a stiff ECM (Figure 2A, cell cycle oscillations). In contrast, plating cells on a stiff ECM in the absence of strong growth signaling does not result in high AKT signaling, migration, Snai1 activation, or EMT (Figure S1A). An alternate route to the same mesenchymal phenotype is to start with isolated cells plated on a stiff ECM in the absence of strong growth signaling. Exposing these cells to growth factors also triggers EMT via the same cascade of signals ( Figure 2B).
Mechanosensitive EMT dynamics is robust to fluctuations in signal propagation, changes in Boolean update, and random errors in model construction The Boolean nature of the model's environmental inputs requires the assumption that the ON state of an input node such as GF_High activated its downstream signaling pathway to saturation. To model intermediate levels of growth signal exposure, we chose to stochastically toggle GF_High ON for a tunable fraction of time steps, which renders the growth signal stochastic. We assume that such intermediate signals correspond to an in vitro concentration range at which the targets of growth signals are in transition from non-responding to saturating (e.g., near the inflection point of an S-shaped response curve), a regime in which the response is most susceptible to biological noise. 77,78 Our stochastic growth input's downstream effects mimic this probabilistic, intermittent response. 79 Similarly, medium (non-saturating) ECM stiffness modeled via a stochastic ON/OFF toggle mimics a cell that tugs on the ECM with imperfect, noisy reliability with respect to downstream Rac1 activation and stress fiber generation. The comparable length scales of stress fibers and collagen fibers (or other ECM components), together with their inhomogeneous spatial organization, further supports a model with a noisy response to moderately stiff ECM. In our model, increasing the ECM stiffness under cells exposed to near-saturating mitogens (GF_High = 90%) shows a rapid increase in proliferation above 50% Stiff_ECM exposure ( Figure S1B, top left), and vice versa (Figure S1B, bottom left). Figure S1B also shows that cells on stiffer matrices spend increasing amounts of time in a mesenchymal state (top row, orange). The response to growth signals is similar, except that we predict a high prevalence of hybrid E/M cells at intermediate mitogen exposure (Figure S1B, bottom row, orange). iScience Article but oscillations such as the cell cycle may be derailed by non-biological update orders generated by asynchronous update. We have previously shown that asynchronous update applied to our cell cycle model generates a mix of behaviors from normal cycles to errors that mimic defects in mitosis (including apoptosis due to mitotic catastrophe). 80 As fully randomized update enriches for errors, we introduced a biased update where a small subset of nodes is updated in order at the start or end of list (see STAR Methods). We used this method to re-run Figure 2 simulations, averaging over a large set of independent dynamical trajectories akin to a plate of noisy cells ( Figure S2). Despite a rapid desynchronization of the cell cycle and PI3K/Akt oscillations across these runs, EMT remains robust to the change in update. Moreover, Figures S1A, S1B, S3A, and S3B portraying the same in silico experiments with different update, are qualitatively similar (note the small fraction of genome-duplication errors introduced by asynchronous update on Figure S3B, top left).
Next, we tested the model's robustness to minor structural errors by generating three distinct ensembles of random mutants ( Figure S4). First, we locked an increasing number of random nodes ON or OFF and averaged the dynamical behavior of the resulting mutant cells. Comparing Figure S4A and S4B shows that both cell cycle progression and mechanosensitive EMT are robust to 2 random knockout/knock-in errors, then struggle to maintain the mesenchymal state with 5 errors per model. For the second ensemble, we removed an increasing number of randomly chosen links ( Figure S4C) showing that the network's dynamics is robust to 5 randomly missing links. Our third ensemble with random errors in gate output show a similar degree or robustness, indicating that missing links have an equivalent effect on the dynamics as flipping the expected output of a node for one combination of inputs. Overall, the ensemble's dynamical behaviors break in a modular fashion with increasing damage, starting with a loss of cell cycle synchrony followed by stabilization of hybrid E/M, then the epithelial state.
Modeling the emergence of hybrid E/M states: Cells at the edge of a monolayer undergo partial EMT in response to stiff ECM exposure Repeating the above simulations with cells placed in cell densities akin to the edge of a confluent monolayer leads to a different outcome. As Figure 2C indicates, exposing cells that maintain adherens junctions with some neighbors but not enough to establish apical/basal polarity respond to stiff ECM exposure by undergoing partial EMT. In this case, Snai1 activation and miR-34 inhibition do not result in high nuclear b-catenin. Instead, cells stabilize in a state that matches the observed transcription factor expression and phenotype of hybrid E/M. Namely, Snai1, Twist1 and Snai2 are expressed, Zeb1 activity is moderate, and miR-34 is repressed but miR-200 is still active ( Figure 2C, EMT module). These cells are migratory and have reduced E-cadherin expression compared to dense monolayers. Yet, they maintain adherens junctions with neighbors and migrate as a group ( Figure 2C, Migration and CIP modules). These behaviors are robust under biased asynchronous update ( Figure S2C), as well as in the presence of random errors in model structure ( Figure S4, right column) Modeling intermediate levels of stiffness also shows that cells with adherens junctions spend an increasing amount of time in a hybrid E/M on stiffer substrates, but do not undergo full EMT (Figures S1C/S3C).
Of interest, we predict that although a sufficiently stiff ECM to sustain stress fiber formation is mandatory for entry into and maintenance of a hybrid E/M phenotype, sustained growth signaling is not necessary. Rather, entry into hybrid E/M only requires a short-lived boost of growth signaling ( Figure 2D) Figure S3C).
To further probe the importance of adherens junctions in partial versus full EMT, we ran a large sample of cells stochastically exposed to a matrix with relatively high stiffness (Stiff ECM = ON in 95% of time-steps) and near-saturating growth factors (95% GF High = ON). In this environment, a randomly sampled cell population will be 40% epithelial ( Figures S1D/S3D Modeling the mesenchymal to epithelial transition: in the absence of strong autocrine signaling, loss of YAP or ERK activity reverses EMT As our model lacks autocrine signals that stabilize the mesenchymal or hybrid E/M state, their stability relies on ongoing feedback between the EMT switch, growth signaling, and migration. To investigate the consequences of breaking this feedback, we simulated the response of isolated mesenchymal cells proliferating on a stiff matrix to two types of environmental change. First, we simulated their re-plating onto very soft ECM, where the maintenance of stress fibers, YAP activity, horizontal cell polarity and migration is no longer possible ( Figure 3A/S5A, synchronous/asynchronous update). This led to a loss of Pak1 signaling, nuclear exclusion of Snai1, and MET. Indeed, even partial matrix softening was sufficient to lock cells into an epithelial state, as mesenchymal entry below 50% Stiff_ECM is rare (Figures S1B/S3B, top rows).
Placing mesenchymal cells into a high cell density environment has a similar effect to soft ECM, as it also triggers MET ( Figures 3B/S5B). This, however, only occurs at densities that forbid cell spreading and stress fiber maintenance. Placing a mesenchymal cell into a neighborhood that resembles the edge of a monolayer does not alter its mesenchymal state ( Figure 3B, initial state). The presence of neighbors alone fails to trigger partial MET to a hybrid E/M state because fully mesenchymal cells lack E-cadherin expression and do not form adherens junctions. Thus, MET only occurs in our model when feedback between the EMT switch and migration is interrupted, and Snai1/Zeb1 activity is lost. Only then do cells re-express E-cadherin and establish adherens junctions.
To highlight the key role of YAP and Pak1 in the maintenance of the mesenchymal state, we modeled partial YAP or Pak1 knockdown in a 95% Stiff_ECM, 95% GF High environment. To do this, we locked YAP or Pak1 OFF for a random fraction of time-steps (tunable from 0% to 100%) and allowed them to respond to their inputs the rest of the time. The rationale for such partial knockdown is that it offers a way to model the incomplete inhibition seen in cells treated with a chemical inhibitor or siRNA silencing, as opposed a genetic knockout. As Figures  The fact that even transient interruption of migration can block Snai1 and flip the EMT switch leads to an unexpected model prediction for proliferating mesenchymal cells. Namely, during mitosis Cdk1/Cyclin B trigger a temporary loss of Pak1 as cells loosen their contact with the ECM, loose their stress fibers and round up to assemble the mitotic spindle ( Figure 2, Migration module in mesenchymal cells). 81 The resulting Pak1 inactivation does not impact the EMT switch directly, as both miR-34 and miR-200 remain repressed by Zeb1 and Snai1. That said, in cells that complete mitosis in the absence of strong growth signals, loss of NF-kB and reactivation of GSK3b compounds the loss of Pak1, all tipping the scales toward the loss of active Snai1. Thus, our model predicts that the loss of mitogenic stimuli leading to cell cycle exit is accompanied by MET ( Figures S6B/S6D). A similar fate awaits hybrid E/M cells exiting the cell cycle (Figure S6C). In contrast, cells that undergo partial EMT with the aid of a short-lived growth signal which does not push them into the cell cycle remain migratory and maintain a hybrid E/M state in low growth factor environments ( Figures 2D/S2D). This behavior explains the biphasic response to increasing growth factors observed with synchronous update (i.e., under simulation conditions where the response to a fixed level of growth stimulus is less noisy). Namely, moderate growth signals appear to increase, but saturating signals decrease the stability of the hybrid E/M state (Figures S1C versus S3C; effect blunted by asynchronous update).
The full phenotype repertoire of our model showcases the mechanosensitive routes to EMT and MET To summarize the model's ability to reproduce the cell phenotypes we expect to see in different microenvironments, we visualized our model's phenotypes (attractors) across all input-combinations ( Figure 3E; for sampling algorithm and complete list of stable phenotypes across all environments see STAR Methods and Table S1). To this end we created a coordinate system in which each axis represents an independent environmental factor. 45 iScience Article picture to represent each stable cell phenotype. For example, apoptosis (blue apoptotic cell) is the only stable state in environments with no growth factors -the leftmost position along the x axis representing growth signals. In environments that also support non-apoptotic attractors, there is more than one phenotype-symbol (e.g., apoptotic cells and green epithelial cells in growth factors with soft ECM; Figure 3E, bottom plane).
As Figure 3E indicates, a mesenchymal phenotype in the absence of external or autocrine drivers of EMT is only stable on stiff ECM, at medium to low cell density, in the presence of strong mitogens (red stars). Weaker, non-proliferative growth stimuli still permit the maintenance of hybrid E/M, provided they can stretch on stiff ECM (moderate/low density, orange triangles). If they become isolated when exposed to strong growth signals, these hybrid E/M cells undergo full EMT to a mesenchymal state (lone red star).
In contrast, both soft ECM and high cell density force cells into an epithelial state (green boxes), whereas loss of all growth signals leads to apoptosis (light blue images).
Our model makes intriguing predictions about the bistable mix of epithelial and hybrid E/M cells that coexist under low growth signaling. Namely, we predict that short spikes of growth signaling that fail to trigger cell cycle entry can push epithelial cells into hybrid E/M at a monolayer's edge. This hybrid E/M state is maintained with a migratory state, until interrupted by mitosis, soft ECM exposure, or high cell density. In contrast, transitions triggered by sustained growth signals trigger cell cycle entry and lead to reversible EMT ( Figure 3E, EMT: red/orange arrows; MET: mustard arrows).

Autocrine TGF-b signaling locks in EMT triggered by growth signaling on stiff ECM and blocks mesenchymal cells proliferation
Next, we set out to model the crosstalk between the biomechanical environment and autocrine TGFb signaling required to sustain the mesenchymal state in vitro. 25,26 To this end we added a TGFb signaling cascade, as well as a secreted TGFb node to the EMT switch ( Figures 4A and S7). In this new module, externally applied or secreted TGFb can bind to TGFb receptors I and II, which in turn activate Smad 2/3/4. 83 Smads drive EMT by promoting nuclear translocation of b-catenin, and inducing Snai2, Lef1, and HMGA2 expression. HMGA2 further induces Snai1, Snai2, and Twist. In addition, we included TGFb 0 s paradoxical effect on the MAPK cascade with PI3K/AKT activation, anti-proliferative p15 induction, and pro-apoptotic BIM/BIK induction and BCL-X L repression ( Figure 4A; see Data S1/File S1 for details). To create autocrine positive feedback, we included TGFb secretion downstream of nuclear b-catenin/TCF4 (repressed by miR-200).
Adding an autocrine loop to a Boolean regulatory model makes the often-unexamined assumption that a single cell is capable of secreting sufficient signal (in our case TGF-b) to drive a saturating response (i.e., pathway is ON rather than OFF). As a result, the effects of autocrine loops can be easily overestimated in Boolean models, indicating irreversible lock-in of EMT in situations where cells in vitro show metastable EMT because of partial autocrine activation. 84 To mitigate this, we included an additional input named Self_Loop. When set to 0, TGFb secreted by an isolated cell is too weak to trigger signaling, whereas Self_Loop = ON represents saturating levels of autocrine TGFb activation. Tuning this input to intermediate values allows us to model a range of autocrine signal strengths in sparsely plated cells. To further account for the possibility that cells with mesenchymal TGFb-secreting neighbors have stronger autocrine signaling, 85  . This hybrid E/M is triggered by short bursts of growth factors that do not persist long enough to trigger cell cycle or lock in the autocrine signaling loop (Figures S11A/S12A versus S11B/S12B). In contrast, short-lived exposure to a stiff ECM either fails to trigger sustained EMT (Figures S11C/S12C) or leads to TGF-b-mediated apoptosis (Figures S11D/S12D) (see next section).  (1) increasing Stiff ECM exposure in the presence of 95% saturating growth stimuli (top two) versus (2) increasing growth factor exposure on 95% stiff ECM (bottom two). Proliferation label: rate of normal cell cycle completion (blue) versus G2 / G1 reset (orange, not observed), aberrant mitosis (green, not observed), or failed cytokinesis (red, not observed) followed by genome duplication, relative to the minimum cell cycle length (21 time-steps), shown as stacked bar charts. EMT label: fraction of time spent in a mesenchymal (orange) versus epithelial (blue) state.
(E) Response of isolated cells on 85% Stiff ECM exposed to 85% saturating growth stimuli to increasing inhibition of autocrine signaling (Self-Loop knockdown). Top: rate of normal cell cycle completion (blue) versus G2 / G1 reset (orange), aberrant mitosis (green, not observed), or failed cytokinesis (red, not observed) followed by genome duplication. Bottom: fraction of time spent in a mesenchymal (orange) versus epithelial (blue) state. Total sampled live cell time: 100,000 steps; synchronous update; initial state for sampling runs: isolated epithelial cell in low mitogens on a soft ECM. See also Figures S7-S14 and Data S1/File S1. iScience Article Next, we weakened autocrine TGF-b from 1 to 0 in isolated cells exposed to an environment that favors but cannot fully force mechanosensitive EMT; namely 85% growth signals and 85% stiff ECM. As Figures 4E and S8E indicate, even a slight weakening of the autocrine loop decreases the time cells spend in a mesenchymal state, restoring proliferation instead. Its full knockout leads to a mix with majority epithelial versus mesenchymal cells (bottom panel).
In summary, our model shows that growth signals and stiff ECM work together to drive EMT in sparsely plated cells, a transition that is subsequently maintained by strong autocrine TGFb signaling.
In contrast to full EMT in sparse cells, partial EMT at the edge of a monolayer does not result in TGFb secretion in our model. As a result, the extended model's behavior is the same as that of the previous one ( Figures 2C, 2D versus S13 and S14). Namely, matrix stiffness and growth signaling induce reversible hybrid E/M and proliferation in cells at the edge of a monolayer ( Figures S13A-B), but sustained migratory hybrid E/M in cells exposed to sub-proliferative bursts of growth signaling ( Figure S13C). We predict that moderate levels of growth factor exposure on stiff ECMs leads to the most time spent in a hybrid E/M state (Figure S14B, bottom), as cells reset from a hybrid E/M to an epithelial state every time they round up for mitosis. Although short bursts are necessary for partial EMT, sustained exposure triggers a conflict between proliferation and migration (hybrid E/M) and requires a new partial EMT transition after each division. Knockdown of E-cadherin breaks the tug of war by favoring full EMT at the expense of proliferation ( Figure S14C).
Next, we set out to text whether altering the biomechanical environment can trigger MET in the presence of strong autocrine TGFb signaling. After verifying that mesenchymal cells remained stable at high cell density (Figures 5A/S15A, left), we exposed them to soft ECM.  S15C). As the density increases, mesenchymal cells shift to an epithelial -hybrid E/M mix generated by crowded cells undergoing full MET, followed by partial EMT whenever space opens ( Figure S16) -an event that disappears in a confluent epithelium. These behaviors are independent of update, though the asynchronous version of our model has lower proliferation and E-cadherin knockdown thresholds ( Figure S14), but higher resistance to cell density-mediated MET ( Figure S15). In addition, the model's behaviors are robust to several simultaneous random mutations ( Figures S9 and S10).

Modeling anoikis resistance: mesenchymal cells lose anchorage dependence
A concerning consequence of EMT in the context of metastatic cancer is loss of anchorage dependence, or resistance to anoikis in mesenchymal cells. 57 To test whether our model can reproduce this, we simulated loss of anchorage with an epithelial cell in low growth factor conditions ( Figure S17A) versus a mesenchymal cell sustained by autocrine TGFb signaling in the complete absence of growth factors ( Figure S17B). As expected, the epithelial cell underwent apoptosis whereas the mesenchymal one did not (Figures S17A-S17C). Overexpression of either AKT or ERK rescued epithelial cells from anoikis ( Figure S17D, left). In mesenchymal cells, however, even complete knockout of AKT was not sufficient to induce anoikis, as TGF-b-mediated ERK activity could maintain survival. In contrast, ERK knockdown could re-sensitize these cells to anoikis ( Figure S17C, right).

External TGFb triggers EMT and cell cycle inhibition on stiff ECM, but apoptosis on soft substrates
TGFb was shown promote EMT on stiff matrixes, but cause apoptosis on soft ECMs in the same epithelial culture. 46 To test whether our model can reproduce both outcomes, we simulated cells at the edge of a monolayer (medium density) on soft versus stiff ECM and exposed them to external TGFb. Indeed, our model epithelial cells underwent apoptosis on soft ECM, an event driven by a combination of BIM overexpression, BCL-X L downregulation, and critically, lack of both ERK and high AKT activity ( Figures 6A/S18A). The latter is caused by weak activation of both mitogenic pathways in the absence of strong focal adhesions and stress fibers. 86,87 In contrast, on stiff ECM TGFb signaling activates both pathways, counteracts its own apoptotic effects, and triggers EMT instead (Figures 6B-C and S18B-C, low/high growth factor env.). As cell density increases, the ability of TGFb to trigger EMT weakens. 59  iScience Article responding to TGFb, even if its receptors are accessible to external TGFb ( Figures 6D/S18D). In a microenvironment with saturating growth signals and 50% external TGFb, most cells on soft ECM have an epithelial phenotype and are prone to apoptosis, but also to occasional cell cycle entry ( Figure 6E, 20% of saturating proliferation rate on 20% stiff ECM with synchronous update). As matrix stiffness increases, the rate of both proliferation and apoptosis drops, giving way to EMT. Our asynchronous update results are qualitatively similar, though the model is substantially more sensitive to both TGFb-induced apoptosis and EMT (Figure S18E). In addition to its robustness to update, the model's behaviors are largely unchanged by mutations up to the lockdown of 5 random nodes ( Figure S19A) or removal of R10 links ( Figure S19B).  To test the effects of weakening autocrine TGFb signaling, we started with cells exposed to near-saturating growth signals and scanned their responses as a function of increasing ECM stiffness ( Figure S20A, y axis) at densities increasing from isolated to that of a monolayer's edge ( Figure S20A, x axis). As Figure S20A illustrates, even a 10% loss of autocrine signal strength in isolated mesenchymal cells (modeled as Self_Loop = OFF for a random 10% of time steps), restored proliferation on stiffer ECMs at all densities, rather than only in hybrid E/M cells capable of maintaining adherens junctions at a monolayer's edge ( Figure S20A, top). Intriguingly, our model predicts a small uptick of cell cycle errors such as slippage from G2 back to G1  iScience Article and failure to complete cytokinesis following telophase, both more common among sparce cells (2-3% of all divisions). An increase of apoptosis with the weakening of autocrine signaling is also predicted on stiff ECMs at very low density, at a slightly higher rate than cell cycle errors. Meanwhile, weakening the autocrine loop severely destabilized the mesenchymal state at very low density, triggering MET on all but the stiffest ECM ( Figure S20A, bottom 3 panels). Next, we repeated the above scan with cells on 95% stiff ECM, scanning their responses as a function of increasing mitogen exposure ( Figure S20B, y axis) at increasing densities ( Figure S20B, x axis). In this case a weaker loop increased the prevalence of cell cycle errors even more ( Figure S20B, top panels), but more importantly, reduced the stability of mesenchymal states in favor of hybrid E/M.
To summarize our results, we recreated the map cell phenotypes our model creates in different microenvironment-combinations, accounting for the availability of external TGFb as a fourth dimension (Figure 7; complete list of detected model attractors, including those with a density-dependent TGFb loop and Trail = ON, are summarized in Table S2 and listed in Table S3). In contrast to the model with no autocrine TGFb signaling ( Figure 3E), the quiescent mesenchymal phenotype is stable in the absence of external TGFb in all environments except on very soft ECMs at very high density -assuming a strong autocrine signaling loop (Figure 7, left panel, red stars). In contrast, ongoing and saturating external TGFb eliminates the hybrid E/M phenotype characteristic of monolayer edges, whereas the epithelial phenotype is only stable at very high density (Figure 7, right panel). Of interest, our model predicts that cells at this high density do not undergo apoptosis in response to TGFb; rather they use basal AKT activation by TGFb as a survival signal even in the absence of all other mitogens, or in the absence of anchorage -though they cannot survive the absence of both (Figure 7, right panel, green boxes).

DISCUSSION
Modeling the signaling networks responsible for EMT has been an ongoing focus for the last decade. 22 iScience Article control partial versus full EMT, with or without transforming signals, are not adequately linked to molecular mechanisms responsible for each outcome. To address this, we first built a 136-node Boolean regulatory network model of partial as well as full EMT, driven by a cell's biomechanical environment in the absence of transforming factors such as TGF-b, Wnt, Notch, Hedgehog or IL6 ( Figure 1C). Instead, EMT in this model is triggered by mitogenic signals on stiff ECM, together with a lack of contact inhibition leading to migration. Rac1/PAK1 activation and PI3K/AKT-mediated NF-kB signaling kickstart EMT by increasing the expression, localization, and stability of SNAI1. In line with most published models, our network can maintain distinct epithelial, hybrid E/M and mesenchymal states (Figures 1B and 3E). Unlike other models, we also track the coordination between EMT, migration, cell cycle progression and apoptosis. With the aid of mechanosensitive inputs such as low density, our model can reproduce EMT driven by strong mitogens such as EGF on stiff ECMs (Figures 2, 4, S1B, and S1C), 52,53 with the observations that loss of E-cadherin can promote EMT (Figures S1D and S14C). 54,55 Our extended 150-node model endowed with a TGFb signaling module ( Figures 4A and S7) reproduces the need for a potent autocrine TGFb loop to maintain the mesenchymal state on soft ECM ( Figure 4B), in the absence of strong mitogens ( Figure 4C), or at high density ( Figure 5A, left). Although this autocrine loop is a common feature to most published EMT models, here we also reproduce TGFb-induced apoptosis on soft ECMs versus EMT on stiff ones ( Figures 6A and 6B), 46 TGFb-induced inhibition of proliferation in mesenchymal cells ( Figure 6C), 56 the lack of TGFb-induced EMT in very dense monolayers ( Figure 6D), 59 and mesenchymal anoikis resistance sustained by autocrine TGFb ( Figure S17). 57,58 This board repertoire of behaviors is a direct result of a modular modeling strategy by which we built the EMT transcriptional network as an extension of a model series that started with the detailed control of cell cycle progression, 82 its coordination with apoptosis, 80 then examined the effects of contact inhibition and loss of anchorage on both phenotypes. 45 In addition to reproducing known cell behaviors, our model offers several mechanistic, experimentally testable predictions.
1) In our model, the density of adherens junctions around a single cell is a biomechanical signal that biases EMT toward hybrid E/M versus a mesenchymal outcome. Thus, we predict that relatively isolated cells in sparsely plated epithelia exposed to EGF on a stiff matrix are more likely to undergo full EMT, whereas cells at the edge of a monolayer or with a similar contact-surface stop at hybrid E/M instead ( Figure 2). This can be tested with single-cell resolution live-imaging experiments tracking junctional E-cadherin localization with b-catenin levels/localization or miR-200 activity (junctional versus cytosolic versus nuclear b-catenin live reporters are descried in 91 ; live E-cadherin probe on an orthogonal imaging channel in 92 ; live miR-200 probes can be designed according to 93 ). Specifically, cells imaged at increasing density on stiff ECM (>100 kPA) in standard growth-promoting media are predicted to generate three distinct subpopulations as a function of junctional E-cadherin: (1) Cells with low nuclear b-catenin (or high miR-200) and high E-cadherin surrounded by neighbors (epithelial), (2) cells with moderate nuclear b-catenin, substantial junctional E-cadherin and some free edges (hybrid E/M), and (3) a high nuclear b-catenin group with no E-cadherin and few/no neighbors (mesenchymal). We further expect their balance to shift from mostly mesenchymal at very low density toward hybrid E/M, then to epithelial in very dense monolayers. In contrast, data showing no correlation between local cell density (single cell neighborhood) and nuclear b-catenin and/or miR-200 expression would refute our predictions and necessitate a reevaluation of our understanding of mechanosensitive EMT.
2) We predict that short, intermittent mitogen bursts that lack the strength to induce cell cycle entry can drive sparsely plated epithelial cells into a metastable migratory hybrid E/M state, which is maintained by uninterrupted migration (Figures 2D, S11A, and S13C). As a corollary, we predict a biphasic response to growth signal exposure where moderate signaling results in the largest fraction of hybrid E/M cells ( Figure S14B), whereas saturating mitogens interrupt migration and reset cells at each mitotic rounding (Figures S11B and S13B). This could be tested with an experimental setup involving transient growth factor (or serum) exposure followed by washout, repeated with increasing exposure lengths and/or growth factor concentrations. 79 Cells would be live-imaged before, during and after this pulse to track their migration speed, nuclear b-catenin levels, and division. Our model predicts that shorter/weaker pulses preferentially generate persistent migratory cells with moderate nuclear b-catenin (higher than the pre-pulse baseline), whereas longer/stronger pulses result in divisions, ll OPEN ACCESS 16 iScience 26, 106321, April 21, 2023 iScience Article followed by an immediate drop in daughter cell migration. We also expect nuclear b-catenin to be highest in cells imaged after exposure to moderate pulse length/strength (imaged around the time of appearance of daughter cells). Here, it is possible that the migration triggered by growth factor pulses is not as persistent in live cells as our deterministic model predicts. Yet, observing some persistent mitigation and/or sharp drops in speed following mitosis would still validate our model. In contrast, observing a loss of migration immediately on growth factor withdrawal would require us to reevaluate the bistability of our migration module.
3) We predict that the biomechanical environment alone can trigger MET by blocking stress fiber formation and muting ERK/AKT signaling, whereas Pak1 hyperactivation can reestablish hybrid E/M ( Figure 3). Furthermore, soft ECM and high cell density cooperatively break the autocrine TGFb loop and trigger MET via contact induced Sprouty 2 activation ( Figure 5). This could be tested by triggering full EMT with TGFb in a sparse population, re-plating them at very high density on a soft ECM (0.5 kPa), then imaging the reappearance of junctional E-cadherin followed by phospho-AKT or phospho-ERK staining. 94 Here we expect to see MET, preferentially in the densest areas of the plate. Repeating the imaging in the presence of Pak1 overexpression and Sprouty 2 knockdown would further test the accuracy of our model's molecular mechanism driving mechano-sensitive MET.

4) In line with previous model predictions and observations about the asymmetry of EMT versus MET, 27,95 our model does not undergo a direct transitions from mesenchymal to hybrid E/M state.
Rather, model environments that favor hybrid E/M over the mesenchymal state induce full MET followed by partial EMT ( Figure S16). This can be tested by live-imaging nuclear b-catenin or miR-200 with junctional E-cadherin in a series of experiments like the one above, at increasing ECM stiffness (0.5 kPA to 100 kPA) and increasing cell densities (sparse to dense enough to forbid stress fiber formation). According to our model, a random selection of mesenchymal cells plated onto moderately stiff ECMs at medium density should undergo an abrupt drop in nuclear b-catenin followed by the re-expression of E-cadherin, appearance of adherens junctions, and a moderate b-catenin re-accumulation at the edges of dense clusters.

5) Although
TGFb-induced EMT requires stiff ECM, we predict that the mesenchymal state can be maintained on soft ECMs by autocrine TGFb ( Figure 4B). The above setup with very low density and soft ECM should show a robust preservation of the mesenchymal state (no drop in nuclear b-catenin), whereas the same experiment in the presence of a TGFb signaling inhibitor (e.g., a dominant negative soluble receptor that weakens the autocrine loop) is expected to show MET.
6) Finally, we predict that at low-to medium cell density and medium-strength growth signals, autocrine TGFb signaling that does not saturate the TGFb pathway (does not reliably guarantee its downstream effects) will generate a dynamic mix of mesenchymal and hybrid E/M cells, the latter of which proliferate and thus continually reset to undergo EMT anew ( Figure S20). This is accompanied by a small uptick in apoptosis and cell cycle errors that generate polyploidy (specifically, slippage from G2 and/or telophase to G1) in the absence of any mutations related to checkpoint control ( Figure S20). This can also be tested with the setup above using stiff ECMs but a range of growth factor levels. Live imaging that tracks key cell state markers (nuclear b-catenin, junctional E-cadherin, miR-200) with division and apoptosis could identify the hybrid E/M / epithelial transition following mitosis, as well as the link between non-saturating autocrine signaling and apoptosis. Finally, testing the DNA content of the cells at the end of these live imaging runs could test our cell cycle error predictions.
Rather than our ability to reproduce known cell behaviors, we believe that the above predictions represent a more important contribution to characterizing the roles of partial versus full EMT in non-oncogenic processes such as wound healing. Testing them could boost the credibility of our model in predicting the effects of targeted interventions. Alternately, failed predictions would point to unexpected gaps, helping us further probe the behavior of cells in different biomechanical environments. Ultimately, correct integration of mechanosensing and transforming signals will be required for building predictive cell population models that follow the 2D or 3D processes of wound healing, the initiation of cancer invasion, or the establishment and growth of metastatic lesions.
As detailed under Limitations of the Study below, it is possible that the TGFb-mediated cell cycle block is too strong in our model, iScience Article mesenchymal cells despite potent autocrine TGFb. Yet there is substantial experimental evidence for the tug of war between proliferative potential and the migratory mesenchymal state. 7 For example, hepatic oval cells capable of regenerating the liver when hepatocyte proliferation is compromised where shown to have a hybrid E/M phenotype. 96 Moreover, cancer stem cells in squamous cell carcinoma and ovarian cancer exist as a heterogeneous mix of proliferative hybrid E/M and mesenchymal states, and they dynamically regenerate this mix from either subpopulation. 97,98 Indeed, in vivo metastatic tumor growth in ovarian and breast cancer is primarily because of hybrid E/M cells. 101,102 We can model this division/migration conflict in environments that can induce and maintain EMT biomechanically, and thus do not rely on uninterrupted autocrine TGFb signaling ( Figure S20). What our model currently does not account for, however, is the connection between the hybrid E/M state and the regulatory network responsible for a stem cell-like state, observed in all forms of EMT (i.e., developmental, regenerative, and pathological). 30,99 Integrating a stemness module, with an explicit accounting for the likely context-dependent role of 'phenotypic stability factors' such as OVOL and GRHL2 known to stabilize the hybrid E/M state and promote stemness, would be an important next step in delineating the similarities and differences between EMT in different tissue and disease contexts. 30,100 In addition to limitations, it is important to delineate the boundaries of our model; cell behaviors linked to EMT that we explicitly do not address. Two of these stand out as key to our ability to integrate our regulatory network into multiscale models of tissue behavior. First, our model does not address the patterning effects of Notch signaling in a population of cells undergoing EMT. Notch not only promotes EMT by inducing Snai1, 39 but also generates distinct patterns between neighboring cells depending on the ligands expressed on neighbors. 37 Cells expressing Delta ligands promote lateral inhibition patterns that alternate Delta High /Notch Low cells with the opposite phenotype. This is central to the endothelial to mesenchymal transition, as it generates alternating patterns of tip cells and stalk cells required for sprouting angiogenesis. 71,101 In contrast, the Jagged ligand promotes lateral induction, leading to a strong clustering of similar cells. 102 Previous computational models of Notch signaling showed that lateral induction via Notch-Jagged can generate intermingled clusters of epithelial, mesenchymal, and hybrid E/M cells akin to the patterning seen in tumors. 37 This pattern formation is heavily influenced by the tissue microenvironment such as inflammatory signals (IL-6, TNFa), 103,104 Wnt and Hedgehog signaling, 105,106 as well as hypoxia, 107 most of which promote EMT and induce Jagged.
The second direction that merits our attention is the relationship between EMT and senescence. The process of senescence makes short-term use of damaged cells by shutting down their ability to divide, but using them to boost survival, proliferation and EMT in their surroundings. 108 Although senescent cells trigger their own immune clearance, their slow accumulation is accelerated by chronic damage such as inflammation or oxidative stress and leads to tissue aging. 109 There is increasing evidence that senescence and the mesenchymal state are mutually exclusive in individual cells. 110,111 Mesenchymal cells are largely protected from senescence, whereas senescent cells do not respond to signals that trigger EMT. 112 During localized healing of an epithelium, this creates a useful distribution of labor: damaged cells go senescent but induce division and EMT in neighbors. These mesenchymal cells are then protected from senescence until they reestablish an intact epithelium. As neither senescent nor mesenchymal cells contribute to barrier formation, this can accelerate dysfunction in aging tissues. 113 Worse, in epithelial tumors this crosstalk protects invasive but damaged cells from senescence, aiding metastasis. 114 Yet, the two phenomena are usually studied in isolation, and the relationship between hybrid E/M cell states and senescence remains unexamined. Building a predictive molecular model of the mechanisms that forbid their coexistence and probe the transition from healing to the pathological roles of EMT and senescence is a key next step that benefits from our current integration of EMT with cell cycle control, apoptosis, and contact inhibition. Ultimately, we envision these single-cell network models as the engines that drive 3D tissue simulations, able to predict effective interventions that can nudge diseased tissues back to health.

Limitations of the study
Our models have several limitations that can be addressed by future refinement; some related to the way our previous contact inhibition model's environment was built, some specific to EMT and TGFb signaling. As detailed in, 45 we use a crude approximation of the spatial asymmetry of Rac1/Pak1 activity at the leading edge and a minimal migration module. Yet, feedback between cell contractility and ECM stiffness -itself modulated by the migrating cell -was shown to render the behavior of invading cells nonlinear, such that intermediate stiffness is optimal for high-speed migration. 115  iScience Article contractility to move faster and invade more effectively than mesenchymal cells. 116 The migration module of our model needs significant refinement to address these phenomena before it is suited to model collective behavior such as wound healing in 2D. Other areas of improvement inherited from 45 hinge on the limits of how well Boolean models capture gradual changes in ECM stiffness and cell density (detailed in Supplementary Methods and 45 ).
Another limitation relates to our Boolean approximation of the biophysical microenvironment. Namely, our model draws explicit distinctions between ECMs that do or do not support sufficient adhesion for survival, then another one between ones that do or do not aid force generation required for both migration and cell cycle entry. Similarly, our cell density nodes draw a sharp distinction between isolated cells that cannot form any adherens junctions and models slightly higher densities with intermittent junction formation. Another boundary is between cells that have room to polarize horizontally versus those that are so dense that they cannot form stress fibers. Although our model can modulate stiffness and/or density within these two regimes by toggling their inputs ON/OFF with a tunable probability, it does not allow overlap between the two regimes. Such artificial boundaries between district cell response regimes could lead to mismatch with experimental data.
Finally, a potential limitation related to our integration of the EMT switch with the cell cycle is that our model's mesenchymal cells do not show anchorage-independent growth. 7,117 Although we reproduce anoikis resistance, our cell cycle control relies entirely on spreading on a stiff matrix. A deeper examination of how cancer cells acquire anchorage-independent growth, and whether there are pre-requisite mutations that could drive such growth in our current model, requires further work. A related potential problem is that our mesenchymal cells with strong autocrine TGFb loop do not divide at all, owing to a TGFb-mediated cell cycle block. Instead, in environments that do not force full EMT (sub-optimal matrix stiffness, non-saturating autocrine TGFb signaling, moderate density, non-saturating growth signals) our model predicts a fluid coexistence of mesenchymal and hybrid E/M cells, the latter of which divide with high frequency but then briefly reset to an epithelial state before undergoing EMT once more. It is possible that our model's TGFb-mediated cell cycle block is too strong, in that fully mesenchymal cells with sufficiently potent autocrine signals to keep them mesenchymal may nevertheless divide.

STAR+METHODS
Detailed methods are provided in the online version of this paper and include the following:

AUTHOR CONTRIBUTIONS
E.S. and E.G. designed the first version of our Boolean EMT switch and started its integration with the model in 45 , leading to the first mechanosensitive version (data curation; investigation, formal analysis, validation). M.H., A.B., and I.Z. tested and refined the model to converge in the version included here (data curation; formal analysis, validation). E.R.R. supervised the project, acquired funding, developed software, analyzed/refined the final models, ran final simulations, and wrote the paper with revision and edits from E.S.

DECLARATION OF INTERESTS
The authors declare no competing interests.

INCLUSION AND DIVERSITY
One or more of the authors of this paper self-identifies as an underrepresented ethnic minority in their field of research or within their geographical location. One or more of the authors of this paper self-identifies as a gender minority in their field of research. One or more of the authors of this paper self-identifies as a member of the LGBTQIA+ community. One or more of the authors of this paper self-identifies as living with a disability.

Defining and refining the model's Boolean regulatory logic
Boolean gates were inferred from literature and/or chosen based on reported effects of molecular input combinations. These gates are further tuned, within the constraints of experimental literature, to a) generate model behaviors that match all known cell behaviors within the scope of the model and b) eliminate phenotypically similar attractor states within identical environments where known expression/activity pattern of the molecules that differ between these attractors single out one of these states as the most biologically accurate. This fine-tuning is most common for complex gates with several inputs, where preserving the known activating/repressive influence of individual inputs, along with known synergy or independence between pairs, is insufficient to fully specify the Boolean gate. In such cases, the model's dynamical behavior across all environments and environment changes is used to further constrain the gate. Namely, most gate options allowed by experimental data specific to its output node can be eliminated because they produce non-biological behaviors in some context (e.g., non-biological stable state, loss of an expected state, incorrect transition dynamics). Fixing these behaviors occasionally necessitates the inclusion of additional regulators (node) or interactions, always drawn from the experimental literature.
During this initial construction phase we use synchronous update, 121 where every node of the network updates synchronously in every time-step. The advantage is that it renders the dynamics of the system completely deterministic, allowing us to account for the precise role of each molecular species at every causal step along a biological process and compare it to experimental data. By building a synchronous Boolean model first, we could identify the molecular causes of behaviors that deviated from known cell dynamics and iteratively fix the model by altering/expanding the logic gates leading to non-biological behavior. Our analysis of the model's robustness to errors indicates that the average dynamical behavior is difficult to alter by a random change to a few logic gates ( Figures S4, S9, and S10).

Storing Boolean models in Dynamically Modular Model Specification (.dmms) files
To keep track of all model metadata such as its modular structure, experimental data justifying each node, link and gate, module structure, node visualization data, etc., we collated all model-specific information into human-editable and machine-readable files (Data S1/Files S2-S3; do not rename). These . iScience Article parses Data S1/File S2, exports the model in BooleanNet format (Data S1/Files S4-S6), and generates a folder of truth tables and graphical information read by ReganLabBooleanSims (C++ code) used for running virtual experiments.

Generating and storing network visualizations
Model files in .dmms format allow us to store node layout and color to ease the iterative visual rendering of large networks. The command generates a .gml graphical layout using the coordinates and node colors stored in the file, if any. This file can be read and organized by the yED graph visualization software, 122 and re-saved as .gml (Data S1/Files S7, S8). Assuming that the network structure is completely unchanged by these edits, dynmod can read and update coordinates and colors in the.dmms file to preserve the layout:

Generating supplementary tables with detailed model justification
The command generates a formatted LaTeX document (.tex) and bibliography file (.bib) with a series of tables that organize the model's nodes by modules and use metadata from the.dmms file to detail the biological role of each node, the molecular mechanism of action for each link, and the assumptions made in constructing large Boolean logic gates. These can be typeset using LaTeX to generate the large, referenced table in Data S1/File S1.

Sampling the state space of a Boolean model to map its synchronous attractors
We used synchronous Boolean modeling to find every stable phenotype and/or oscillation (attractor) of the model. 45,80,82 To do this we use a stochastic state space sampling procedure adapted from, 123 implemented in dynmod and ReganLabBooleanSims. First, we implemented a noisy version of synchronous Boolean dynamics, in which each regulatory node is affected by a small amount of noise, p noise , in every time-step, resulting in the incorrect output. 124 We then use the noisy dynamics to aid our sampling procedure by starting the network from a random initial condition and simulating a time-course of N series noisy time-steps. As the model generates its noisy dynamical trajectory, the algorithm pauses at each state it visits to find the attractor basin this state would fall into if the dynamics were to continue in a deterministic fashion. In addition, the version implemented in ReganLabBooleanSims also scans the immediate neighborhood of this state by enumerating every state the system could reach from the current one via a single node-state flip and identifying their attractor membership (deterministic dynamics). 82 This allows the algorithm to access parts of the state space the noisy dynamics might never go near and find smaller basins. As a result, the algorithm in ReganLabBooleanSims is slow on random Boolean networks with large numbers of small basins. By contrast, our model's attractor basins, representing robust phenotypes, are typically large and rapidly found. The full algorithm descried in 82 tracks the visitation probability of each state, basin, and transition (not used here). To run the noisy trajectory from multiple random initial conditions, our code automatically detects all environmental input nodes, partitions the full state space of the model into sub-spaces corresponding to each unique environmental node state-combination, and samples each subspace from N rnd random initial conditions. First, we used dynmod for an exhaustive attractor search with the parameters:

OPEN ACCESS
ReganLabBooleanSims outputs a series of pdfs that show all detected attractors organized by the external environments they are stable in (e.g., Data S1/File S9). These are marked by a numerical ID, as well as a quick visual indicator of their phenotype. The numerical IDs are not preserved between independent sampling runs but remain unchanged if ReganLabBooleanSims can read the results it saved from a prior sampling from the model folder (this also saves repeated samplings of unchanged models). The phenotype indicators are orange circle = quiescent; red oval = cell cycle; orange circle with dark border = apoptotic; orange circle with vertical line = quiescent, duplicated genome; E/H_Migr/M = EMT module in epithelial/ hybrid/mesenchymal state. Lines indicate transitions in response to environmental change. This rudimentary attractor visualization can be used to specify initial phenotypes and environments for further experiments described below. Figures 3E and 7 offer a simplified illustration of the attractor repertoires of our models by omitting all apoptotic attractors in conditions with ECM = 0 and Trail = 1 ( Figure 3E), or Self_Loop = 0 and Trail = 1 (Figure 7). In addition, we omit all non-cycling live attractors with double DNA content (4N), but otherwise identical to their counterpart with 2N DNA (described in 80,82 ). Finally, omission of the Trail = 1 condition from Figure 7 also hides 24 quiescent mesenchymal attractors that appear resistant to Trail exposure (not discussed).

Automated isolation of a subnetwork from a larger (multi-switch) Boolean network
ReganLabBooleanSims can automatically isolate and examine the dynamical behavior of any network module, such as the EMT switch ( Figures 1A and 1B). To this end, we have previously developed an algorithm that defines the Boolean gates of nodes when they lose some of their incoming connections. 82 The main goal was to optimally preserve the regulation of a node by its remaining inputs. Briefly, whenever a subset of inputs is removed from a Boolean gate, the algorithm assumes that they are frozen into either an ON or an OFF state. To preserve the dynamical influence of the remaining nodes, ReganLabBooleanSims finds one of the 2 k possible combinations of frozen inputs such that: a) all remaining input nodes are functional,

Biased asynchronous update
Asynchronous update schemes change the state of one node at a time and use this new state as they update the node's targets. Random order asynchronous models, the scheme adopted here, update every node in every time-step, but they do so sequentially in a random order re-shuffled before every step. Asynchronous update schemes are favored in biological modeling, as they simulate the unfolding of the same regulatory process along many slightly different paths, each with different likelihood, 125 mimicking the stochasticity present in vitro. Moreover, asynchronous update eliminates potential artifacts of synchronous modeling, such as behaviors that rely on perfect and deterministic coordination of parallel signals. That said, they also generate biologically non-realistic sequences of molecular events by failing to follow up on the effects of short-lived signals that live cells reliably respond to. As a result, we use a hybrid update scheme developed in 80 where the update order of most, but not all nodes is random (biased asynchronous update). Specifically, we update 12 of the 150 nodes either first or last, depending on their correct state, as follows (for detailed justification for each node, see 80

Time courses with reversible environmental change
To test the model cell's responses to reversible changes in their environment, we choose a specific phenotypic state in a desired initial environment as an initial condition, then specify the length of time we wish to flip one environmental variable. ReganLabBooleanSims shows the model's state or behavior for 50 timesteps, flips the environment node for the designated length of time, then flips it back and further follows the network's dynamics for a total of 400 time-steps. These in silico experiments for synchronous and biased asynchronous update are specified in the Virtual_Experiment_List.txt file as (bias rules are currently hard coded in ReganLabBooleanSims): Modeling an ensemble of mutant networks (robustness to network errors) To test the robustness of our models to random mutations and/or errors in construction, we generated three distinct types of mutant network ensembles. In the first ensemble, each network has n node randomly selected nodes permanently locked ON or OFF (this too is random). In the second ensemble we remove n link randomly selected links. Boolean gates for the targets of each such link are derived using the - Figure 1B Modules iScience Article procedure used when isolating modules; namely the portion of the original gate where the removed input is locked in its least canalizing value is preserved. For the third ensemble, we choose n gate random nodes, and we flip a single random output value of their response logic table. This alters their response for a single, specific combination of input values.
The following virtual experiment line added to Virtual_Experiment_List.txt generates an expression/activity time course for all three ensembles (1000 mutant networks each) starting from an initial state matching that of attractor ID 31, responding to 100 time-steps of high growth factor exposure (followed by its reversal), with 15 random errors/network in each ensemble.

Modeling non-saturating inputs and partial knockdown/overexpression within a Boolean framework
To test the effect of non-saturating, intermediate environmental inputs we simulated the model while overriding one or more input node in each time-step with a stochastically generated ON/OFF value of a set probability, thus tuning the cell's average ''signal exposure'' between 0 (none) to 1 (saturating). Similarly, we model partial or complete knockdown/hyperactivation of a node by forcing it OFF/ON with a set probability in each time-step, otherwise leaving it to obey the Boolean rule that normally governs its behaviors. 45,80 This offers a way to model partial inhibition/hyperactivation seen in cells or tissues treated with a drug (e.g., a chemical inhibitor), siRNA silencing, or a pool of a constitutively active proteins alongside the endogenous one. All these result in partial inhibition/saturation of the molecule, leading to blunted downstream effects. In contrast, full knockdown only mimics gene-edited null cell lines where both gene copies (or all homologues represented by a single Boolean node) were removed. Furthermore, as the intermediate levels we model represent cellular concentrations at which the targets our knocked down molecule influences are near the inflection point of S-shaped response curves (i.e., near their activation threshold), the effects of intermediate activity in real cells become highly noisy. 77,78 Our method of randomly allowing a node to affect its targets or be locked OFF/ON mimics an intermittent signal that can stochastically activate its downstream cascade, but the signal flowing through the cascade is unreliable. This is a coarse method to model the coupled effects of losing some but not all signal propagation and the resulting noisy response.
The following Virtual_Experiment_List.txt lines generate a stochastic time course of the model's dynamics in non-saturating environments, with or without partial node knockdown/hyperactivation (ReganLabBoo-leanSims expects a complete list of input node settings, as shown below): Boolean network modules representing distinct cellular regulatory functions Overview of previously published regulatory modules 1. Growth factor signaling (green node clusters on Figure 1C, details in 80 ) is dynamic module tracking basal survival signaling through PI3K/AKT, the mitogen-induced MAPK cascade (lime green), cell cycle linked oscillations in PI3K/AKT signaling (bright green), and mTORC1 activation (mustard).

Automated phenotype detection
The code package ReganLabBooleanSims uses an automated process to assign specific biological phenotypes to any state or time-course of the model and monitor phenotypic transitions 45 ; namely healthy vs. erroneous cell cycle progression, apoptosis, as well as contact inhibited, migratory, epithelial, mesenchymal and hybrid E/M states. The signatures of these biological states are user-defined as follows and read by ReganLabBooleanSims from a file named VirtualExp_Phenotypes_for_STATS.txt (Data S1/ File S10).

Phenotype statistics in non-saturating stochastic environments
To generate phenotype statistics in response to non-saturating environmental inputs, we ran time courses of T stats = 50,000 time-steps in which each input node was randomly toggled ON/OFF in each time-step with an input-specific ON-probability. The simulations tracked the number of cell cycle events, marking cycles completed without error ( Figure S1B, blue bars), the number of genome duplication even iScience Article from G2 (orange bars), the number of premature metaphase-anaphase transitions that did not involve completion of the mitotic spindle followed by genome duplication (green bars), and the number of genome duplication events in the absence of a cytokinesis step between telophase and the next S-phase (red bars). In addition, the simulation counted and stopped at each apoptotic event. Time courses that resulted in apoptosis before time T stats were restarted until a minimum of T stats steps of live-cell dynamics were sampled. In addition, the simulation tracked the average time spent in each reversible phenotype listed above by counting time-steps in which the model's state satisfies the constraints that define each phenotype of interest. To generate phenotype statistics with incomplete knockdown or overexpression of target molecule (s), we combined the non-saturating stochastic inputs described above with a similar stochastic locking of the target molecule OFF or ON with a tunable probability p KD (knockdown) or p OE (over-expression), respectively. In time-steps where the molecule was not locked ON or OFF, it followed the internal Boolean regulatory influences of the rest of the network as if it was unperturbed.
All phenotype statistics are generated by ReganLabBooleanSims via in silico experiment(s) in VirtualExp_ Phenotypes_for_STATS.txt, as illustrated below. Once the statistics are completed and written to file, ReganLabBooleanSims auto-generates the Python code required to show the results as stacked bar graphs such as Figures S1B-S1D. The simplest of these runs generates results as a function of a single, increasing input variable: Revealing the effects of a partial knock-out often requires choosing the right environment where the effect can manifest. To this end, adding an additional (J_Ecadherin KD) parenthesis at the end of NonSatura-ting_Stats_Scan_1Env_fnKDOE commands generates a different output. Namely, ReganLabBooleanSims generates a grid of bar graphs, all as a function of increasing knockdown. The horizontal direction along the grid shows results for increasing input settings (e.g., GF_High exposures from 5% to 95% below). The assayed phenotypes are organized vertically.
Finally, ReganLabBooleanSims can generate results two increasing, non-saturating environmental inputs (e.g., Stiff_ECM and CellDensity_Low) as a 2D grid of bar graphs (1 page for each phenotype). The x axis within each bar graph corresponds to increasing KD/OE of the node in the last parenthesis (e.g., - Figure S1B; code generates a figure panel with bar graphs of all phenotype statistics at increasing GF_High from 0% to 100% (x axis).