Identifying (un)controllable dynamical behavior in complex networks

We present a technique applicable in any dynamical framework to identify control-robust subsets of an interacting system. These robust subsystems, which we call stable modules, are characterized by constraints on the variables that make up the subsystem. They are robust in the sense that if the defining constraints are satisfied at a given time, they remain satisfied for all later times, regardless of what happens in the rest of the system, and can only be broken if the constrained variables are externally manipulated. We identify stable modules as graph structures in an expanded network, which represents causal links between variable constraints. A stable module represents a system “decision point”, or trap subspace. Using the expanded network, small stable modules can be composed sequentially to form larger stable modules that describe dynamics on the system level. Collections of large, mutually exclusive stable modules describe the system’s repertoire of long-term behaviors. We implement this technique in a broad class of dynamical systems and illustrate its practical utility via examples and algorithmic analysis of two published biological network models. In the segment polarity gene network of Drosophila melanogaster, we obtain a state-space visualization that reproduces by novel means the four possible cell fates and predicts the outcome of cell transplant experiments. In the T-cell signaling network, we identify six signaling elements that determine the high-signal response and show that control of an element connected to them cannot disrupt this response.


Boolean Stable Motifs
Here, we expand upon the concepts discussed in [22,30]. A Boolean network is a discrete-time system composed of variables or nodes Xi (t) that can take values in {0, 1} along with update functions fi = fi ({Xj}). At each timestep, a subset of {Xi (t)} is updated according to X * i = fi ({Xj}), where X * i denotes the updated value of Xi (t). The two most commonly considered update schemes are synchronous update, in which all nodes are updated at each time-step, or general asynchronous update, in which only one randomly selected node is updated at each timestep.
The expanded network of a Boolean system is built by first constructing a pair of virtual nodes for each Boolean variable corresponding to its two states. In a slight abuse of notation, we denote these nodes Xi and ¬Xi. To proceed, we note that the Boolean update functions can be written in disjunctive normal form (DNF) (1.1) The negation ¬fi of each Boolean function can also be written in DNF. Each term of the form M j m=1 X lm Qj q=1 ¬Xp q is called a prime implicant of fi and either corresponds to an existing virtual node, or else we identify the term with a composite node and draw edges from each factor X lm and ¬Xp q to the composite node. From each node in the expanded network that corresponds to a prime implicant of fi (or ¬fi), we draw an edge to the virtual node Xi (or ¬Xi). This scheme implies the following: if there is an edge in the expanded network from M j m=1 X lm Qj q=1 ¬Xp q to Xi, then whenever M j m=1 X lm Qj q=1 ¬Xp q = 1 holds, Xi will always update to 1 in the next update step. Similarly, if there is an edge from M j m=1 X lm Qj q=1 ¬Xp q to ¬Xi and M j m=1 X lm Qj q=1 ¬Xp q = 1 holds then Xi will update to 0 (i.e., ¬Xi becomes 1). Thus in either case, the edge is both a maintenance and driving edge. As such, the Boolean expanded network, as described above, is an expanded network of the type described in the main text. It is also noteworthy that the update functions fi can be reconstructed from the expanded network.
A stable motif in the expanded network is a subgraph S with the following properties: (i) all nodes in S have a parent (regulator) in S or else have a self-loop, (ii) S does not include any contradictory nodes (e.g., both Xi and ¬Xi), (iii) S includes all the parents (regulators) of included composite nodes, and (iv) S has no proper subgraphs satisfying the first three properties (note that condition (iv) implies that S is strongly connected). We call subgraphs satisfying the first three conditions stable modules. Such network components correspond to values of the Boolean variables that are stable, independent of the rest of the system.
As an example, consider the following Boolean network: For convenience, we give the negated update functions as well: We can now construct the expanded network, as the update rules and their negations are written in DNF in Equations 1.2 and 1.3. By inspecting these rules, as written, we see that in addition to the virtual nodes A, B, C, ¬A, ¬B, and ¬C, we will need to construct composite nodes A ∧ B, ¬A ∧ ¬C, ¬A ∧ C, and ¬B ∧ C, for a total of ten nodes. We draw edges to the composite nodes from each of their factors (e.g., there are edges from A and B to A∧B). We also draw edges to the virtual nodes according to the prime implicants of the update rules (e.g., there are edges from A and C to B, and edges from ¬A ∧ C, and ¬B ∧ C to ¬C). The completion of this procedure yields the full Boolean expanded network.
We graphically represent the network given by Equations 1.2, construct its expanded network, and identify stable modules and stable motifs in Figure 1.1. This network has two attractors. The first is the steady state in which all variables are true. The second is the one in which A is false, and B and C oscillate. If A is ever false, it stays false. Plugging A = 0 into the update rules for B and C gives fB = C, fC = ¬C, so C oscillates. The timing of when C is true and when it is false depends on the update scheme, but the update scheme does not change the fact that C oscillates. Node B follows C; thus it oscillates as well. If A is true at some late time, it is not oscillating (because A = 0 self-stabilizes). If A = 1 is stabilized, it can only be because B = 1 is also stabilized, which requires C = 1 as well. This can also be seen from the fact that the only solution of the steady state equations A = A ∧ B, B = A ∨ C, and C = A ∧ B ∨ ¬C is A = B = C = 1. Because these attractors are characterized only by the value of A, there can only be two attractors (one for each value). We further note that the stable motifs that characterize these attractors are independent of update scheme (this is a general feature of stable motifs in Boolean systems [30]). There are two stable motifs, each of which is a stable module, and node C can be added to one of the motifs to create a larger stable module. The first consists of only ¬A, and this stable motif captures the behavior that if A ever becomes false, it remains false. The other two stable modules both contain A, B, and A ∧ B. One also contains C. The A, B, and A ∧ B motif indicates that if A and B are both true, they both remain true, independent of C (e.g., if C = 0 were held fixed by some external control, A = B = 1 would remain stable). The stable module that also contains C identifies A = B = C = 1 as a steady state of the system.

Dynamical Properties of Biomolecular Systems
The dynamics of biomolecular systems are almost always bounded (i.e. the within-cell concentrations of proteins or molecules do not increase without bound). Furthermore, most biomolecular interactions and direct regulatory relationships are monotonic; it is infrequent that a regulator such as a transcription factor or enzyme directly inhibits certain processes (transcription or reactions) and activates others. In many cases, apparent non-monotonic relationships arise from the effects of as-of-yet unknown mediators. This situation is likely in biomolecular networks, in which regulatory relationships are often not mapped to the level of elementary reactions. For example, in considering the regulatory relationshipẋ = w 3 − w − x, we might at first believe that w is a non-monotonic regulator of x, when in reality, a more accurate model of the system might beẋ = y − w − x,ẏ = α w 3 − y , which consists of an incoherent feed-forward loop involving monotonic regulation. This example is especially believable when α 1, because y will respond extremely quickly to changes in w.
There is typically only one additional variable required per instance of non-monotonic regulation. We note that even when a hidden variable cannot be biologically motivated, it may still be introduced without loss of fidelity in order to make the relationships in the model monotonic.
For these reasons, we focus on ODE systems that take the formẋ where Fi is continuous, monotonic in each of its arguments, and bounded.

Expanded Networks and Stable Modules
In this section, we present the key definitions and results of the main text in the more formal style of dynamical systems theory. Throughout this section, we consider a dynamical system (T, M, Φ), where Φ is a monoid action of T on a state space M . That is, T is a set representing the evolution variable (e.g., time), and we place an evolution operator + on T that is associative and for which T has an identity element, e. The space M plays the role of the phase space or state space. The monoid action Φ describes the evolution of the system by mapping a point x ∈ M and a "time" t ∈ T to a new point in M , which is denoted by so that Φ is consistent with the notion of time evolution. It can be given explicitly (as in Boolean networks, where Φ is the time-explicit update function) or implicitly (as in ODE systems, in which Φ is given by the solutions to the ODEs for each initial condition). In the main text, we focus on the case in which T is one-dimensional and M is finite-dimensional, however, the contents of this section apply generally.
To apply our methods, we must have a rule for determining whether or not one "time" follows another. To do this we define a total preorder on a submonoid S of T denoted by and satisfying e s for all s ∈ S, where e ∈ S is the identity. A total preorder is, by definition, a transitive relation for which either or both x y or y x hold.
Of particular interest here is the case of autonomous first order ODE systems, in which T is R + (+), M is a box in R N , Φ is the flow of the defining differential equations, and we take S = T . We also consider Boolean systems, in which T is Z + (+), M is {0, 1} N , Φ is the Boolean update function, and we again take S = T .
The objects of primary importance in the main text are expanded networks, which are constructed for sets of statements regarding the values of variables. These statements may be viewed as subsets of M . To define these networks requires that we have a more general notion of what it means for a trajectory to be confined to a region for a particular duration. We thus define the following sets: Each element of the persistence set is a time for which the trajectory starting at r is confined to R. Note that e is always an element of P (R, r), so the persistence set is never empty. With the above shorthand, we can formally define the maintenance relation as follows: Definition 2. For non-disjoint subsets R1 and R2 of M , we say R1 maintains R2 under S (denoted R1 → R2) if for each r ∈ R1∩R2, P (R1, r) ⊆ P (R2, r) holds with equality holding only when P (R1, r) and P (R2, r) are equal to S. 2 with its two attractors (blue is true, red is false, and the gradient represents oscillation) (ii) the full expanded network of Equations 1.2 (iii) the three stable modules that exist in the network. The "NOT A" stable module (motif) captures the behavior that if A ever becomes false, it remains false. The other two stable modules identify the steady-state attractor of the system.
Informally, the above definition means that R1 maintains R2 if any trajectory starting in R1 ∩ R2 that is confined to R1 for a certain duration is confined to R2 for a strictly longer duration. The maintenance relation induces a directed graph on any collection of regions of M . We next develop the notion of expanded networks in the special case in which there are no composite nodes. We then prove the stabilization results of the main text for this special case. We later show that the case in which composite nodes are present follows as a simple extension of that result.
Definition 3. The composite-free expanded network for a collection R of regions in M is the directed graph G with vertex set R and edges given by the maintenance relation.
We denote the intersection of all regions in a subgraph G of G by I (G).
This leads to the definition of a stable module for a composite-free expanded network and to a key result about composite-free expanded networks.
Definition 5. A finite, source-free subgraph G of a composite-free expanded network G with I (G) nonempty is called a (composite-free) stable module of G with trap space I (G).
Theorem 6. The trap space, I (G), of a (composite-free) stable module G of G is closed under S, i.e., for any s ∈ S, Φ s (I (G)) ⊆ I (G) holds.
Proof. Because G is source-free, each node in G is reachable from some cycle in G. Suppose C1 → C2 → ... → CN forms a cycle in G. Consider any r ∈ I (G). Because no P (Ci, r) can be its own proper subset it must be that P (Ci, r) = S holds for all regions Ci that are part of any cycle. Because R1 → R2 implies P (R1, r) ⊆ P (R2, r), it follows that P (Ri, r) = S holds for all Ri ∈ V (G).
In order to fully align with previous work in Boolean systems, we require a formal notion of the "composite nodes" discussed in the main text. To motivate the definition of these composite nodes here, consider the following remark: Remark. If Ra i → R b i holds for finite sequences ai and bi The case in which P i Ra i , r = φ is vacuous. Hence, take P i Ra i , r non-empty.
Note that P i Ra i , r is equal to i P (Ra i , r). This implies that P i Ra i , r is a subset of P i R b i , r . Because the sequence of maintenance relations is finite, we may choose Ra k such that P (Ra k , r) ⊆ P (Ra i , r) holds for all Ra i , and thus P i Ra i , r = P (Ra k , r) holds. We can similarly choose R bm such that P i R b i , r = P (R bm , r) holds. If P (Ra k , r) = P (R bm , r), then P (Ra m , r) = P (Ra k , r) = S follows from the fact that P (Ra k , r) ⊆ P (Ra m , r) and Ra m → R bm both hold.
With this in mind, we introduce the notion of composite nodes as a shorthand for considering intersections of regions. Suppose that i R b i → R holds for some R, but that each R b i is already under separate consideration. By the previous remark, i R b i , and thus R, will be maintained if each R b i is maintained. If i R b i is considered explicitly as a region, then the intersections of regions that maintain each R b i must also be considered explicitly and separately (and the regions that maintain these must then be considered, etc.). To avoid this, we introduce the composite node C ({R b i }): With this definition, we are now able to define the expanded network as discussed in the main text: Definition 8. For a given collection R of regions in M and collection C of composite nodes with factors in R, the expanded network G is the directed graph on the nodeset V (G ) = C ∪ R with edges placed according to the following two rules: i) if R ∈ R is maintained by V ∈ V (G ), then (V, R) is an edge of G ii) if C ∈ C has a factor R ∈ R, then (R, C) is an edge in G .
Definition 9. Let G be a subnetwork of an expanded network G . G is source-free if every node of G has a parent node in G. G is consistent if the intersection of all regions, including those represented by composite nodes, in G is non-empty (i.e., I (G) is non-empty). G is complete if the factors of each composite node in G are also in G; if G is complete then G is also an expanded network. A sourcefree, consistent, and complete expanded (sub) network is called a stable module.
Corollary 10. If an expanded network G is a stable module, then I (G ) is closed under S.
Proof. Let G have composite nodes C and regions R. Let R be the set of all possible intersections of regions in R and let G be the composite-free expanded network constructed on R. Note that because G is complete and consistent, I (G) = I (G ) is non-empty. Consider any Because G is source-free, each R b i has at least one parent region Ra i for which Ra i → R b i holds, and i Ra i ∈ R maintains i R b i . Therefore, the composite-free expanded network G constructed on R is source-free and consistent. Thus G is a composite-free stable module, implying that I (G) = I (G ) is closed under S.
It is important to note that the expanded network for a dynamical system is not unique. Rather, a collection of regions and composite nodes must first be specified. After this is done, the expanded network is uniquely defined, but in practice, this network need not be fully identified. It suffices to find the edges between regions and composite nodes that form a stable module, which yields a trap space. The identification of additional edges cannot invalidate the closure of the resulting trap space under S, although these added edges may lead to the discovery of smaller trap spaces contained within the initial trap space. Therefore, trap spaces can be discovered by identifying stable modules without explicitly constructing an expanded network that the stable module is embedded within.

Worst case systems
To motivate the use of "worst-case" subsystems, as used in the algorithm presented in the main text for biological ODE systems, consider a family of dynamical systems {(T, M, Φ f ) : f ∈ F } with a shared submonoid S. We allow f to be a discrete or continuous index (i.e., F can be any set). If the maintenance relations of a stable module G hold for all members of the family, then straightforward application of the stable module theorem indicates that the corresponding trap space of the stable module is closed under each Φ S f . If a total preorder can be placed on F such that f1 f2 and Φ S f 1 (I (G)) ⊆ I (G) imply that I (G) is closed under Φ S f 2 , then it suffices to show that Φ S f min (I (G)) ⊆ I (G) holds for any fmin satisfying fmin f for all f ∈ F . This is the concept underlying the selection of a "worst-case" system. The actual system under consideration together with all externally controlled versions of that system are taken to form {(T, M, Φ f ) : f ∈ F }, and the worst-case system is one described by the monoid action Φ S f min . In our applications, F is the set of fixed regulatory effects and their values, and f1 f2 only if for all initial conditions r in a candidate stable module, Φ S f 1 (r) has a smaller minimum distance to each of the candidate thresholds than Φ S f 2 (r) does (informally, Φ f 1 is "closer to violating" the candidate stable module). Rather than explicitly constructing {(T, M, Φ f ) : f ∈ F }, we rely on the properties of monotone systems to identify the ODEs governing Φ f min directly, as is described in the main text.

Sign-Consistent Networks
Following Sontag [20], sign-consistent networks are those to which we can assign σi = ±1 to each node of the associated regulatory network such that the product σiσjis 1 whenever nodes i and j are connected by a positive edge, and −1 whenever they are connected by a negative edge. Such networks represent systems for which any effect of one variable on another is unambiguously monotonic, implying a lack of negative feedback loops and incoherent feed-forward loops. These notions are equivalent: in a graph that lacks negative cycles and feed-forward loops, starting from any node of each weakly connected component, one may assign σi = 1 and perform a breadth-first traversal ignoring edge direction, assigning subsequent σi values by alternating sign whenever a negative edge is used.
Consider the subnetwork of the regulatory network which corresponds to a candidate stable module, Sc (which, as described in the main text, is source free, compositeclosed, and consistent). The nodes of this subnetwork are the nodes characterized by the statements that make up the virtual nodes of Sc, and the edges of the subnetwork determine the maintenance or driving edges among virtual nodes of Sc. Importantly, edges among the nodes of this subnetwork that do not contribute to the stable module are not included in this subnetwork. It is not difficult to realize that a subnetwork defined this way is sign-consistent. Indeed, variables represented by virtual nodes of the form xi < T α i can be assigned σi = −1 and all others can be assigned σi = 1. This is possible because the only sources of negative regulation in the subnetwork arise from nodes of this first type, which, due to the consistency of Sc, always contribute pairs of negative edges to Sc and the associated modified regulatory subnetwork.

Example: Identification of a Stable Module
Consider the systeṁ By inspection, we note that the following structure is present in the expanded network, for unspecified thresholds: We therefore conjecture the existence of a stable module that includes this structure and (y > Ty) → (z > Tz). The only edge in the regulatory network incident to a variable represented in this candidate stable module that is not represented in this candidate stable module is the regulation of z by u. This is a negative edge because H (y) < 1, and makes the network sign-inconsistent. To search for valid thresholds for this candidate stable module, we fix the effect of this edge at the "worst case" value. In this example, the worst case corresponds to the value of u that most greatly inhibits z. Because u will be bounded by sup x = sup H (w) = 1, that value is u = 1. This yields the following modified system:ẇ Note that this is a monotone, sign-consistent system because we are holding the sign-inconsistent regulation fixed. We now find the steady states of this modified system following the MIOS procedure outlined in [20,29]. The steady states correspond to the roots of the function feedback function which is obtained by considering the regulation of (e.g.,) z by its inputs when the regulatory effect of z (on w) is held fixed (i.e., w will approach H (z), so x will approach H (H (z)), and so on until z approaches k (z) + z). There are three roots of this equation in the interval (0, 1). We consider the largest root, which we call r and find numerically to be r = 0.68, though the following discussion is equally valid with either other root used in place of r. The modified subsystem has a steady state w = H (r) , x = H (H (r)) , y = H (H (H (r))) , z = r, (5.5) and so (z > r) → (w > H (r)) → (x > H (H (r))) → (y > H (H (H (r)))) → (z > r) is a stable module (we may, by inspection of the original system, extend the module to include (x > H (r)) → (u > H (r))).
We interpret this result as follows. If w, x, y > H (r) and z > r hold at any time, then these will continue to hold for all time. This is true even if u is externally controlled, provided u remains smaller than one. If we are unsatisfied with this restriction on u, we may instead take the "worstcase" value for u to be infinity (rather than one), in which case the resulting characteristic has only one root, 1/2, in (0, 1); this corresponds to the stable motif z > 1 2 → w > 1 2 → x > 1 2 → y > 1 2 → z > 1 2 , which has a similar interpretation.

Example: Comparison of Phase-Portrait and Stable Module Analyses
We consider the system given in Equations 6.1: where H (x) = x 2 x 2 +1/16 . In this example, z = z (t) is not specified beyond being constrained to the interval [0, 1]. In this section, we will analyze this system using a traditional visual approach (phase portraits) and compare with the results obtained using the methods of the main text. To study this system using phase portraits requires a straightforward adjustment to account for the presence of the unknown function z: we examine only the projection of each vector of the vector field onto the xy-plane, with the understanding that the z-component of system trajectories varies in some unknown way. To visualize the resulting three-dimensional phase portrait, we produce slices of constant z. In Figure 6.1 we show cross-sections at equal intervals, and also at z = 0.5125, as this is approximately the value at which the number of intersections of x and y nullclines in each slice changes (these intersections do not necessarily correspond to equilibria, as z is not constant).
We can interpret this series of phase-portrait crosssections as a collection of phase portraits for systems in which z is held fixed as a constant parameter. Despite the lack of information regarding z, we can extract several conclusions from these figures. For example, we see that there is always an intersection of nullclines in which x is low (below about 0.1) and y is high (above about 0.95). One may hypothesize that x being low and y being high is robust to variations in z and test this hypothesis through careful examination of the vector field in a box bounding all of these qualitatively similar nullcline intersection points. Alternatively, one may follow our proposed methods.
The first step is to identify the causal relationships between variables. For example, if y is less than or greater than some threshold for sufficiently long, this will drive x above or below 1 − H (y). In the language of virtual nodes, this statement implies that any virtual node y ≶ T has an edge to x ≶ 1 − H (T ) in the expanded network. It is therefore reasonable to study candidate stable modules of the form y ≶ T1 ↔ x ≷ T2. Note that the direction of the inequalities is opposite for the two statements; this is a general feature of stable modules supported by mutual inhibition feedback loops. We consider that z (t) varies between 0.5125 and 1, corresponding to the region of bistability for fixed z (t) (Figure 6.1). We can analyze this case by considering the "worst-case" system for each of our two candidate stable modules (one candidate has y > T1 and x < T2, while the other has the inequalities reversed). The construction of the worst-case system relies on the observation that an increase (decrease) in z increases (decreases) y, while decreasing (increasing) x.
In one stable module candidate, x is bounded below while y is bounded above. A decrease in z perturbs x and y toward violation of these bounds. Therefore, the worst-case system for the first candidate stable module The steady state of this worst-case system with the lowest value of x (and thus the highest value of y) is given by x ≈ 0.16 and y ≈ 0.58. Therefore, one stable module for this system is x < 0.16 ↔ y > 0.58. Note that this corresponds to an intersection of nullclines in the z = 0.5125 panel of Figure 6.1. In the other candidate stable module, x is bounded below and y is bounded above. In this case, the worst-case system is characterized by a large value of z, i.e., we set z to its maximum value, one. We then find the stable module x > 0.96 ↔ y < 0.052 from the steady state of this worst-case system with the largest value of x (and thus lowest value of y). These two stable modules describe positive invariant sets (i.e., trap spaces) that are valid for as long as z remains in the interval [0.5125, 1]. We may uncover less restrictive (i.e., bigger) trap spaces by considering the other steady states of the two worst-case systems we have already constructed; any steady state of a worst case system provides valid bounds for a candidate stable module. In this section, we have used traditional techniques alongside the methods of the main text to analyze the system of Equations 6.1. We wish to emphasize that traditional visual methods, such as phase portraits and bifurcation diagrams become cumbersome as the system dimension becomes large, whereas the difficulty in solving the inequalities required for identifying thresholds for each candidate stable module scales linearly with the subsystem dimension. The crux of the method of our approach is to identify a simple "worst-case" subsystem that can be easily analyzed, and to pick out properties of the worst-case subsystem that are preserved within the real system. In this way, we can study a large, complex system by bounding the behavior of its simpler constituent subsystems. In addition, we provide the full names of the variables for convenience in