Parity and time-reversal elucidate both decision-making in empirical models and attractor scaling in critical Boolean networks
Jordan C. Rozum, Jorge Gómez Tejeda Zañudo, Xiao Gan, Dávid Deritei, Réka Albert
Parity and time-reversal elucidate decisions in high-dimensional state space - application to attractor scaling in critical Boolean networks
Jordan C. Rozum , Jorge Gómez Tejeda Zañudo , Xiao Gan , Réka Albert Department of Physics, The Pennsylvania State University, University Park, PA 16802, USA Eli and Edythe L. Broad Institute of MIT and Harvard, Cambridge, MA, 02142, USA Department of Medical Oncology, Dana-Farber Cancer Institute, Harvard Medical School, Boston, MA, 02115, USA Network Science Institute and Department of Physics, Northeastern University, Boston, MA 02115, USA Channing Division of Network Medicine, Department of Medicine, Brigham and Women’s Hospital, Harvard Medical School, Boston, MA 02115, USA Department of Biology, The Pennsylvania State University, University Park, PA 16802, USA
Abstract
We present new applications of parity inversion and time-reversal to study the emergence of complex behavior from simple dynamical rules. Consideration of emergence has become increasingly important in engineering, medicine, and social science. Here, our focus is a common class of discrete, nonlinear, and stochastic models that is often applied in these fields. We describe a new parity-based encoding of causal relationships in such systems that facilitates several new theorems regarding the controllability of complex systems and their behavior under widely used state-space projection techniques. Together with a new time-reversal construction, this encoding (called the “parity expanded network”) allows one to quickly identify discrete analogs of stable and unstable manifolds, as we illustrate by studying decision-making in models from the systems biology and statistical physics literature. These applications underpin a novel attractor identification algorithm implemented for Boolean networks under stochastic dynamics. The algorithm’s performance advantages over existing methods have enabled the resolution of a longstanding open question of how the number of attractors scales with network size, and in particular whether the scaling matches what is observed biologically. Our results are based on the exact (or near-exact) enumeration of attractors in a Boolean network, and not on numerical simulations or on the sampling of state space. We consider critical random Boolean networks more than 80 times larger than ever before studied under stochastic update (N=16,384). We provide strong evidence of a power law scaling with exponent 0.12±0.05 (95% CI [0.04,0.22]), a quarter of the previously conjectured value, and an order of magnitude smaller than the analytically derived upper bound. Crucially, the scaling is also much slower than what has been derived from experimental data. This suggests that additional mechanisms are at play that enable organisms to exhibit their comparatively diverse array of phenotypes.
Popular Summary
Emergence is the phenomenon of complex behaviors arising from the collective behavior of simple parts. Understanding emergence is central to many problems across scientific disciplines, including those involving phase transitions, cell differentiation, or opinion spreading. Boolean networks are a common framework for studying emergence that integrate agent-based modeling, nonlinear dynamics, and network analysis. Their statistical properties have been of interest in statistical mechanics, while their value as an empirical tool is often applied in biology. Here we introduce a new approach to finding the long-term behaviors (attractors) of these systems. We adapt classic methods in mechanics to define new parity and time reversal constructions and describe their relationship with existing state-space projection methods. We apply our methods to random Boolean networks, standard models of the genetic interactions that determine a cell’s function. Our methods are successful in stochastically evolving networks approximately 80 times larger than ever before considered. We discover that the average number of attractors increases very slowly with the size of the network -- the scaling exponent is ten times smaller than the previous best upper bound, and less than one fifth the value estimated experimentally. This suggests that regulatory networks sustain many more stable phenotypes than expected by chance in standard models, inviting investigation into the biology underlying this surprising phenotypic diversity. The primary insights of this work rely on the properties of a Boolean system’s “expanded network”, which has generalizations for multi-level and ODE systems. Thus, we demonstrate a general principle: a system’s relationship to its time reversal and state-space inversion constrains its repertoire of emergent behaviors.
Introduction
Many complex systems in the natural, social or technological realm exhibit emergent behavior, i.e., collective dynamics arising from the interaction of entities governed by simple rules [1–5]. Examples include phase transitions [6,7], flocking [8], consensus formation [9], and spontaneous synchronization of oscillators [10]. Modeling frameworks that are frequently used to study the collective behavior of individuals include nonlinear dynamics [11,12], agent-based models [13], cellular automata [14], and network models [5,15,16]. Boolean models sit at the intersection of these approaches. They assign a time-varying binary variable to each system entity, represented as a node in a network of interactions. They exhibit diverse long-term dynamics (attractors) that represent collective behavior (for example, consensus of individuals) and they describe the evolution toward an attractor (for example, consensus formation from an initially disordered state). Boolean modeling of electronic circuits is well-known among physicists [17], but many other familiar models can also be viewed as Boolean models. The quenched (zero-temperature) Glauber model, a dynamic variant of the Ising model, considers the dynamics of each atom’s two possible spin orientations under the influence of its neighbors [18,19]. Another example category, widely studied in statistical mechanics, includes models of spreading binary opinions through social networks (reviewed in [5,9]). The McCulloch & Pitts neural network model introduces a propositional logic of all-or-none neuronal activation [20]; the Hopfield model also considers two activities for each neuron and assumes a complete network of interactions [21]. Boolean models are especially well suited to elucidate system-level decision-making, i.e., robust commitment towards one of the dynamical attractors in a multi-stable system. This has made their use especially wide-spread in biology. They were introduced by Stuart Kauffman and René Thomas [22,23] as prototypical models for gene regulatory networks that underlie cell fate decisions (such as those that happen during cell differentiation). A large body of research has since shown that the attractors of Boolean models correspond to cell fates or stable patterns of cell activity (such as the cell cycle). Boolean models have become a common tool in systems biology due to their ability to integrate and encode current knowledge of a biological process, fill any gaps of knowledge with hypothesized interactions, and predict the behavior of the system under untested circumstances, for example novel perturbations [24]. They are frequently used to study cell differentiation processes (e.g. T cell specialization [25,26]), developmental processes (e.g. patterning during embryogenesis in
Drosophila melanogaster [27–29]), and cancer (e.g. metastatic reprogramming [30–32], prediction of targeted therapies [33,34]). Model predictions in a variety of systems were verified experimentally [32,35–42]. Alongside models of specific systems, analysis of the expected collective behaviors exhibited by generic Boolean models has also proven insightful. Ensembles of Boolean models (Random Boolean Networks) have been studied by physicists for decades (reviewed in [43–45]). These ensembles exhibit many features of biological cells, including stability against perturbations and plausible scaling laws for the number and size of attractors with the system size. Despite Boolean models’ discrete nature and apparent simplicity, it is nontrivial to connect dynamical properties of decision-making to the underlying interaction network, and brute-force exploration of their state spaces is not generally feasible. A typical Boolean model of a biological process with a few dozen variables has tens of billions of states. Genome-scale models can have thousands of variables, resulting in too many states to encode in the observable universe. This challenge has motivated decades of research analyzing Boolean dynamics without exhaustive state-space searches [46]. One approach is to analyze how interaction topology constrains dynamics. René Thomas and collaborators demonstrated that a positive feedback loop is necessary for the emergence of multistability and a negative feedback loop is necessary for sustained oscillation [23,47]. These conditions apply not only in Boolean systems (reviewed in [48]) but also in more general discrete dynamical systems as well as continuous systems (reviewed in [49]). While the body of research regarding how network structure constrains dynamics has proven invaluable, it is important to note that multiple Boolean systems are compatible with each interaction network. Ambiguity can be eliminated by defining a network whose graph structure unambiguously represents the update functions that govern the time evolution of each variable. One such network representation, the expanded network (also called the logical or prime-implicant hypergraph) [28,50–52], defines two virtual nodes for each entity, denoting the two possible values of its binary variable. Connections among virtual nodes encode the update functions. The structure of these auxiliary networks has been used to help identify the attractors of a system [31,35,53,54] and to determine control strategies to drive the system to a desired attractor [36,55]. The most parsimonious of these control strategies involve determining the so-called driver node(s) of self-sustaining circuits [56,57]. In this way, fixing the state of a few driver nodes ensures the system's convergence into a target space from any initial condition. Expanded network structures have also proven useful in studying the behavior of complex oscillations in Boolean models [54]. Although initially developed in the context of Boolean modeling, the underlying concept of the expanded network is not a purely Boolean phenomenon and has been generalized to multi-level discrete and ODE analogs [58,59] where it has been used to study the robustness and parameterization of feedback loops [60,61]. Here, we go beyond the previous use of the expanded network and characterize each virtual node with a binary activity and impose a parity structure that facilitates the proof of new theorems about the effects of perturbations on system trajectories. With these additions, the expanded network is a complete representation of the Boolean dynamical system. In addition, we describe for the first time the time reversal of stochastically asynchronous Boolean systems. Using parity and time reversal transformations in tandem, we developed a new algorithm to efficiently identify all attractors of large-scale Boolean systems. We apply the algorithm to answer the long-standing question of how quickly the number of attractors in asynchronous Random Boolean Networks increases with network size.
Overview of Boolean modeling
Constructing a Boolean model usually starts with the synthesis of the modeled system’s interaction graph. An interaction graph is a signed directed graph whose nodes are the entities of a system. Positive edges indicate an activating influence, while negative edges indicate an inhibitory influence (unsigned edges are also possible). In a Boolean model, the activity of each entity 𝑖 is characterized by a variable 𝑋 that can take one of two values: 1 (“active”) or 0 (“inactive”). Each 𝑋 updates its value according to the output of an update function 𝑓 which maps every system state 𝑋 to either or . The entities upon whose variables 𝑓 depends are the regulators of 𝑖 (as encoded in the interaction network). Any Boolean update function 𝑓 can be specified explicitly as a truth table or algebraically using logical operators (e.g., “ 𝑛𝑜𝑡 ”, “ 𝑜𝑟 ”, “ 𝑎𝑛𝑑 ” with symbolic representations ¬ , ∨ , and ∧ , respectively). One useful representation is the disjunction of the function’s prime implicants (irreducible sets of regulator states that result in 𝑓 (𝑋) = 1 ) [62]. Alternative representations, such as Boolean threshold functions, are able to efficiently express special classes of Boolean functions, but are not generally applicable [63–66]. There are several schemes for determining the timing of variable updates. In the traditionally used synchronous update all variables are updated simultaneously at each time-step, making the system’s dynamics deterministic. The assumption of synchronicity is not suitable for systems that exhibit multiple time scales or stochasticity [67,68]. Instead, we apply stochastic asynchronous update , in which at each step a single variable is randomly chosen to update its value (each variable must have a non-zero update probability). This update scheme contains the frequently studied uniform update probability as a special case. Updating asynchronously removes spurious oscillations that arise from unrealistic perfect synchrony [67,69]. Once the update functions are determined and an update scheme is selected, the Boolean system is fully specified. Each Boolean system induces a state transition graph (STG) with nodes that represent all possible system states and with directed edges from one node (system state) to another when the parent state can be updated in one time-step to attain the child state. For a Boolean system with 𝑁 entities the STG has nodes. Under stochastic asynchronous update, there are between and 𝑁 edges leaving any node of the STG. Importantly, the STG is the same for any choice of nonzero update probability for each node. The STG can have as many as 𝑁 edges (i.e., for 𝑓 (𝑋) = ¬𝑋 , 𝑖 = 1. . 𝑁 ). Whenever a state is such that the output of a variable’s update function leaves that variable unchanged, that state has a self-loop in the STG; these are often left implicit (hidden). Source nodes in the STG are states that are only attainable as initial conditions and are termed “Garden of Eden states” [70,71] . The attractors of a Boolean system are the terminal strongly connected components of the STG (i.e., they have no edges that exit the component). They are divided into two types: point attractors (also called fixed points or steady states), which contain only one state, and complex attractors, which contain more than one state. Informally, attractors of a Boolean system represent the long-term behaviors of the system, much like they do in ODE systems. Point attractors and the simplest type of complex attractor (i.e., a cycle in which each state transition involves the change of a single variable) are independent of the timing of variable updates [69,72]. The Kauffman 𝑁 − 𝐾 random Boolean network (RBN) model [22] has been frequently studied in systems biology and statistical physics for over 50 years. In this model, each of 𝑁 nodes receives 𝐾 edges from randomly selected nodes, and the (quenched) update functions are chosen such that each of input combinations yields an output of with probability 𝑝 . These models traditionally use synchronous update. Their attractors exhibit many features of biological cells, including stability against random external perturbations and plausible scaling laws for the number of attractors and attractor cycle lengths with the number of nodes 𝑁 [45,73]. Specifically, tuning the indegree 𝐾 or the activation bias 𝑝 can produce an order-to-chaos transition at (in the thermodynamic limit, when 𝑁 → ∞ ) [43,74]. In the ordered regime, almost all nodes quickly attain a stationary state and small, transient perturbations tend to dissipate. In the chaotic regime, on the other hand, the number of fluctuating nodes is proportional to 𝑁 and perturbations grow exponentially fast. At there is a so-called critical regime in which, on average, small perturbations persist indefinitely but remain small. RBNs with more heterogeneous topologies [75] or alternative distributions of update functions [65,66,76–78] also exhibit these regimes. The networks of greatest interest for biological systems were conjectured to lie near the critical regime [73]. More recent work analyzing the dynamical behavior of gene regulatory networks [79–82] and studying evolutionary models of gene regulatory networks [83,84] has provided further support for this conjecture. There has been a rich history of research on the scaling of the average number of attractors with network size in the biologically relevant critical regime (see [45] for a review of this work). Based on the relationship between the number of cell types and DNA content, Kauffman originally conjectured that the average number of attractors grows with the number of genes as a power law with exponent [22,73]. In the half-century since the introduction of the 𝑁 − 𝐾 ensemble, numerical and theoretical analyses by multiple research groups have shown that the scaling of the number of attractors depends on the updating scheme and is not necessarily a power law [85–91]. Under synchronous update, the average number of attractors for the
𝐾 = 2, 𝑝 = 0.5 critical Kauffman networks grows faster than any power law (in contrast to earlier conjectures) [90]. Under stochastic asynchronous update, the number of attractors grows as a power law, with an upper bound given by 𝑁 @AB ≈ 𝑁 D.EF [87]. The existence of a power law scaling in the asynchronous case is supported by a study of stable synchronous attractors (which are preserved under a perturbation of node update synchrony) in networks of up to
𝑁 = 40 nodes [67]. This study reported an exponent of and a low average number of stable synchronous attractors ( for
𝑁 = 30 ). Another relevant result is that in
𝑁 − 𝐾 ensembles with 𝑝 = 0.5 the average number of point attractors is 1, regardless of the value of 𝑁 or 𝐾 (implying that for large 𝑁 at least a small number of nodes oscillate in most attractors) [92]. Importantly, for the 𝐾 = 2 critical Kauffman networks under stochastic asynchronous update, the exact value of the power law exponent is still unknown.
A new framework: The expanded network through the lens of parity
We use the parity transformation of a Boolean system to encode the update functions as a network whose structure provides insights into the system’s dynamics. The parity transformation acts on a Boolean system by the change of variables 𝑋 ↦ ¬𝑋 for all variables 𝑋 , resulting in the inversion of the transformed system. Under this change of variables, the new update functions are the negations of their original counterparts. Formally, we define the parity transformation by its action on a Boolean system's update functions, i.e., as mapping the original system with update functions 𝑓 to the parity transformed system governed by update functions ¬𝑓 . As we shall see, this induces further transformations on any structure derived from the update rules (e.g., the state transition graph), and so, in a slight abuse of notation, we speak of the parity transformation as acting on all these structures simultaneously. The Boolean parity transformation is analogous to the spatial inversion of a physical system in which all coordinate variables are mapped to their new values by negation. This analogy is especially appropriate when considering its effect on the state transition graph of the system. The parity transformation corresponds to a relabeling of the nodes of the state transition graph so that all 1s become 0s and vice versa (see Figure 1). One may view the node labels in the state transition graph as spatial coordinates (so that states lie on the vertices of a unit hypercube, with the all-off state at the origin and the all-on state at the opposite corner). In this layout, the parity transformation of the Boolean system corresponds exactly to a spatial inversion of this hypercube through its center (see Figure 1). A new definition of the expanded network using parity-inversion provides a complete and concise structural description of system dynamics
The expanded network (also called the logical or prime-implicant hypergraph) [28,50–52] was introduced as an auxiliary network constructed from the Boolean update functions. The expanded network nodes represent Boolean literals (e.g., 𝑋 and ¬𝑋 ) and its hyperedges (generalized edges that connect sets of nodes) represent prime implicants of the update functions. Here we introduce a new definition of the expanded network that uses parity-related concepts and highlights its role as an invariant of the parity transformation. This new definition allows one to study the expanded network as a dynamical system in its own right and greatly facilitates our formal proofs by treating the active and inactive states on equal footing. Each node 𝐼 of the expanded network, called a virtual node , is an ordered pair 𝐼 = (𝑛(𝐼), 𝑠(𝐼)) consisting of a system entity, in this context denoted 𝑛(𝐼) , and a value 𝑠(𝐼) , which is either the constant or the constant . In a Boolean model there are two virtual nodes associated with each system entity 𝑖 , namely (𝑖, 1) and (𝑖, 0) ; we call this pair of virtual nodes contradictory. A set of virtual nodes that does not contain any contradictory pair is called consistent . The activity variable for 𝐼 , written 𝜎 O , is the Boolean literal 𝑋 A(O) if 𝑠(𝐼) = 1 or ¬𝑋 A(O) if 𝑠(𝐼) = 0 ( 𝜎 O can be written agnostically as ¬𝑋 A(O) 𝑥𝑜𝑟 𝑠(𝐼) ). The value of 𝜎 O in a state 𝑋 is written 𝜎 O (𝑋) . For example, if we consider a virtual node 𝐼 = (𝑖, 0) , then 𝜎 O (𝑋) = 1 holds when 𝑋 is in state 𝑋 . A set of virtual nodes 𝑆 is active (in a system state 𝑋 ) if all its members are active (i.e., if 𝜎 O (𝑋) = 1 for all 𝐼 ∈ 𝑆 ). The update function 𝐹 O for each virtual node’s activity is inherited from the update function for 𝑛(𝐼) , i.e., 𝐹 ( (𝑋) =𝑓 (𝑋) and 𝐹 ( (𝑋) = ¬𝑓 (𝑋) are the update functions for (𝑖, 1) and (𝑖, 0) ; note that the activity of contradictory virtual nodes must always be updated together and have opposite states. A hyperedge connects a set of parent virtual nodes 𝑆 = {𝐼 T , 𝐼 D , . . . , 𝐼 V } to a target virtual node 𝐽 if ∧ O∈Y 𝜎 O is a prime implicant of 𝐹 Z . Pictorially, we represent hyperedges with more than one parent using intermediary “composite nodes”, which correspond to “ 𝑎𝑛𝑑 ” gates. For an example of a Boolean system and its expanded network, see Figure 1B,C. Figure 1: The relationship between the expanded network and parity transformation illustrated on a three-variable Boolean system. The state transition graphs are shown in panel A; each system state is represented as the triple 𝑋 [ 𝑋 \ 𝑋 ] . Each state with fewer than three outgoing edges (state transitions) also has a self-loop, which is omitted for visual clarity. The interaction network and update functions are indicated in panel B and the expanded network in panel C. The system has two attractors, shown in panels D and E, with blue nodes active and grey nodes inactive. These two attractors (D and E) are indicated in green and purple in panel A, respectively. States in the state transition graphs of panel A are arranged so that the individual variables define coordinate axes. In this arrangement, the states form the corners of a cube. The parity transformation reflects each attractor through the center of this cube. The expanded network in panel C has two parts, or layers: regular virtual nodes (A,1), (B,1) and (C,1), whose activity updates according to the usual update functions, and negated virtual nodes (A,0), (B,0) and (C,0) that use the parity transformed update functions. The filled black circle represents the hyperedge from the set {(A,0), (B,0)} to (C,0) and indicates the “ 𝑎𝑛𝑑 ” operation in the update function of the virtual node (C,0). Positive regulation (black arrows) stays within a layer of the expanded network, while negative regulation (red arrows) crosses between layers. In any state of the system, half of the expanded network is active. For example, attractor 1 (panel D) has the virtual nodes (A,1), (C,1), (B,0) active (indicated in blue). Attractor 2 (panel E) has the virtual nodes (A,0), (B,1), (C,1) active. The active nodes of any system state (A,C in attractor 1) are the active regular nodes of the expanded network, and the active negated nodes correspond to the nodes that are active in the corresponding state of the parity transformed system (B in attractor 1). The above definitions recapture the expanded network as defined in [52] by combining information from the original Boolean system with that of its inversion to represent the update functions in a purely structural form. We can highlight the role of the parity transformation by partitioning the nodes of the expanded network according to whether each virtual node 𝐼 is “regular” (i.e., 𝑠(𝐼) = 1 ) or “negated” (i.e., 𝑠(𝐼) = 0 ). The activity update functions 𝐹 O for the regular nodes are given by the update functions 𝑓 of the original system, whereas the activity update functions for the negated nodes are given by the parity-transformed update functions. We can thus view the expanded network as having two layers: one with activities that update according to the functions of the original system, and one with activities that update according to the functions of the inversion of the Boolean system (see Figure 1C). The parity transformation interchanges these layers, but leaves the structure of the expanded network invariant. Because each Boolean literal is represented in one layer or the other, Boolean negation is no longer needed to describe the update functions; negative influence instead manifests as inter-layer hyperedges. In this view, inhibition is strongly related to the interplay between a Boolean system and its own inversion. Indeed, one can show that if a Boolean system lacks negative feedback loops and has no paths of opposite sign between any two nodes, then there is a change of variables that leaves the parity layers of the expanded network disconnected from one another; this follows from Theorem 19 of [48] and the observation that inter-layer hyperedges in the expanded network can only arise from inhibition. This change of variables is achieved by inverting a proper subset of the variables. Stable motifs are subgraphs of the expanded network that describe trap spaces
The expanded network transforms many of the dynamical properties of a Boolean system into purely structural features. Arguably the most important of these structural features are the stable motifs [52], which correspond to specific states of generalized positive feedback loops. Stable motifs belong to the broader class of stable modules [59], which determine trap spaces in the dynamics -- regions of the state-space that, once entered by a trajectory, cannot be exited (a system’s attractors are a notable special type of trap space). Stable motifs were originally defined for expanded networks viewed as bipartite graphs, rather than for hypergraphs. We recover this previously defined concept in the parity-transformation view of the expanded network by way of an alternative definition. A stable module 𝑀 is a non-empty sourceless sub-hypergraph of the expanded network such that 𝑀 does not overlap (does not share virtual nodes) with its image under parity. A stable motif is a stable module that does not contain any smaller stable module; note that this implies that a stable motif is strongly connected. For example, in Figure 1 the virtual nodes (A,1) and (B,0) form a stable motif, as do the virtual nodes (A,0) and (B,1). Because stable motifs correspond to positive feedback loops, they necessarily have an even number of inter-parity layer edges in the expanded network (this implies that for any given stable motif, there is a change of variables for which all virtual nodes of that stable motif correspond to “ON” states). Because stable modules (and thus also stable motifs) are sourceless, every virtual node in a stable motif 𝑀 can be maintained in its active state by other virtual nodes in 𝑀 , meaning that the activity of 𝑀 is self-sustaining. To see this precisely, recall that a set of virtual nodes is active if all its members are active. As such, a stable motif 𝑀 is active in a system state 𝑋 if each variable to which a virtual node of the stable motif refers takes the value described by the virtual node. Once 𝑀 is activated, it cannot be inactivated except via direct override of its virtual node activities by external controls (as opposed to inactivation by external control of upstream pathways). Thus, 𝑀 describes a control-robust trap space in which the values of certain variables are stationary [50,55]. The trap spaces corresponding to the activity of stable modules are exactly those considered by [53], with larger trap spaces (more states) corresponding to smaller stable modules (fewer constrained variables) and vice-versa. One can verify in Figure 1 that the stable motif {(𝐴, 1), (𝐵, 0)} (which describes the trap space containing attractor 1) does not overlap with its parity transformed set {(𝐴, 0), (𝐵, 1)} . In this example, the parity transformed set also forms a stable motif (which describes the trap space containing attractor 2), but in general this set would only be a conditionally stable motif (as defined in [54]. New theorems regarding driver nodes in the expanded network
Boolean models offer a natural representation of large perturbations. For example, keeping a variable in a fixed state simulates a loss-of-function or constitutive activation of the respective entity [23]. This type of analysis has been used, for example, to pinpoint key proteins in pathological cell processes [55,56,58]. A related problem of great interest is the identification of interventions that ensure convergence to a desired attractor [93,94]. In this section we translate important results on the properties of driver nodes (interventions that confine the system's trajectories [57]) into the new parity-inversion definition of the expanded network. In addition, we leverage the parity properties of the expanded network to prove new results about driver node sets, and their relation to attractors. In particular, we show that if a virtual node,
𝐼 = (𝑖, 𝑠) , leads to activation of a stable motif, then any attractor that does not have that stable motif active also must have 𝐼 inactive, and so 𝑋 cannot oscillate in that attractor. Recall that a pair of nodes 𝑋, 𝑌 is called consistent if they are not contradictory (we define (𝑖, 1) and (𝑖, 0) as contradictory). We say that a set of mutually consistent (non-contradictory) virtual nodes 𝑆 drives a virtual node 𝑌 ∉𝑆 if fixing all virtual nodes
𝑋 ∈ 𝑆 to be active yields a controlled system for which all attractors have 𝑌 active. We say that 𝑆 drives 𝑋 ∈ 𝑆 if 𝑆 drives all the parents of a hyperedge targeting 𝑋 . The set of all virtual nodes driven by 𝑆 is called the domain of influence (DOI) of 𝑆 , written 𝐷𝑂𝐼(𝑆) [57]. Calculating
𝐷𝑂𝐼(𝑆) can be difficult in general, so we focus instead on a commonly used and easily calculated subset of the driving relation: the logical domain of influence (LDOI). A set 𝑆 of virtual nodes logically drives a virtual node 𝑌 (which may or may not be in 𝑆 ) if there is a path in the expanded network hypergraph from a subset of 𝑆 to 𝑌 with all virtual nodes in the path consistent with 𝑆 (this requires that 𝑌 is consistent with 𝑆 ). The set of all 𝑌 logically driven by 𝑆 is called the logical domain of influence of 𝑆 , written 𝐿𝐷𝑂𝐼(𝑆) . When considering the LDOI or DOI of a set
𝑆 = {(𝑖, 𝑠)} of size one, we write
𝐿𝐷𝑂𝐼(𝑖, 𝑠) or 𝐷𝑂𝐼(𝑖, 𝑠) to avoid excessive notation. Informally,
𝐿𝐷𝑂𝐼(𝑆) corresponds to the variable values that become fixed after percolating 𝑆 through the update functions and simplifying algebraically. As demonstrated in [57], if 𝑆 is a subset of 𝐿𝐷𝑂𝐼(𝑆) , then
𝐿𝐷𝑂𝐼(𝑆) contains a stable motif. If
𝐿𝐷𝑂𝐼(𝑆) contains a stable motif 𝑀 , we say that 𝑆 (logically) drives 𝑀 . Here we make several new connections between driver sets and the attractor repertoire of a system. Formal statements and proofs are given in Appendix A. Our main result in this section states that if 𝑆 is active in some part of an attractor, then 𝐷𝑂𝐼(𝑆) is also active somewhere in that same part (Theorem 1, Appendix A). This result illustrates how the domain of influence of a set of virtual nodes can constrain which states must coexist within an attractor. Such considerations are important in constructing Boolean models in which certain system configurations should correspond to different long-term qualitative system behaviors (e.g., phenotypes). Two corollaries of this result allow one to study the conditions under which an attractor avoids activating a particular set of stable motifs. This problem is of interest, for example, in biology, where one or more stable motifs may correspond to a diseased state of the system; the goal in this case is to identify drug targets that avoid stabilizing the diseased state. We will also make use of these results in later sections to enumerate a system’s attractors. The two corollaries of Theorem 1 require us to introduce the concept of contradiction boundary . If a virtual node 𝑌 is not consistent with a set of virtual nodes 𝑆 , but there is a path from a subset of 𝑆 to 𝑌 with all path nodes except 𝑌 consistent with 𝑆 , then we say that 𝑌 is in the contradiction boundary of 𝑆 . In this case the LDOI of 𝑆 will contain a parent node of 𝑌 , but it will not contain 𝑌 . If 𝑆 is internally inconsistent (i.e., if it contains contradictory virtual nodes), then we say that the contradiction boundary of 𝑆 is the set of contradictory virtual node pairs in 𝑆 . The contradiction boundary captures whether a set of virtual nodes 𝑆 will eventually activate a state that contradicts 𝑆 . If the contradiction boundary of 𝑆 is non-empty, we say that 𝑆 is self-negating . The concept of contradiction boundary is important because of the following property: if 𝑆 is self-negating, then there is no attractor of the system in which all of 𝑆 is permanently active. The first of the two corollaries (Corollary 1, Appendix A) can be viewed as a compatibility condition for a stable motif’s activity and the activity of a virtual node that drives it. In the context of system control, it states that if fixing 𝑋 is sufficient to eventually activate 𝑀 , then the oscillation of 𝑋 is also sufficient to activate 𝑀 . It can also be viewed as a consistency condition for attractors: if 𝑋 = 𝑠 leads to activation of 𝑀 , then we cannot have 𝑋 = 𝑠 in an attractor (even transiently) in which 𝑀 is inactive. This result provides a powerful way to identify circumstances in which no stable motif activates. Specifically, we collect all single-node drivers of all stable motifs into a set 𝛥 , and test whether or not ¬𝛥 = {¬𝐼: 𝐼 ∈ 𝛥} is self-negating (whether or not its contradiction boundary is empty); here ¬𝐼 = ¬(𝑖, 𝑠) =(𝑖, ¬𝑠) . If it is self-negating, then it follows that at least one element of ¬𝛥 is not permanently active in each attractor, and thus (by Corollary 1) at least one stable motif must eventually stabilize. The formal statement of this result is given as Corollary 2 in Appendix A. Time-reversal of Boolean systems and stable motifs of the time-reversed system as unstable “Gardens of Eden”
A second set of foundational results we introduce in this work is the construction of the time reversal of an asynchronous-update Boolean system -- despite the system’s inherent stochasticity. We use the time reversal transformation to help identify discrete analogs of unstable manifolds in the state space. The time reversal of a Boolean system governed by update functions 𝑓 is called the time reversed system and is governed by the update functions 𝑓 , where 𝑓 = ¬𝑓 (𝑋 = ¬𝑋 ) . Like the parity transformation, we formally define the time reversal as a transformation of the Boolean system given by its effect on the update functions. It similarly induces transformations of structures that are derived from these functions. For example, the state transition graph of a system is related to its time reversed counterpart by reversing the direction of all edges (see Figure 2A). Thus, one may follow the evolution of the time reversed system on the STG by following edges in the reverse direction. Equivalently, and in analogy to concepts in solid-state physics, one may imagine that all but one of the nodes in the STG are occupied by walkers who may not share an STG node. If the walkers randomly follow the edges of the STG, then the position of the unoccupied “hole” evolves according to the possible trajectories of the time reversed system Figure 2: Illustration of the time reversal transformation on the example from Figure 1. Panel A shows the original system’s STG alongside the STG of its time reversal. The time reversal has the effect of reversing the direction of each edge in the STG. The interaction network and update functions of the system are depicted in panel B alongside its expanded network in panel C. Red interaction network edges indicate inhibition, while black edges indicate activation. The virtual nodes and hyperedges of the system’s two stable motifs are highlighted in blue and green in the expanded network, and the regions of state space in which they are active are highlighted in the system’s STG on the left in panel A. The blue region includes the point attractor 011 and the green region includes the point attractor 101. The same information is presented in panels D and E for the time reversed system; here, stable motifs and their corresponding nodes in the time reversed STG are highlighted in orange and yellow. Note that the system states highlighted in orange (000 and 001) viewed in the original system’s STG form a subgraph that has no in-component and the state 000 is a Garden of Eden state. The same is true for the system states highlighted in yellow. In this example, the sign of the regulation (inhibition vs activation) is reversed under the time reversal. This holds in general for all interactions except self-regulation, which does not change sign under time reversal. The concept of Garden of Eden states, which are source nodes of the STG [71], can be generalized as subgraphs of the STG that do not have incoming edges; we call these subgraphs Garden of Eden spaces. They are analogous to unstable manifolds of ODE systems in the sense that no trajectory can enter these spaces from the outside. Any trap space of a system, and in particular any of its stable motifs, is a Garden of Eden space in the time reversed system, and vice versa. For example, the states marked in blue in Figure 2A form a trap space of the original system and a Garden of Eden space of its time-reversed counterpart while the states marked in yellow form a trap space of the time-reversed system and a Garden of Eden space of the original system. An important consequence of this time-reversal-based mapping between trap spaces and Garden of Eden spaces is that no attractor of a system can cross the boundary of any of the system's trap spaces or Garden of Eden spaces. This observation is especially helpful in eliminating states from consideration when searching for attractors via direct STG construction, or in reducing the number of relevant initial conditions for study, as we will see in the next section. Stable motif succession determines state-space decision-making and attractors
As a Boolean system’s state transition graph (STG) has nodes and up to 𝑁 edges (for stochastic asynchronous update), for systems with many entities (large 𝑁 ) it is impractical to use the state transition graph to determine the system’s attractor repertoire. Stable motif, or trap space, based attractor identification methods are often more effective [50,52]. In the iterative approach of [52], a system’s stable motifs are identified, and one is selected to “lock in”. The system is then reduced under the assumption that the system’s state is confined to the region described by the stable motif, resulting in a reduced network . Rephrased in our framework, each stable motif is selected in turn and the expanded network is simplified under the assumption that the motif’s virtual nodes are active, resulting in a reduced Boolean system. We use the notation 𝑅𝑒𝑑(𝑓, 𝑀) for the reduced system that results from a Boolean system with update function 𝑓 after assuming that stable motif 𝑀 is active. The process is repeated for each reduced system until all possible permutations of stable motif activation are explored. The result is a succession diagram, 𝛴 , which is a directed acyclic graph whose nodes are the unions of stable motifs used to obtain each reduced system. Considering two nodes 𝑃 and 𝑄 of the succession diagram 𝛴 , there is an edge from 𝑃 to 𝑄 if there is a stable motif 𝑀 of the reduction 𝑅𝑒𝑑(𝑓, 𝑃) of 𝑓 by 𝑃 such that 𝑄 = 𝑃 ∪ 𝑀 . For example, in Figure 3 the reduction by the blue stable motif (panel D) contains the pink stable motif, thus the succession diagram (panel E) contains an edge from the blue stable motif to the union of the blue and pink stable motifs. Each path in the succession diagram 𝑝 = (𝑝 T = {}, 𝑝 D , . . . , 𝑝 A = 𝑃) from the empty set to a node 𝑃 ∈ 𝛴 defines a stable motif history 𝑀 = 𝑝 \𝑝 for 𝑖 = 0. . 𝑛 − 1 and a corresponding sequence of reduced systems 𝑅𝑒𝑑(𝑓, 𝑝 ) . For each sink node of the succession diagram (i.e., whenever a reduced system has no additional stable motifs) there is guaranteed to be an attractor in which all reduced variables are fixed in their reduced states. If any variables remain unreduced, then one or more of the unreduced variables oscillate in the attractor. By repeating the reduction for every possible choice of stable motifs, one can construct a list of attractors of the Boolean system. For an example of this process, see Figure 3. The succession diagram serves as a summary of the decisions in the system dynamics that lead to successively more restrictive nested trap spaces. Each node of the succession diagram corresponds to a region of the state space in which the denoted stable motifs are active. Each branch point in the succession diagram represents potential choices to be made; which choice is ultimately selected by the system depends on various factors including the stochastic update order. It follows from the nestedness of these trap spaces that the system cannot transition between regions that are not connected by a path in the succession diagram, i.e., the succession diagram encompasses the entire repertoire of decisions the system is capable of making. The close relationship between the succession diagram and branch points in the dynamics is illustrated in Figure 9 (in Appendix B) by constructing the full state transition graph of the example from Figure 3. In the section “Application to decision-making in empirical biological network models” below, we illustrate how the succession diagram can be used to analyze the complex state space decision-making in systems biology models. Figure 3: Outline of the iterative stable motif reduction process on a simple example. Panel A depicts a Boolean system’s interaction network and update functions. The corresponding expanded network is shown in panel B; there are three stable motifs, highlighted in blue, green or purple. In panel C, the stable motif corresponding to 𝑋 [ = 𝑋 \ = 1 is selected and the effect of maintaining its activity is highlighted in blue; in particular, it leads to 𝑋 ] =𝑋 q = 1 . The variables that are unfixed ( 𝑋 r and 𝑋 s ) then form a bistable system whose expanded network is shown in panel D. The reduced system has two stable motifs, highlighted in pink or purple; the latter was also a stable motif of the original system. All possible sequences of stable motif selection and reduction are summarized in a succession diagram (panel E). Each node of the succession diagram is a set of virtual nodes that contains the stable motifs that were selected for use in the reduction. The colors of the edges in the succession diagram indicate which stable motif is selected to get from one reduction to the next. The process terminates when no stable motifs remain, and the attractors of the maximally reduced systems are identified as attractors of the unreduced system (highlighted in yellow, orange, and brown). One must take special care to consider the possibility of oscillations that avoid activation of stable motifs. For example, consider the system shown in Figure 4 𝑓 [ (𝑋) = 𝑓 \ (𝑋) = ¬𝑋 [ ∧ ¬𝑋 \ ∨ 𝑋 ] ; 𝑓 ] (𝑋) = 𝑋 [ ∧ 𝑋 \ . This system contains only one stable motif, {(𝐴, 1), (𝐵, 1), (𝐶, 1)} , which corresponds to the system’s sole point attractor 𝑋 [ = 𝑋 \ = 𝑋 ] = 1 . Previous stable motif or trap-space based methods [52,53] would correctly identify this point attractor by finding the corresponding stable motif. This system, however, contains an additional attractor in stochastic asynchronous update. In the second attractor, 𝑋 ] remains in the 0 state, while 𝑋 [ and 𝑋 \ oscillate; this second attractor is not identified by previous iterative stable motif reduction methods. For a detailed discussion of the nature of these oscillations and their implications for the completeness of the attractor repertoire identified by iterative stable motif reduction, see [55]. Such oscillations motivate us to propose more robust automated methods that can identify oscillations that fail to activate stable motifs. These methods make practical what previously was impractical: the identification of the attractor repertoire of ensembles of Boolean systems, including the identification of complex attractors. Overview of the attractor identification algorithm
The following subsections describe the theoretical foundations of each step of our proposed attractor identification algorithm, which builds upon the maximal trap space identification algorithm of [53] as implemented in the PyBoolNet Python library (https://github.com/hklarner/PyBoolNet/), and the recursive reduction approach of [52]. We give a visual overview of the algorithm in the conclusion of this section. In this work, we focus on the mathematical results underpinning the library; details of the implementation and the code itself are available at https://github.com/jcrozum/StableMotifs/. Our iterative stable motif approach to attractor identification closely follows that of [52] in which stable motifs are recursively used to produce reduced Boolean systems (and corresponding reduced logical hypergraphs) until no additional stable motifs remain (see Figure 3). Our approach has one key difference, however: at each stage in the iteration, we identify whether the reduced system permits any complex attractors that do not activate any additional stable motifs. We call such attractors motif-avoidant and call reduced systems with motif-avoidant attractors terminal . The remainder of this section outlines this process in detail.
Identifying terminal reduced networks using time reversal and drivers
In general, testing for terminality requires analysis of the system’s state transition graph. This approach, however, is computationally expensive. To avoid excess computation in special cases in which terminality can be determined without state-space simulation, we conduct a series of simple tests using necessary or sufficient conditions for terminality for every reduced system
𝐺 = 𝑅𝑒𝑑(𝑓, 𝑃) for
𝑃 ∈ 𝛴 . The first and simplest of these tests is to examine the out-degree of 𝑃 . If 𝑃 has out-degree zero, then 𝐺 does not contain stable motifs and is trivially terminal. A second test is supplied by Corollary 2 (see Appendix A). We consider the set 𝛥 of all virtual nodes that individually drive any stable motif of 𝐺 . If 𝐺 is terminal, then by Corollary 2, all members of 𝛥 must be inactive in some attractor of 𝐺. For example, if
𝑋 = 0 and
𝑌 = 1 each drive a stable motif of 𝐺 , then there must be an attractor with 𝑋 = 1 and
𝑌 = 0 fixed; otherwise 𝐺 would not be terminal. A second example is given in Figure 4: as virtual node (C,1) is a driver of the stable motif marked in green, it must be inactive in any motif-avoidant attractor of the system. Before sampling the STG to identify the attractors of 𝐺 in which all of 𝛥 is inactive, we first attempt to rule out the existence of any such attractor in 𝐺 . We do this by considering the negated set ¬𝛥 = {¬𝐼: 𝐼 ∈ 𝛥} and its LDOI. In particular, if 𝐺 is a terminal reduced system, then it has an attractor in which all virtual nodes in ¬𝛥 are active, or said algebraically, ∧ O∈¬w 𝜎 O (𝑋) = 1 holds for each state 𝑋 in the attractor. This implies that the nodes in the LDOI of these virtual nodes are active as well, i.e. ∧ O∈¬w (∧ Z∈xqyO(O) 𝜎 Z (𝑋)) = 1 holds for each state 𝑋 in the attractor as well. Furthermore, if each variable value in ¬𝛥 is stably fixed in some attractor, then the update values of those variables are also fixed. This is equivalent to requiring that ∧ O∈¬w 𝐹 O (𝑋) = 1 holds true for each state in the attractor. Thus, combining these various algebraic expressions, we find that 𝐺 can only be non-terminal if the expression 𝑅 =∧
O∈¬w (𝜎 O (𝑋) ∧ 𝐹 O (𝑋) ∧ (∧ Z∈xqyO(O) 𝜎 Z (𝑋))) is equal to in at least one state that has all stable motifs of 𝐺 inactive. Furthermore, all states belonging to motif-avoidant attractors of 𝐺 satisfy 𝑅 = 1 . In the example of Figure 4, we find that there is only one single-node driver of the system’s sole stable motif, namely (𝐶, 1) drives the stable motif marked in green in panel A. Thus, the negated driver set is ¬𝛥 = {(𝐶, 0)} and
𝐿𝐷𝑂𝐼(𝐶, 0) is empty, so we find
𝑅 = 𝜎 (],T) (𝑋) ∧ 𝐹 (],T) (𝑋) = ¬𝑋 ] ∧ (𝑋 [ ∨ 𝑋 \ ) . Finally, we note that if 𝑋 is in a stable motif 𝑀 h of the time-reversal of 𝐺 , it can only be in an attractor of 𝐺 if 𝑀 h has an empty contradiction boundary in 𝐺 (i.e. if 𝑀 h is not self-negating when viewed as a set of virtual nodes in the expanded network of 𝐺 ). Thus, when considering states that may be in an attractor of 𝐺 that does not activate any of its stable motifs, we can ignore states belonging to stable motifs of 𝐺 , belonging to self-negating stable motifs of 𝑇𝑅(𝐺) , and states for which 𝑅 is false. We call the collection of all such states the terminal restriction space of 𝐺 . For example, the terminal restriction space of the system in Figure 4 is the uncolored portion of its state transition graph. In practice, focusing on the terminal restriction space can reduce the number of variables that must be simulated by a considerable amount. For example, during the attractor identification of the T-LGL network of [95], when checking the terminality of a reduced system with 36 variables, we find that the terminal restriction space only involves 5 variables, resulting in a reduction of the search-space by a factor of over two billion. (The analysis of a simplified model of this same system, due to [96], is presented in Appendix C). We have also found many examples of random Boolean networks ( 𝐾 = 2, 𝑝 = 0.5, 𝑁 = 4096 ) in which the terminal restriction space of the initial network is of size zero; a state-space reduction of BTF| . Figure 4: An illustration of methods to refine the portion of the STG that may contain stable motif avoidant attractors using the system 𝑓 [ (𝑋) = 𝑓 \ (𝑋) = ¬𝑋 [ ∧ ¬𝑋 \ ∨ 𝑋 ] ; 𝑓 ] (𝑋) = 𝑋 [ ∧ 𝑋 \ as an example. Panel A presents the system’s expanded network and its sole stable motif (in green). This stable motif has one driver set, {(C,1)} (bold outline). By Corollaries 1 and 2 (Appendix A), any stable motif avoidant attractor must have 𝑋 ] = 0 fixed and must satisfy 𝐹 (],T) = 1 , i.e., either 𝑋 [ or 𝑋 \ must be zero. Panel B shows the expanded network of the time reversed system, with its three stable motifs highlighted in purple, blue, and orange. Panel C shows the STG of the system. The states corresponding to each of the stable motifs in panels A and B are highlighted in matching colors. States with (C,1) active are highlighted with a bold outline; none of these states can be part of a stable motif avoidant attractor. The three stable motifs of the time reversal partition the STG into five subgraphs based on which time reversal stable motifs are active: none, purple, blue only, orange only, or both orange and blue. These subgraphs, highlighted by color, are Garden of Eden spaces of the forward-time system. No attractor of the forward-time system can cross between these regions. Any stable motif avoidant attractor of the system must reside in the uncolored, unbolded portion of the STG. Indeed, the states 100, 000, and 010 form a stable motif avoidant attractor in which 𝑋 ] = 0 is fixed while 𝑋 [ and 𝑋 \ oscillate. Simulation and further reduction of the relevant state space
If we are unable to rule out the possibility that the reduced system 𝐺 is terminal, we conduct a state-space search for attractors that do not activate any stable motifs in 𝐺 . Because we are only concerned with a subset of the state-space, we are able to avoid simulating some of the states for an increase in computational speed. We achieve this by keeping track of a set of system states in the state transition graph of 𝐺 , 𝑆𝑇𝐺(𝐺), that cannot be part of any motif-avoidant attractor. We denote this set of system states 𝐾 . Initially, 𝐾 is the terminal restriction space of 𝐺 . We take each state of the state transition graph that is not in 𝐾 , i.e. each state 𝑋 ∈ 𝑉(𝑆𝑇𝐺(𝐺))\𝐾 and compute its successors using the Boolean update rules. If 𝑋 has a successor in 𝐾 , then 𝑋 and all its ancestor states in 𝑆𝑇𝐺(𝐺) (i.e., all states from which we have previously found a path to 𝑋 ) are added to 𝐾 . Similarly, 𝑋 and its ancestors are added to 𝐾 if there is a stable motif 𝑀 h of 𝑇𝑅(𝐺) active in 𝑋 , but inactive in one of its successors. The subgraph 𝑆𝑇𝐺(𝐺)\𝐾 of 𝑆𝑇𝐺(𝐺) necessarily contains any attractors that are in 𝐺 but none of its network reductions. These attractors are identified as the terminal strongly connected components in 𝑆𝑇𝐺(𝐺)\𝐾 . There are two major computational benefits to this approach to simulating a reduced STG. First, each stable motif of a reduced system 𝐺 constrains at least one variable, so the size of the state space is reduced by at least a factor of two for every stable motif in 𝐺 . Second, for a given STG, the task of finding attractors is equivalent to the task of finding terminal strongly connected components, which scales linearly with the number of edges in the STG [97]. As such, reducing the size of the directed graph that must be searched for attractors generally results in drastically improved performance. The degree to which performance is increased depends upon the final size of the set 𝐾 , which itself depends on the stable motif structure of 𝐺 and of its time reversal. State space reduction by node deletion projection
In some cases, even after the reduction described in the previous subsection the terminal restriction space is too large to simulate in a reasonable amount of time. In these cases, we can obtain information about the number of motif-avoidant attractors by studying a simplified version of the system. Specifically, we use the network reduction method of [98,99], in which variables that do not self-regulate are iteratively “deleted” by substituting their update functions into their successors’ update functions (see Figure 5). This method, which we call deletion projection , provides a projection map, 𝜋 , that is proven to preserve point attractors. It was also shown that for every complex attractor of the projected Boolean system the STG of the original network contains a strongly connected component. This nevertheless might not be a complex attractor because its out-component is not necessarily empty in the original network. We combine the concepts of stable motifs and driver nodes with the deletion projection method of [98,99] to investigate the terminality of a motif-reduced Boolean system. We show that the domain of influence of a set of virtual nodes that specifies a state for every variable in the projected system has a domain of influence in the unprojected system that specifies exactly one state, called a “representative state” in [91] (Lemma 1, Appendix A). Combining this result with Theorem 1 (Appendix A) leads to our main result of this section: the activity of stable motifs in attractors is preserved by the deletion projection (Theorem 2, Appendix A). Figure 5: Deletion projection example. Each panel shows the interaction network, expanded network, and STG of a system before projection (top), after projecting out one variable (middle), and after projecting out two variables (bottom). The attractors of each system are highlighted in green and purple such that attractors project onto attractors of the same color. Each attractor has one active stable motif, highlighted in the same color (green or purple, respectively). As seen here, the projection procedure preserves point attractors exactly. In general, the number of complex attractors (in this case zero) serves as an upper bound on the number of complex attractors in the unprojected system (see [98] for a more thorough discussion of oscillations in projected systems). One may view the action of the projection on the STG as contracting the top four nodes (100, 101, 110, and 111) to a single node (1--), with representative state 101 (transitions among these four nodes ultimately lead to 101); a similar view can be taken of the bottom four nodes. In accordance with Theorem 2, the stable motifs ({(A,1),(B,0)} and {(A,0),(B,1)}, which project to the self-activating virtual nodes (A,1) and (A,0), respectively) are preserved in the sense that, for example, the attractor with (A,1) active in the bottom panel corresponds to an attractor in which the stable motif {(A,1),(B,0)} (which projects onto (A,1)) is active. This property here implies that the value of 𝑋 [ is sufficient to determine which attractor the system attains. In cases for which other simulation methods are computationally prohibitive, we resort to the deletion projection method and application of Theorem 2 in order to investigate terminality of a motif-reduced Boolean system. We find the attractors of the projected system and test whether each projection attractor is inconsistent (in the sense of Theorem 2) with all stable motifs. If none of the projection attractors are inconsistent, then the reduced system from which the projection was constructed is necessarily not terminal. Otherwise, the deletion reduction has motif-avoidant attractors. If the motif-reduced system has a motif-avoidant attractor, it necessarily corresponds to one or more of the motif-avoidant attractors of the deletion projection. In this way, the number of motif-avoidant attractors of the deletion projection places an upper bound on the number of motif-avoidant attractors of the motif-reduced system. As a conclusion to this section we present a visual summary of the attractor identification algorithm in Figure 6. We also present a “by hand” example of the approach on a five-variable system in Appendix D. Figure 6: Flowchart illustrating the procedure for generating a succession diagram. Steps are tagged with colored circles indicating which figures provide additional detail for each step. The method for constructing a succession diagram is similar to that presented in Figure 3, but with the additional step of searching for stable motif avoidant attractors using the methods summarized in Figures 4 and 5, which rely on special properties of driver nodes and Boolean time reversal (Figure 2). In brief, the approach is to iteratively consider every possible sequence of stable motif activation (of any length, including zero). The system is simplified under the assumption that all stable motifs in the sequence are active. We then determine whether there are attractors that do not activate any stable motifs of the reduced system. If such an attractor exists, or if the reduced system does not contain stable motifs, the sequence of stable motifs is called terminal. The set of terminal sequences is in surjective correspondence with the attractors of the system. Application to decision-making in empirical biological network models
In this section we illustrate how our methods can be applied to validated Boolean models of biological networks to analytically connect regions of state space to decision-making (points-of-no-return) in their dynamics and to subnetworks (stable motifs) in the underlying interaction network. The biological network models we focus on are a model [96] of a type of white blood cell cancer (T-LGL leukemia) and a model of the mammalian cell cycle [100]. In Appendix C we apply our methods to the T-LGL leukemia model. In particular, we identify and characterize the Garden of Eden spaces, illustrate an informative partitioning of the state space, and fully describe the attractor basins. In the subsection below, we apply our methods to the Phase Switch of the cell cycle model.
Decision-making in the cell-cycle network dynamics
The Phase Switch is a tri-stable molecular circuit that drives mitosis. Its three steady states mark three stages of the cell cycle: G1, G2, and the spindle assembly checkpoint (SAC). Under various biologically relevant conditions, such as coupling to other molecular circuits, the Phase Switch oscillates between these stages in the order they occur in the cell cycle [54] -- here, however, we study this switch in isolation. Further details of this model, including the update functions and stable motifs, are given in Appendix E. Our focus is on how details of a system’s decision-making capacity can be explored by analyzing the system’s time reversal in conjunction with the succession diagram. We emphasize insights gained from the stable motifs of the time-reversed system. In principle, more granular initial condition tracing is possible via the full succession diagram of the time-reversed system. Figure 7 shows how the succession diagram and Garden of Eden spaces provide a concise summary of the irreversible commitments in state space; this is especially useful when the state transition graph is too large or too dense to analyze directly (compare panels B and E). The compressed state space illustration of Panel E unites information from the original and time-reversed system. It is equally or more informative than the full state transition graph (panel B). The compressed succession-diagram representation has the benefit that it can be constructed for very large networks whose state transition graph cannot be built. Indeed, as we will demonstrate in the next section, we can build succession diagrams for networks whose state transition graphs have more nodes ( ≈ 2
BTTT ) than can physically exist in the observable universe ( ≈ 2 |TT , assuming one state per Planck volume and disregarding graph edges). An important feature of this analysis is that decisions can be ascribed to specific subgraphs of the interaction network. These subgraphs correspond to the strongly connected components of the parity expanded network that do not have any parity-invariant subgraphs. These can provide important biological insights. For example, from the succession diagram (Figure 7C), we see that the P1 stable motif cannot be active during the spindle assembly checkpoint (SAC) phase. Driver node analysis of the Cdk1-Wee1 feedback associated to the P1 stable motif (Figure 7D) indicates that knockout of Cdk1 alone blocks the spindle assembly checkpoint. Following P1 commitment, the choice between G1 and G2 is decided by the activity of the Cdc25A-CyclinA feedback loop (P0) in competition with the Cdh1-CyclinA feedback loop (P6) when pAPC is inactive (P5).
Figure 7: Parity and time-reversal analysis elucidates state-space decision-making in the Phase Switch model of [100]. Appendix E provides further details about the model, including its update functions. The interaction network is presented in panel A. The full state transition graph is indicated in panel B. The three point attractors of this model are the three large nodes highlighted in blue, yellow, and red, and correspond to the G1, G2, and spindle assembly checkpoint (SAC) phases of the cell cycle. The basins of attraction exclusive to each attractor are highlighted in the matching color. Grey nodes have paths to multiple attractors; the path taken depends upon the stochastic update order. The stable motifs of the time reversed system identify unstable “Garden of Eden” regions of the state space (highlighted in green) which cannot be (re)entered from the outside. A darker shade of green indicates the overlap of multiple Garden of Eden spaces. The three large green nodes are the “Garden of Eden” states of the system, and are obtained analytically as the fixed points of the time-reversed system. Panel C shows the stable motif succession diagram of the system, with stable motifs of the network and reduced network denoted by labels P0-P6. The colors indicate, by the same scheme as in panel B, which attractors are possible after commitment to each stable motif. Panel D shows the subnetworks and node states associated with these stable motifs (blue: active/on, grey: inactive/off). These subnetworks and their states are obtained as strongly connected components of the expanded network that do not contain parity-invariant subgraphs. In panel E, we illustrate how the succession diagram and garden spaces together describe (analytically) the possible decisions the system can make during its dynamics. Ovals represent specific subsets of the state space. The green oval stands for “Garden of Eden” spaces, with the internal diamonds representing generic states in several categories of these spaces: the darkest green diamond represents Garden of Eden states (fixed points of the time reversed system); the green diamond represents states in the overlap of multiple stable motifs in the time reversed system; the light green diamond represents states in a single stable motif of the time reversed system; and finally, the grey diamond represents states that do not lie in any stable motif of either the forward or backward time system. That states of this last type belong in a Garden of Eden space follows from the definition of stable motifs and the lack of motif-avoidant attractors in this system. Progression from darker to lighter green represents the commitment to exiting Garden of Eden spaces. Grey ovals represent overlapping spaces determined by the stable motif written inside each oval. The diamond symbols in these overlaps mark relevant trap spaces, e.g. the yellow diamond marks the trap space in which motifs P1, P5 and P6 are all active. Each edge corresponds to an irreversible commitment to a smaller trap space. For example, the yellow-blue diamond marking the intersection of the P1 and P5 spaces has two edges: the edge to the yellow diamond indicates a transition to the intersection of the P1, P5 and P6 motif regions, which determines an irreversible commitment to G2 (entry into the exclusive basin of attraction to G2), and the edge to the blue diamond indicates a transition to the intersection of the P0, P1 and P5 motif regions, which determines an irreversible commitment to G1 (entry into the exclusive basin of attraction to G1). Circular symbols indicate the attractors and highlight their position in the narrowest trap space. The graph structure in panel E is isomorphic to the graph structure in panel C, demonstrating that the stable motif succession diagram encapsulates the trajectory of the system in state space.
Random Boolean network results
We study the scaling of the average number of attractors of ensembles of
𝐾 = 2, 𝑝 = 0.5
Kauffman networks under stochastic asynchronous update. This quantity provably scales as a power law bounded by 𝑁 @AB [87]. In this section, we obtain the current best numerical estimate of the scaling exponent value (see Supplemental Material and the Random Boolean Network Application folder at https://github.com/jcrozum/StableMotifs/ for details on ensemble generation and analysis). We generate ensembles for increasing network size 𝑁 (from size 𝑁 = 2 to size
𝑁 = 4,096 ) and apply our attractor identification algorithm to construct the succession diagram for each network and to determine or bound its number of attractors. For greater than 96% of the more than unique networks generated we are able to exactly enumerate the attractors, and this fraction never falls below ( for 𝑁 = 2,048 ) for any value of 𝑁 considered here. For all but ( of the total) of the networks without an exact attractor count we can constrain the number of attractors by using the deletion projection and/or counting the number of maximal stable modules (see Supplemental Material for details and implementation). Of these thirteen networks, the plurality (5) are of size , corresponding to of the networks generated for this size. For each of these remaining thirteen, we use the trivial lower bound of for the number of attractors and choose an upper bound larger than that of the most attractor-rich networks of the same size (see Supplemental Material for details). This makes these networks outliers without introducing undue sensitivity. The scaling of the average number of attractors 〈𝐴〉 as a function of network size 𝑁 is shown in Figure 8. Figure 8: Summary of RBN attractor scaling fits. Symbols indicate the measured number of attractors and the lines represent fits of the form 〈𝐴〉 = 𝑎 + 𝑏𝑁 (cid:134) using nonlinear least-square fitting (see Supplemental Material for implementation details). The fits yield intercept 𝑎 = −0.38 , 𝑎 = 0.44 , and 𝑎 = −0.16 for the exact counts, upper bounds, and lower bounds, respectively. The exponents of the fits ( 𝑐 in 〈𝐴〉 = 𝑎 + 𝑏𝑁 (cid:134) ) are reported in the legend; see the main text or Supplemental Material for bootstrapped confidence intervals. For reference, the sizes of previously considered asynchronous RBN ensembles ([101,102]) and several genetic networks from biology [103–107] are annotated on the horizontal axis. By fitting power law 〈𝐴〉 = 𝑎 + 𝑏𝑁 (cid:134) we find that the exponent is 𝑐 = 0.12 ± 0.05 (one standard deviation; CI [0.04,0.22] ) when fitting only the networks for which we are able to exactly enumerate the attractors (blue circles and curve on Figure 7). The attractor upper bounds (orange symbols and curve) yield an upper bound on the scaling exponent of 𝑐 = 0.20 ± 0.05 ( CI [0.11,0.30] ). To ensure that the considered networks are sufficiently large, we analyze the scaling of the number of maximal stable modules in networks as large as 𝑁 = 16,384 . Because maximal stable modules correspond to disjoint trap spaces in the state space, their number serves as a lower bound on the number of attractors. We merge these lower bounds with the lower bounds for the rest of our ensemble. In practice, the lower bounds are very often in 1-1 correspondence with the number of attractors, as is supported by the good agreement between the exact count of attractors and the lower bounds on the attractor counts. In particular, the scaling of the lower bounds ( 𝑐 = 0.13 ± 0.04 , CI [0.06,0.21] ) is consistent with the attractor scaling for 𝑁 ≤ 4096 and continues at least until
𝑁 = 16,384 , providing strong support that we have probed sufficiently large networks. These networks are of comparable size to many frequently studied genomes (e.g.,
𝑁 ≈ 4,000 for the
E. coli genome, while for the Human genome,
𝑁 ≈ 20,000 ), so any scaling effects beyond
𝑁 = 16,384 are unlikely to be of biological interest. All the scaling exponents are well below the initially conjectured square root scaling ( 𝑐 = 0.5 ) and the theoretical maximum of 𝑐 =𝑙𝑛 4 ( ≈ 1.39 ). For additional details regarding the fitting procedure and error estimation, see Supplemental Material. Overall, we obtain the best current estimate of the exponent of 𝐾 = 2, 𝑝 = 0.5
Kauffman networks under stochastic update. Our analysis represents an 80-fold increase in network size over previous exact or near-exact enumeration analysis of asynchronous
𝑁 − 𝐾
RBNs [101,102]. We find a significantly lower exponent than the original exponent found by Kauffman and those identified in other stochastic updating schemes [87][101,102].
Discussion
Two central questions in complexity science are “what emergent behaviors could a complex system exhibit?” and “what controls a complex system’s selection of one emergent behavior or another?”. While interesting questions from a purely theoretical perspective, they also have important applications in engineering, social science, and medicine. The past several decades have seen a growing number of collaborations between biologists, computer scientists, mathematicians, and physicists that approach these questions using Boolean networks. This work continues that interdisciplinary tradition, using geometric intuitions from physics to prove new results in the mathematics of Boolean dynamics, which we have applied to develop improved computational methods for analyzing complex systems common in biology and other sciences. Coming full circle, these methods have yielded new results in the statistical mechanics of random Boolean networks. Our methods allow efficient identification of the subspaces of the state space where robust commitments happen, and connect these decision-making spaces to subnetworks in the underlying interaction network. Time reversal identifies a previously unexplored type of decision-making: the commitment to exit a “Garden of Eden” space, eliminating the option to ever return to that space. We associate each decision with the stabilization or de-stabilization of strongly connected subnetworks that do not contain any parity-invariant parts. The parity symmetry of Boolean systems has led us to propose a new definition of the expanded network that allows us to view the dynamics as a process on the expanded network itself. It explicitly presents the parity transformation implicitly underlying previous definitions. The expanded network, much like the state transition graph, is parity-invariant (up to node relabeling) and completely encodes the system’s dynamics, but its number of nodes grows linearly, rather than exponentially, with the system dimension. By treating a system and its parity inversion simultaneously, we can study its response to large knock-out and knock-in (constitutive activity) perturbations within the same framework, never needing to distinguish the two at the abstract level. This stands in contrast to the similar ideas first presented in [28], upon which we build. In this earlier construction, the expanded network is a static structure that comprises three types of nodes: regular virtual nodes, negated virtual nodes, and composite nodes. The use of three types of nodes and lack of dynamical structure required careful case-splitting and multiple representations of the system in order to study system perturbations using the expanded network (see e.g., [57]). In contrast to the earlier, static, network, our construction here is a dynamical object with a uniform node set, greatly simplifying its use as a complete description of the system’s dynamics. Rather than distinguish nodes that correspond to on and off states at a definitional level, we can distinguish these nodes by their behavior under a global parity transformation. This view elucidates the role of the expanded network as a parity invariant and of its stable motifs as strongly connected components that have no parity-invariant subgraphs. These insights have enabled our formal proofs of novel results that relate network reduction methods, driver sets, and stable motifs -- results which we have leveraged to develop a fast attractor-finding method. We have presented for the first time the time reversal of a stochastic update Boolean system and demonstrated its usefulness in analyzing the forward-in-time dynamics. Just as stable motifs of a system describe stable spaces in the dynamics, stable motifs in the time reversal of that system describe unstable spaces. This observation is especially helpful in eliminating states when searching for attractors via direct STG construction, or in reducing the number of relevant initial conditions for study. In addition, it demonstrates an important property: the activity of any stable motif of a system or its time-reversal is constant for any attractor. In this way, time-reversal and parity elucidate the “attractor-conserved” quantities of a system’s dynamics. These attractor-conserved quantities are only conserved within attractors and may initially vary. Nonetheless, it is surprising to note that in an inherently stochastic system there is a well-defined notion of time-reversal that yields asymptotic conservation laws. Also surprising is that the time-reversal transformation, though defined globally, is implemented via 𝑓 = ¬𝑓 (𝑋 = ¬𝑋 ) and can thus be applied one variable at a time. On the state transition graph, this has the effect of reversing only edges corresponding to changes in the affected variables. One avenue for future investigation is to consider the applications for such transformations. These partial time-reversals are in analogy to the previously discussed partial parity transformation that allows any stable motif to be realized within a single parity layer [48]. By combining new results stemming from the parity view of the expanded network with our novel time-reversal construction, we developed a method for fast attractor identification in stochastic update Boolean systems. We employed this new method to explore the scaling in attractor number for asynchronous 𝑁 − 𝐾
Kauffman networks. Using these methods, we probed the power law scaling of the
𝐾 = 2 critical ( 𝑝 = 0.5 ) Boolean networks under asynchronous update. Our new techniques allowed us to find or bound the number of attractors in these networks for sizes larger than ever before considered. The power law scaling we observe is much lower (by a factor of about 10) than the theoretical maximum of 𝑙𝑛4 , and also lower than the originally conjectured 0.5 scaling exponent. The low average number of attractors we find ( 〈𝐴〉 ≈ 4 for networks with
𝑁 ≈ 4,000 ) is consistent with previous results on the average number of point attractors and stable synchronous attractors in this ensemble [67,92]. In numerical studies of power law scaling, it is always possible that asymptotic behavior lies outside the simulated region. In this case, however, the scaling phenomenon is of interest primarily as a model of finite-size genomes, and indeed our analysis covers network sizes that are comparable to the size of the human genome. The relatively slow growth of the average number of attractors compared to (i) the originally conjectured 0.5 scaling exponent and (ii) the current cell type scaling estimate of 0.70 (0.88/1.26 = 0.70, based on the experimental data from [45]) has several possible explanations that suggest follow-up investigations. One possibility is that timing-specific attractors contribute significantly to the cell-type scaling, implying that gene-regulatory synchrony plays a crucial role in a cell’s ability to differentiate. Alternatively, the scaling of the attractor number under stochastic update might vary between critical RBN ensembles (e.g. it may differ in ensembles using canalyzing functions or threshold functions), and some of these other ensembles could more accurately reproduce the observed cell type scaling. Analysis of the attractor number scaling in other RBN ensembles using our approach should help answer this question. It is also possible that gene regulatory networks of living organisms have evolved to increase the number of robust attractors, a process which is not fully captured by these ensembles of RBN. The methods developed here can be readily applied to the numerous published Boolean models of biological systems to elucidate their full attractor repertoire. Our framework can also bring further insight into a variety of models that could be reformulated as Boolean models. For example, the quenched Glauber dynamics set on a network [108], models of binary opinion propagation [5], or the Hopfield model [21] can be expressed with threshold Boolean functions. The stable motifs of these systems, and correspondingly the trap spaces of their dynamics, can be identified as particular instances of strongly connected subgraphs. Indeed, in the Watts model of opinion propagation the percolation of an opinion depends on the existence of a strongly connected subgraph of early adopters, who can be influenced by a single neighbor to adopt the opinion [109]. We expect that future adaptation of our methods to these models will be able to reveal rare attractors (metastable states). Apart from the direct application to Boolean models we have emphasized here, time reversal and parity play a role in describing the fundamental logical relationships between entities in complex systems more generally. In this view, the logical parity inversion and time reversal of a system describe a coarse-grained and discretized version of the dynamics, which in turn provides insight into the dynamics of more detailed models (see e.g. [47,61] for further discussion). While the extent to which our key results generalize beyond the stochastic Boolean systems presented here remains an open question, we are encouraged by preexisting analogs of expanded networks in multi-level systems [58] and ODEs [59], as well as by results connecting logic-based models to ODEs [47,60,61,110]. Though our focus here is at the level of interaction logic, our results suggest a new approach to analyzing complexity: studying the relationship of a complex system to its logical parity inversion and time reversal to constrain the system’s repertoire of emergent behaviors.
Acknowledgements
This work was supported by NSF grants PHY 1545832, MCB 1715826, IIS 1814405 to R.A. and the Stand Up to Cancer Foundation/The V Foundation Convergence Scholar Award to J.G.T.Z. (D2015-039).
Appendix A: Formal statements and proofs of key results
Theorem 1 and its corollaries relate an attractor to the domains of influence of driver sets that activate in said attractor.
Theorem 1 : If an attractor 𝒜 contains a state 𝑋 in which the set of virtual nodes 𝑆 is active, then 𝒜 contains a state for which 𝐷𝑂𝐼(𝑆) ∪ 𝑆 is active.
Proof of Theorem 1 : By definition, 𝒜 contains all system states reachable from 𝑋 by all (stochastic asynchronous) update orders. In particular, it contains the set 𝒜 ’ of all states reachable by update orders in which the variables described by 𝑆 are never updated and remain fixed at their values in 𝑋 . Consider the modified system obtained by replacing the update functions 𝑓 for 𝑋 by the constant function 𝑓 (𝑋) = 𝑠 for all (𝑖, 𝑠) ∈ 𝑆 . By definition, the set 𝐷𝑂𝐼(𝑆) consists of all virtual nodes (𝑗, 𝑢) ∉ 𝑆 for which 𝑋 (cid:144) = 𝑢 is fixed in all attractors of the modified system. It may also contain elements of 𝑆 , but this is irrelevant when considering 𝒜 ’ because 𝑆 is active in all states in 𝒜 ’ by definition. In particular, 𝐷𝑂𝐼(𝑆) ∪ 𝑆 is active in all attractor states of the modified system. Because 𝑆 is active in both the modified system and along the paths from 𝑋 to the states in 𝒜 ’, it follows that all paths that start from 𝑋 in the modified system can reach 𝒜 ’. Therefore, 𝒜 ’ contains an attractor of the modified system and thus a state for which 𝐷𝑂𝐼(𝑆) ∪ 𝑆 is active.
Corollary 1 : If
𝐷𝑂𝐼(𝑖, 𝑠) contains a stable module (or motif), 𝑀 , then in every attractor either 1) 𝑀 is active, or 2) 𝑋 is fixed in state ¬𝑠 (or both). Proof of Corollary 1 : If 2) does not hold in an attractor 𝒜 , then 𝒜 contains a state 𝑋 in which 𝑆 = {(𝑖, 𝑠)} is active. Thus, by Theorem 1, 𝒜 contains a state in which 𝐷𝑂𝐼(𝑖, 𝑠) , and 𝑀 in particular, activates. The permanency of the activation of stable modules (and motifs) implies that if 𝑀 is active in any state of 𝒜 , it is active in all states of 𝒜 . Therefore, 1) and 2) cannot both fail to be true. Corollary 2 : Let 𝛥 be a set of single-node drivers of a collection of mutually consistent stable modules (or motifs). If there exists an attractor 𝒜 in which none of these stable modules (or motifs) are active, then ¬𝛥 is not self-negating (its contradiction boundary is empty). Corollary 2 follows immediately from Corollary 1, as discussed in the main text. Lemma 1 and Theorem 2 formalize the relationship between deletion projection of [98,99] and stable motif analysis. Lemma 1 : Let 𝑓 be the update functions of a stochastic asynchronous Boolean system and 𝑓 be a deletion projection of 𝑓 with projection map 𝜋 . Let 𝑋 (cid:145)(cid:146)(cid:147) be any state of the reduced system and 𝑆 be the set of virtual nodes in the expanded network of 𝑓 that are (1) also in the expanded network of 𝑓 , and (2) active in 𝑋 (cid:145)(cid:146)(cid:147) . Then 𝑆 has the property that 𝐷𝑂𝐼(𝑆) ∪ 𝑆 is active in exactly one state 𝑋 of the unreduced system. Proof of Lemma 1 : The entire set of virtual nodes
𝐷𝑂𝐼(𝑆) ∪ 𝑆 is active in the attractors of the modified system in which 𝑆 is kept always active (this follows from the definition of domain of influence). Insofar as we consider only transitions that do not involve updates of the nodes constrained by 𝑆 , the dynamics are unchanged by modifications of the update rules of these variables. In particular, consider the modified system constructed in the proof of Theorem 1 by replacing the update rules 𝑓 for 𝑋 by the constant function 𝑓 (𝑋) = 𝑠 for all (𝑖, 𝑠) ∈ 𝑆 . The deletion projection procedure does not depend on the form of these modified update functions because their corresponding variables are not deleted; thus in this modified system, the deletion project takes the form 𝑓 (𝑋) = 𝑠 . Then, following the convention that constant variables are removed, the deletion projection completely reduces the system to the 0-dimensional system, via the projection map that sends every 𝑋 to 𝑠 . The point attractor correspondence demonstrated in (Naldi et al. 2011) then implies that the original system also has exactly one attractor and that it is a point attractor satisfying 𝑋 = 𝑠 for all (𝑖, 𝑠) ∈ 𝑆 . Because the point attractor specifies a value for every node, it follows that 𝐷𝑂𝐼(𝑆) ∪ 𝑆 is active in exactly one state. This lemma is closely related to various results in [91] regarding what are therein called “representative states”. Lemma 1 tells us that fixing the state of the variables specified by S (which is given by the state of all variables of the projected system) will result in a steady state in the unprojected system. It leads to a key result about the relationship between the deletion projection method and stable motifs: the activity of stable motifs in attractors of the projected system is indicative of their activity in attractors of the original system and vice versa. Theorem 2 makes this correspondence precise.
Theorem 2 : Let 𝑓 be the update functions of a stochastic asynchronous Boolean system that has a stable module (or motif) 𝑀 and an attractor 𝒜 . Further, let 𝑓 be a deletion projection of 𝑓 with projection map 𝜋 . Let 𝒜 (cid:145)(cid:146)(cid:147) be an arbitrary attractor of the deletion projection with 𝒜 (cid:145)(cid:146)(cid:147) ⊆ 𝜋(𝒜) . The following implications hold: 1. If 𝑀 is active in 𝒜 then each state of 𝒜 (cid:145)(cid:146)(cid:147) can be lifted to a state with 𝑀 active (i.e., 𝑀 is active in some state of 𝜋 hD (𝑋 (cid:145)(cid:146)(cid:147) ) for every state 𝑋 (cid:145)(cid:146)(cid:147) ∈ 𝒜 (cid:145)(cid:146)(cid:147) ) and 2. If any state 𝑋 (cid:145)(cid:146)(cid:147) in 𝒜 (cid:145)(cid:146)(cid:147) can be lifted to a state with 𝑀 active (i.e., if 𝑀 is active in some element of 𝜋 hD (𝑋 (cid:145)(cid:146)(cid:147) ) ), then 𝑀 is active in 𝒜 . Proof of Theorem 2 : First, note that 𝒜 (cid:145)(cid:146)(cid:147) exists by Theorem 1, part 3 in [98]. Next, the two parts of the theorem are proved separately. 1. Assume 𝑀 is active in 𝒜 . Consider any state 𝑋 (cid:145)(cid:146)(cid:147) ∈ 𝒜 (cid:145)(cid:146)(cid:147) . Because 𝒜 (cid:145)(cid:146)(cid:147) ⊆ 𝜋(𝒜) , there exists 𝑋 ∈ 𝒜 such that that 𝜋(𝑋) = 𝑋 (cid:145)(cid:146)(cid:147) . Because 𝑀 is active in 𝒜 , it is in particular active in 𝑋 . Therefore, 𝑋 is a state in 𝜋 hD (𝑋 (cid:145)(cid:146)(cid:147) ) in which 𝑀 is active. 2. Assume 𝒜 (cid:145)(cid:146)(cid:147) contains a state 𝑋 (cid:145)(cid:146)(cid:147) such that 𝑀 is active in some state 𝑋 ∈ 𝜋 hD (𝑋 (cid:145)(cid:146)(cid:147) ) . Consider the set 𝑆 of 𝑑𝑖𝑚 𝑓 virtual nodes in the expanded network of 𝑓 that are active in all states that map to 𝑋 (cid:145)(cid:146)(cid:147) under 𝜋 . By Lemma 1, 𝐷𝑂𝐼(𝑆) ∪ 𝑆 is active in exactly one state. Because 𝑀 and 𝑆 are both active in 𝑋 and 𝑀 is describes a trap space, it follows that 𝐷𝑂𝐼(𝑆) ∪ 𝑆 cannot contradict 𝑀 . Because 𝐷𝑂𝐼(𝑆) ∪ 𝑆 is active in only one state, this implies that 𝑀 is active in 𝐷𝑂𝐼(𝑆) ∪ 𝑆 . Then from 𝑋 (cid:145)(cid:146)(cid:147) ∈ 𝒜 (cid:145)(cid:146)(cid:147) ⊆ 𝜋(𝒜) , it follows that 𝒜 contains at least one state in 𝜋 hD (𝑋 (cid:145)(cid:146)(cid:147) ) , and thus it contains a state in which 𝑆 is active. Thus, by Thm. 1, 𝒜 contains a state in which 𝐷𝑂𝐼(𝑆) ∪ 𝑆 , and 𝑀 in particular, is active. Because 𝑀 describes a trap space, its activation is irreversible. Therefore, 𝑀 is active in all of 𝒜 . Appendix B: Connecting stable motif succession and Garden-of-Eden spaces to decision-making in state space
We consider the example introduced in Figure 3; it is a six-variable system with update functions 𝑓 [ (𝑋) = 𝑋 \ 𝑓 \ (𝑋) = 𝑋 [ 𝑓 ] (𝑋) = 𝑋 [ ∨ ¬𝑋 q 𝑓 q (𝑋) = 𝑋 ] 𝑓 r (𝑋) = 𝑋 \ ∧ 𝑋 s 𝑓 s (𝑋) = 𝑋 r The stable motifs and stable motif succession diagram of this example are indicated in Figure 3. The stable motifs of the time-reversed system, and thus the “Garden of Eden” spaces of the original system, are 𝑋 [ = 0 , 𝑋 \ = 1 and 𝑋 [ = 1 , 𝑋 \ = 0 . Figure 9 illustrates that the state transition graph of the system can be partitioned into Garden-of-Eden spaces (those that have a red background) and trap spaces in which one or more stable motifs are locked-in (indicated by green, blue, purple and pink background colors, which match the colors in Figure 3). The state transition graph indicates that the system irreversibly leaves the Garden-of-Eden spaces and commits to increasingly restricted trap spaces. A particular example of a possible trajectory is to start from a state in the top right group of states (red box with A=0, B=1), enter the bottom right group of states (blue box with A=B=1), then enter the group of states with pink background (E=F=1), then converge into the attractor marked in yellow. This trajectory corresponds to the top branch of the stable motif succession diagram of Figure 3E. Figure 9: The state transition graph of the example from Figure 3. The state nodes are labeled according to the values of the variable 𝑋 [ through 𝑋 s , in alphabetical order. The states that make up the system’s three attractors are highlighted in yellow, orange, and brown. State transitions are represented by grey dashed arrows, and the transitions between various groupings of nodes, or between nodes within a group, are summarized by solid black arrows. The states are grouped according to the nodes of the succession diagram in Figure 3E, with matching colors (green, purple, pink, and blue). Note that each of these state groups is a trap space in the sense that there are no transitions out of any of the stable motif regions. In addition, three additional groups are highlighted in red: the states in which 𝑋 r = 0 and 𝑋 s = 1 ; the states in which 𝑋 [ = 0 and 𝑋 \ = 1 ; and the states in which 𝑋 [ = 1 and 𝑋 \ = 0 . These three groups are the states in which the time reversed system’s three stable motifs are active. As such, these groups are “Garden of Eden” spaces. Indeed, there are no state transitions into any red group from outside that group. Appendix C: Analysis of a T-LGL model using time-reversal
To illustrate the new insights that the expanded network framework and the time-reversal system can bring, we analyze a simplified network model [96] of a type of white blood cell cancer (T-LGL leukemia) shown in Figure 10. This model considers the relationships between five proteins (denoted “cer”, “disc”, “fas”, “flip”, and “s1p”) and the process of cell death (apoptosis, denoted “apo”). The governing update functions are as follows: 𝑓 (cid:150)(cid:151)(cid:152) (𝑋) = 𝑋 (cid:147) ∨ 𝑋 (cid:150)(cid:151)(cid:152) 𝑓 (cid:134)(cid:146)(cid:145) (𝑋) = ¬𝑋 (cid:150)(cid:151)(cid:152) ∧ 𝑋 (cid:154)(cid:150)(cid:153) ∧ ¬𝑋 (cid:153)D(cid:151) 𝑓 (cid:147) (𝑋) = (¬𝑋 (cid:150)(cid:151)(cid:152) ∧ 𝑋 (cid:154)(cid:150)(cid:153) ∧ ¬𝑋 (cid:154)@ ) ∨ (¬𝑋 (cid:150)(cid:151)(cid:152) ∧ 𝑋 (cid:134)(cid:146)(cid:145) ) 𝑓 (cid:154)(cid:150)(cid:153) (𝑋) = ¬𝑋 (cid:150)(cid:151)(cid:152) ∧ ¬𝑋 (cid:153)D(cid:151) 𝑓 (cid:154)@ (𝑋) = ¬𝑋 (cid:150)(cid:151)(cid:152) ∧ ¬𝑋 (cid:147) 𝑓 (cid:153)D(cid:151) (𝑋) = ¬𝑋 (cid:150)(cid:151)(cid:152) ∧ ¬𝑋 (cid:134)(cid:146)(cid:145) Previous analysis of the model indicated that it has two attractors, one corresponding to a T-LGL (cancerous) state and one corresponding to the normal behavior (apoptosis). We find that the system has three stable motifs, indicated in Figure 10A. Two of these are mutually consistent and describe trap spaces that contain the cancerous state, while the third contains the apoptotic state. The time-reversed system contains two stable motifs (Figure 10B), which together partition the state space into four regions (Figure 10C). The virtual node sets that define two of these regions (R1 and R2) are self-negating, meaning that these regions cannot contain attractors. The stable motif {(apo,1)} is active in all states of region R4. The LDOI of this stable motif specifies that fixing apo=1 drives the rest of the variables into their inactive states (see the edges originating from the blue virtual node in Figure 10A). Thus, region R4 contains the point attractor corresponding to apoptosis, which is the larger bold-outlined dark-blue node in Figure 10E. The stable motifs colored in red in Figure 10A are active in the subset of region R3 that corresponds to irreversible commitment to the T-LGL point attractor. The flow between the four regions can be readily determined (see thick dashed arrows in Figure 10E). The flow between regions indicates that all system decisions about which attractor to select occur in regions R1 and R3. Any trajectory in R1 may result in the initiation of the apoptosis process ( 𝑋 (cid:150)(cid:151)(cid:152) = 1 ) and thus exit R1 via R2, which ultimately leads to the apoptotic attractor. Otherwise, trajectories eventually exit into region R3. In region R3, there are two possibilities: one of the TLG-L stable motifs activates, or the system escapes the region into R4 and the apoptotic attractor. This second scenario is possible because the restriction of the system to R3 (by setting 𝑋 (cid:150)(cid:151)(cid:152) = 0 ) gives rise to a so-called conditional stable motif [54] {(cer,1),(fas,1),(s1p,0)}, that is self-sustaining in R3, but ultimately drives (disc,1), which in turn drives (apo,1) and exit from R3 into R4, where it is no longer self-sustaining. The four states in which {(apo,0),(cer,1),(fas,1),(s1p,0)} is active (colored in light blue in Figure 10E) are therefore in the exclusive basin of attraction of the apoptosis attractor (i.e. the set of states from which every trajectory leads to the apoptosis attractor).Thus, by combining the inter-region flow with the irreversibility of stable motifs and the properties of driver sets, we recapitulate the previously obtained exclusive basins of the two attractors (in light blue and dark blue for the apoptosis attractor and in red for the cancerous attractor in Figure 10E). In summary, analysis of the expanded network and of the time-reversed system allows an efficient and informative partitioning of the state space and identifies which of the partitions may and which may not contain an attractor. Figure 10: Time reversal and domain of influence analysis of the simplified T-LGL leukemia network of [96]. Panel A depicts the expanded network of the system with its stable motifs highlighted in red and blue, respectively (the two partially shaded red nodes and the dotted hyperedges connected to them indicate that a stable motif is formed even when either (fas,0) or (flip,1) is omitted). Panel B depicts the system’s time reversal, with its stable motifs highlighted in green and brown, respectively. The stable motifs of the system’s time reversal allow us to partition the state-space into four regions according to which time-reversal stable motifs are active: R1) 𝑋 (cid:134)(cid:146)(cid:145) ∧ 𝑋 (cid:153)D(cid:151) = 1 , 𝑋 (cid:150)(cid:151)(cid:152) = 0 , R2) 𝑋 (cid:134)(cid:146)(cid:145) ∧ 𝑋 (cid:153)D(cid:151) = 1 , 𝑋 (cid:150)(cid:151)(cid:152) = 1 , R3) 𝑋 (cid:134)(cid:146)(cid:145) ∧ 𝑋 (cid:153)D(cid:151) = 0 , 𝑋 (cid:150)(cid:151)(cid:152) = 0 , and R4) 𝑋 (cid:134)(cid:146)(cid:145) ∧ 𝑋 (cid:153)D(cid:151) = 0 , 𝑋 (cid:150)(cid:151)(cid:152) = 1 . The properties of Garden of Eden spaces (which include R1, R1 ∪ R2, and R1 ∪ R3), imply that no attractor lies in more than one of these regions. Any attractor in R1 or R2 has {(cer,1),(s1p,1)} active, but this set of virtual nodes is self-negating (its LDOI is depicted in yellow in panel C, with the white nodes with yellow border indicating the contradiction boundary). In any attractor in R3, either {(s1p,0),(apo,0)} or {(cer,0),(apo,0)} (or both) is active. The first of these is self-negating, but the LDOI of the second of these (depicted in yellow in panel D) fixes all variables, implying that R3 contains a single attractor: the point attractor corresponding to the highlighted virtual nodes. By a similar argument, R4, defined by the activity of either {(s1p,0),(apo,1)} or {(cer,0),(apo,1)} contains a single point attractor and no oscillations. The partitioning of the (forward time) STG according to the four regions is depicted in panel E, with each region colored according to which time-reversal stable motifs are active therein. States in the STG are colored in dark blue or red according to whether the stable motifs of the corresponding color are active in each state. The two point attractors are highlighted by bold outlines and are larger than other nodes. For the T-LGL attractor (red), the exclusive basin of attraction coincides with the activation of the corresponding stable motifs. In the apoptosis attractor (blue), the exclusive basin of attraction contains all states in which the apoptosis stable motif is active (dark blue), as well as four additional states (light blue) that can be identified by considering stable motifs of the system obtained by fixing 𝑋 (cid:150)(cid:151)(cid:152) = 0 (see text). Grey states may lead to either attractor depending on the stochastic update order. Solid black arrows indicate intra-region state transitions. Dashed arrows indicate inter-region transitions; the light grey arrows connect states in different regions, while the thick dashed lines summarize these state transitions as region to region transitions. A self-loop on a region indicates that it contains an attractor. Notably, the thick dashed lines illustrate the instability property of Garden of Eden spaces: once a trajectory exits the union of R1 and R2, for example, it can never return. Appendix D: Example of generating the stable motif succession diagram
We illustrate the attractor identification algorithm on a five-variable example with update functions 𝑓 [ (𝑋) = (¬𝑋 [ ∧ ¬𝑋 \ ) ∨ 𝑋 ] 𝑓 \ (𝑋) = (¬𝑋 [ ∧ ¬𝑋 \ ) ∨ 𝑋 ] 𝑓 ] (𝑋) = (𝑋 [ ∧ 𝑋 \ ) ∨ 𝑋 q 𝑓 q (𝑋) = 𝑋 r 𝑓 r (𝑋) = 𝑋 q As indicated by colors in the expanded network (top panel of Figure 11), this Boolean system has three stable motifs. Two of these arise from the E-D positive feedback loop and are mutually exclusive. Stable Motif 1 (blue) is formed by the virtual nodes (𝐷, 1) and (𝐸, 1) , either of which can serve as a driver node of the motif. Stable Motif 3 (purple) is formed by the virtual nodes (𝐷, 0) and (𝐸, 0) , either of which can serve as a driver node. Stable Motif 2 (green) is formed by the virtual nodes (𝐴, 1) , (𝐵, 1) , and (𝐶, 1) , the last of which is its driver node. The set of single node drivers is thus 𝛥 = {(𝐷, 1), (𝐸, 1), (𝐷, 0), (𝐸, 0), (𝐶, 1)} , and so the negated single driver node set is ¬𝛥 ={(𝐷, 1), (𝐸, 1), (𝐷, 0), (𝐸, 0), (𝐶, 0)} , which is self-negating by virtue of being self-contradictory. Therefore, at least one of the three stable motifs must lock in. Following the locking in of Stable Motif 1, the expanded network reduces completely, and an "all-on" point attractor is obtained (Attractor 1). If instead Stable Motif 3 locks in, the example of Figure 4 is obtained (Reduced Network 3); this system has one motif-avoidant oscillation (Attractor 3) and one point attractor (Attractor 2) that is reached if the sole stable motif (Stable Motif 3a) of Reduced Network 3 is activated rather than avoided. Finally, we may consider that Stable Motif 2 activates in the original network, yielding the bistable switch shown in Reduced Network 2. This reduced network has two stable motifs (2a and 2b), either of which results in a fully reduced network and yields Attractor 1 or Attractor 2, respectively. Reduced Network 2 has negated single driver set ¬𝛥 = {(𝐷, 1), (𝐸, 1), (𝐷, 0), (𝐸, 0)} , and thus has no motif-avoidant attractors (i.e., either Stable Motif 2a or 2b must eventually activate). Figure 11: Summary of example in Appendix D. The expanded network of the system is shown in the topmost panel, and below it are the three stable motifs it contains (labeled Stable Motif 1-3). After selecting a stable motif, the effects of its activation are depicted below it as a reduced network (labeled Reduced Network 1-3, respectively). This pattern is repeated iteratively until the resulting reduced network contains no further stable motifs, indicating that a minimal trap space has been uncovered. Each minimal trap space is typically small and therefore identifying the attractor(s) it contains is trivial. Attractors are listed along the bottom of the figure, with the binary strings indicating the states of the variables, in alphabetical order, that are fixed in each attractor; an ‘X’ in this string indicates oscillation. At each stage in the reduction process, one must consider whether there are so-called motif-avoidant attractors that exist in a reduced network but fail to activate any of that network’s stable motifs. One such case occurs in Reduced Network 3: Attractor 3 involves the oscillation of two of the nodes in such a way that Stable Motif 3a never activates.
Appendix E: Details of the cell cycle Phase Switch model of Deritei et al. (2016)
The model of the cell cycle Phase Switch we consider here is an 11-variable network of interactions among proteins that govern mammalian cells’ entry and exit from mitosis (cell division). The model was introduced in [100] and further analyzed in [111]. The update functions are given below; we use the entity name and the name of its corresponding variable interchangeably (e.g., CyclinA should be understood to mean 𝑋 ](cid:156)(cid:134)@ in the functions below). f Cdc25C = CyclinA or (CyclinB and Cdk1) f
Cdh1 = not CyclinA and (not CyclinB or not Cdk1) f
Cdk1 = Cdc25C and (CyclinA or CyclinB) and (Cdk1 or not Wee1) f
CyclinA = (CyclinA or Cdc25A) and not (pAPC and Cdc20 or Cdh1 and UbcH10) f
CyclinB = (not pAPC and not Cdh1) or (not Cdc20 and not Cdh1) f pAPC = (pAPC and Cdc20) or (CyclinB and Cdk1) f
Cdc20 = pAPC and not Cdh1 and not Mad2 f
UbcH10 = not Cdh1 or (UbcH10 and Cdc20) or (UbcH10 and CyclinA) or (UbcH10 and CyclinB) f
Mad2 = (not pAPC or not Cdc20) and CyclinB and Cdk1 f
Wee1 = not CyclinA and not CyclinB or not Cdk1 f
Cdc25A =CyclinA and not Cdh1 This system has three point attractors, which correspond to distinct phases of the cell cycle, namely the G1 (first growth) phase, the G2 (second growth) phase, and the spindle assembly checkpoint (SAC). The G1 attractor is made up by Cdc25C = Cdk1 = CyclinA = CyclinB = pAPC = Cdc20 = UbcH10 = Mad2 = Cdc25A = 0, Cdh1 = Wee1 = 1. The G2 attractor is made up by Cdc25C = CyclinA = CyclinB = UbcH10 = Wee1 = Cdc25A = 1, Cdk1 = Cdh1 = pAPC = Cdc20 = Mad2 = 0. The SAC attractor is made up by Cdc25C = Cdk1 = CyclinB = pAPC = UbcH10 = 1, CyclinA = Cdh1 = Cdc20 = Wee1 = Cdc25A = 0. The stable motifs of this system and of its reduced networks are labeled as follows P0 := {(CyclinA, 0), (Cdc25A, 0)} P1 := {(Cdk1, 0), (Wee1, 1} P2 := {(Cdc20, 0), (CyclinB, 1), (Cdc25C, 1), (Cdk1, 1), (Cdh1, 0), (Mad2, 1)} P3 := {(Cdc25, 0), (Cdk1, 0)} P4 := {(CyclinB, 0), (Cdh1, 1)} P5 := {(pAPC, 0)} P6 := {(CyclinA, 1), (Cdh1, 0)}
Supplemental material
Supplemental Material 1: Ipython Notebook and input data for attractor scaling fits.
References Boccara N. Modeling Complex Systems. Springer Science & Business Media; 2010. 2. Schweitzer F. Brownian Agents and Active Particles: Collective Dynamics in the Natural and Social Sciences. Springer Science & Business Media; 2007. 3. Shalizi CR. Methods and Techniques of Complex Systems Science: An Overview. Topics in Biomedical Engineering International Book Series. 2006. pp. 33–114. doi:10.1007/978-0-387-33532-2_2 4. Sayama H. Introduction to the Modeling and Analysis of Complex Systems. Open SUNY Textbooks; 2015. 5. Barrat A, Barthelemy M, Vespignani A. Dynamical Processes on Complex Networks. 2008. doi:10.1017/cbo9780511791383 6. Vespignani A. Modelling dynamical processes in complex socio-technical systems. Nature Physics. 2012. pp. 32–39. doi:10.1038/nphys2160 7. Mora T, Bialek W. Are biological systems poised at criticality? J Stat Phys. 2011;144: 268–302. 8. Ramaswamy S. The Mechanics and Statistics of Active Matter. 2010 [cited 6 Jul 2020]. doi:10.1146/annurev-conmatphys-070909-104101 9. Castellano C, Fortunato S, Loreto V. Statistical physics of social dynamics. Reviews of Modern Physics. 2009. pp. 591–646. doi:10.1103/revmodphys.81.591 10. Motter AE, Timme M. Antagonistic Phenomena in Network Dynamics. Annu Rev Condens Matter Phys. 2018;9: 463–484. 11. Strogatz SH. Nonlinear Dynamics and Chaos. 2018. doi:10.1201/9780429492563 12. Tyson JJ, Novak B. A Dynamical Paradigm for Molecular Cell Biology. Trends Cell Biol. 2020;30: 504–515. 13. Railsback SF, Grimm V. Agent-Based and Individual-Based Modeling: A Practical Introduction, Second Edition. Princeton University Press; 2019. 14. Toffoli T, Margolus N. Cellular Automata Machines: A New Environment for Modeling. MIT Press; 1987. 15. Barabási A-L. Network Science. Cambridge University Press; 2016. 16. Newman M. Networks. 2010. doi:10.1093/acprof:oso/9780199206650.001.0001 17. Vollmer H. Introduction to Circuit Complexity. Texts in Theoretical Computer Science An EATCS Series. 1999. doi:10.1007/978-3-662-03927-4 18. Glauber RJ. Time - Dependent Statistics of the Ising Model. J Math Phys. 1963;4: 294–307. 19. Godrèche C, Luck JM. Metastability in zero-temperature dynamics: statistics of attractors. Journal of Physics: Condensed Matter. 2005. pp. S2573–S2590. doi:10.1088/0953-8984/17/24/014 20. McCulloch WS, Pitts W. A logical calculus of the ideas immanent in nervous activity. The Bulletin of Mathematical Biophysics. 1943. pp. 115–133. doi:10.1007/bf02478259 21. Hopfield JJ. Neural networks and physical systems with emergent collective computational abilities. Proceedings of the National Academy of Sciences. 1982. pp. 2554–2558. doi:10.1073/pnas.79.8.2554 22. Kauffman SA. Metabolic stability and epigenesis in randomly constructed genetic nets. J Theor Biol. 1969;22: 437–467. 23. Thomas R. Boolean formalization of genetic control circuits. J Theor Biol. 1973;42: 563–585. 24. Montagud A, Traynard P, Martignetti L, Bonnet E, Barillot E, Zinovyev A, et al. Conceptual and computational framework for logical modelling of biological networks deregulated in diseases. Brief Bioinform. 2019;20: 1238–1249. 25. Mendoza L. A network model for the control of the differentiation process in Th cells. Biosystems. 2006. pp. 101–114. doi:10.1016/j.biosystems.2005.10.004 26. Naldi A, Carneiro J, Chaouiya C, Thieffry D. Diversity and Plasticity of Th Cell Types Predicted from Regulatory Network Modelling. PLoS Computational Biology. 2010. p. e1000912. doi:10.1371/journal.pcbi.1000912 27. Sánchez L, Thieffry D. A Logical Analysis of the Drosophila Gap-gene System. Journal of Theoretical Biology. 2001. pp. 115–141. doi:10.1006/jtbi.2001.2335 28. Albert R, Othmer HG. The topology of the regulatory interactions predicts the expression pattern of the segment polarity genes in Drosophila melanogaster. J Theor Biol. 2003;223: 1–18. 29. Thieffry D, Chaouiya C, Sánchez L. From Gradients to Stripes: A Logical Analysis of Drosophila Segmentation Genetic Network. Bioinformatics of Genome Regulation and Structure II. pp. 379–390. 2006 doi:10.1007/0-387-29455-4_36 30. Cohen DPA, Martignetti L, Robine S, Barillot E, Zinovyev A, Calzone L. Mathematical Modelling of Molecular Pathways Enabling Tumour Cell Invasion and Migration. PLOS Computational Biology. 2015. p. e1004571. doi:10.1371/journal.pcbi.1004571 31. Steinway SN, Zanudo JGT, Ding W, Rountree CB, Feith DJ, Loughran TP, et al. Network Modeling of TGF Signaling in Hepatocellular Carcinoma Epithelial-to-Mesenchymal Transition Reveals Joint Sonic Hedgehog and Wnt Pathway Activation. Cancer Research. 2014. pp. 5963–5977. doi:10.1158/0008-5472.can-14-0225 32. Udyavar AR, Wooten DJ, Hoeksema M, Bansal M, Califano A, Estrada L, et al. Novel Hybrid Phenotype Revealed in Small Cell Lung Cancer by a Transcription Factor Network Model That Can Explain Tumor Heterogeneity. Cancer Res. 2017;77: 1063–1074. 33. Eduati F, Jaaks P, Wappler J, Cramer T, Merten CA, Garnett MJ, et al. Patient --