Dimension Reduction of Large AND-NOT Network Models

Boolean networks have been used successfully in modeling biological networks and provide a good framework for theoretical analysis. However, the analysis of large networks is not trivial. In order to simplify the analysis of such networks, several model reduction algorithms have been proposed; however, it is not clear if such algorithms scale well with respect to the number of nodes. The goal of this paper is to propose and implement an algorithm for the reduction of AND-NOT network models for the purpose of steady state computation. Our method of network reduction is the use of"steady state approximations"that do not change the number of steady states. Our algorithm is designed to work at the wiring diagram level without the need to evaluate or simplify Boolean functions. Also, our implementation of the algorithm takes advantage of the sparsity typical of discrete models of biological systems. The main features of our algorithm are that it works at the wiring diagram level, it runs in polynomial time, and it preserves the number of steady states. We used our results to study AND-NOT network models of gene networks and showed that our algorithm greatly simplifies steady state analysis. Furthermore, our algorithm can handle sparse AND-NOT networks with up to 1000000 nodes.


Introduction
Boolean networks (BN) have been used successfully in modeling biological networks, such as gene regulatory networks [1][2][3][4][5] and provide a good framework for theoretical analysis [6,7]. However, the analysis of large networks is not trivial. For example, even the problem of finding or counting steady states has been shown to be hard [8][9][10]. Even comprehensive sampling of the phase space is of limited use, once a model contains 50 or 100 nodes.
In order to simplify the analysis of such networks, several model reduction algorithms have been proposed [11][12][13]. However, it is not clear if such algorithms scale well with respect to the number of nodes. These reduction algorithms are based on using "steady state approximations" to remove nodes in a BN. More precisely, to remove a node i in a Boolean network, f = (f 1 , . . . , f m ) : {0, 1} n → {0, 1} n , one assumes that the i-th variable is in steady state and replaces all instances of the i-th variable by its Boolean function. For example, we can reduce the BN f (x 1 , x 2 , x 3 ) = (¬x 3 , x 1 ∨ ¬x 3 , x 1 ∧ ¬x 2 ), by making the substitution x 3 → f 3 = x 1 ∧ ¬x 2 ; then, we obtain the reduced BN h(x 1 , x 2 ) = (¬(x 1 ∧ ¬x 2 ), x 1 ∨ ¬(x 1 ∧ ¬x 2 )). This process can be repeated iteratively without changing the number of steady states.
There are two important aspects in the reduction of BNs. One is the representation of the Boolean functions (e.g. Boolean operators, polynomials, binary decision diagrams, truth tables), and the other is the way in which the reduced network is simplified to ensure that the wiring diagram is consistent with the Boolean functions (e.g. Boolean algebra, polynomial algebra, substitution). It is in these two aspects where algorithms can stop being scalable. For example, although polynomial algebra makes the manipulation of Boolean functions very systematic, the polynomial representation of simple Boolean functions can be large. For instance, storing x 1 ∨ x 2 ∨ . . . ∨ x k and ¬x 1 ∧ ¬x 2 ∧ . . . ∧ ¬x k in polynomial form grows exponentially with respect to k. On the other hand, although using Boolean operators can be more intuitive and efficient at representing Boolean functions, their simplification also grows exponentially with respect to the number of variables.
The reduction algorithm in this paper is tailored specifically to the computation of steady states of AND-NOT networks and takes advantage of the sparsity typical of gene regulatory networks. AND-NOT networks are BNs where the functions are of the form y i1 ∧ y i2 ∧ · · · ∧ y ir where y ij ∈ {x ij , ¬x ij }. We focus on AND-NOT networks because they have been shown to be "general enough" for modeling and "simple enough" for theoretical analysis [14][15][16]. Also, synthetic AND-NOT gene networks can be designed by coupling synthetic AND gates (e.g. [17]) and negative regulation. Also, AND-NOT functions are a particular case of nested canalizing functions, which have been proposed as a class of BNs for modeling biological systems [18][19][20][21][22].
Our dimension reduction algorithm for AND-NOT networks has two important properties: First, it preserves all steady state information; more precisely, there is a one-to-one correspondence between the steady states of the original and reduced network. Second, it runs in polynomial time.
As in previous reduction methods, the main idea of our algorithm is that one can use steady state approximations without changing the number of steady states; however, there are some key differences. First, the only reduction steps that are allowed are those that result in a reduced AND-NOT network. Second, since we are using AND-NOT networks only, we can make additional reductions that cannot be done with other networks. It is important to mention that AND-NOT networks are completely determined by their wiring diagrams. This is important for two reasons: First, we can store AND-NOT networks efficiently using their wiring diagrams and thus avoid the problem that the polynomial representation has. Second, we can state all reduction steps and simplification of the reduced network at the wiring diagram level and thus avoid the problem that the Boolean representation has.  given by: is an AND-NOT network. For example, f 1 = j∈{2,4} x j ∧ j∈{5} ¬x j .
For example, it is easy to check that 000011 is a steady state of the AND-NOT network in Example 2.2.
Definition 2.4. The extended wiring diagram of an AND-NOT network is defined as a signed directed graph G = (V G , E G ) with vertices V G = {0, 1, . . . , n} (or {0, x 1 , . . . , x n }) and edges E G given as follows: (i, j, +) ∈ E G ((i, j, −) ∈ E G , respectively) if x i is a positive (negative, respectively) regulator of f j . If f j = 0, then (0, j, +) ∈ E G . Positive edges are denoted by -◮ and negative edges by -•. We will refer to the extend wiring diagram as simply wiring diagram.
For example, the wiring diagram of the AND-NOT network in Example 2.2 is shown in Figure 1.

Reduction of AND-NOT Networks
3.1. Reduction Steps and Algorithm. As mentioned in the Introduction, the idea is to assume that nodes are in steady state and remove them from the network by replacing the variable by the corresponding AND-NOT function. At the wiring diagram level, the idea is to remove nodes and insert edges so that the sign of the edges are "consistent". For example, a path i-◮ j-• k should become i-• k after removing node j; and i-• j-• k should become i-◮ k after removing node j. The actual rules for doing this depend on the the properties of the node being removed and the incoming and outgoing edges. Figure 2 shows the steps at the wiring diagram level. We claim that each of these reduction steps do not change the number of steady states and that the one-to-one correspondence is algorithmic. The proofs follow directly from basic properties of Boolean algebra, so we only give the idea behind each reduction step. • Reduction Step R0. Here node i does not have any outgoing edges, so this node does not contribute to the number of steady states and can be removed. Note that given a steady state of the reduced AND-NOT network, the steady state of the original network can be found simply by inserting (in the i-th entry) x i = f i . Note that this reduction step is also valid for general BNs. • Reduction Step R1. Here we have f i = 0; and we remove node i by replacing x i with f i = 0. For example, if i-◮ j, then f j = x i ∧ w j for some AND-NOT function w j . By replacing x i with 0 we obtain f j = 0 ∧ w j = 0; that is, we add the edge 0-◮ j and remove all other incoming edges of j. On the other hand, if i-• k, then f k = ¬x i ∧ w k for some AND-NOT function w k . By replacing x i with 0 we obtain f j = ¬0 ∧ w k = w k ; that is, the edge i-• k is removed and all other edges towards k remain present. Note that given a steady state of the reduced AND-NOT network, the steady state of the original network can be found simply by inserting (in the i-th entry) x i = 0. We also notice that this reduction step is also valid for general BNs, but not at the wiring diagram level (the wiring diagram of the reduced network depends on the actual Boolean functions). • Reduction Step R2. Here we have f i = 1; and we remove node i by replacing By replacing x i with 1 we obtain f j = 1 ∧ w j = w j ; that is, the edge i-• j is removed and all other edges towards j remain present. On the other hand, if i-• k, then f k = ¬x i ∧ w k for some AND-NOT function w k . By replacing x i with 1 we obtain f j = ¬1 ∧ w k = 0; that is, we add the edge 0-◮ k and remove all other incoming edges of k. Note that given a steady state of the reduced AND-NOT network, the steady state of the original network can be found simply by inserting x i = 1. We also notice that this reduction step is also valid for general BNs, but not at the wiring diagram level. • Reduction Step R3, R4. For R3 we have f i = x j ∧ x k ∧ w i for some AND-NOT function w i , and a node j with Boolean function f j = ¬x k ∧ w j for some AND-NOT function w j . If we are at a steady state, then we have two cases, either In either case x i = 0, so by assuming that f i = 0 we are not changing the steady states of the AND-NOT network. That is, we add the edge 0-◮ i and remove all other incoming edges of i. The reduction step R4 is analogous. It is important to mention that this reduction step is not valid for general BNs. • Reduction Step R5. Here we have that f i = ¬x j ∧ ¬x k and f j = x k ∧ w j for some AND-NOT function w j . If we are at a steady state then we have two cases, either In either case we have x i = ¬x k , so by assuming that f i = ¬x k we are not changing the steady states. It is important to mention that this reduction step is not valid for general BNs. • Reduction Step R6, R7. For R6 we have that f i = ¬x k ; and we remove node i by replacing for some AND-NOT function w j2 , then we obtain f j2 = ¬x k ∧ w j2 , which is an AND-NOT function as well. Note that given a steady state of the reduced AND-NOT network, the steady state of the original network can be found simply by inserting x i = ¬x k . The reduction step R7 is analogous. We notice that this reduction step is also valid for general BNs, but not at the wiring diagram level. Also, the reduction is no longer valid if i has more incoming edges (the reduced network would not be an AND-NOT network). • Reduction Step R8. Here we have that all outgoing edges of i are positive. and we remove node i by replacing then we obtain f j = ¬x k1 ∧ x k2 ∧ w j , which is an AND-NOT function. Note that given a steady state of the reduced AND-NOT network, the steady state of the original network can be found simply by inserting x i = ¬x k1 ∧ x k2 . We also notice that this reduction step is also valid for general BNs, but not at the wiring diagram level. It is important to mention that the reduction is no longer valid if i has any negative outgoing edge. • Reduction Step R9. Here we have a circuit with positive edges only. If we are at a steady state, we have two cases, either x i1 = 0 or x i1 = 1. If x i1 = 0, then it follows that x i2 = 0 and working forward we obtain that x i1 = x i2 = . . . = 0. Similarly, if x i1 = 1, we obtain that x i1 = x i2 = . . . = 1. Thus, by collapsing this circuit into a single node we do not change the number of steady states. Note that given a steady state of the reduced AND-NOT network, the steady state of the original network can be found simply by inserting (x i2 , x i3 , x i4 , . . .) = (x i1 , x i1 , x i1 , . . .). Note that this reduction step is no longer valid if one of the edges in the circuit is negative. This reduction step is also valid for general BNs (removing one node at at time), but not at the wiring diagram level.
It is important to mention that reduction steps R0 − R8 cover the possible reductions where we only need to look at incoming and outgoing edges of a node i. Other reduction steps could be considered by looking upstream and downstream of a node; for example, one can generalize R3 and R4 to include longer feedforward loops (e.g.
However, their detection becomes computationally expensive. Reduction step R9 is included because such circuits can be detected in linear time [23].
The actual algorithm is given below. The idea is to iteratively apply the reduction steps until the network is no longer reducible (every time a reduction step is used, new reducible nodes can appear). Note that there are many orders in which one can apply the reduction steps, and in some cases they can result in different reduced networks (with the same number of states). Based on the performance of preliminary simulations, the order given below was chosen.
Algorithm. Input: AND-NOT network G. Output: List of steady states.
(3) Use R1, R2 to remove from G the nodes in Z. (4) Go to (1). (5) Use R3, R4 to find new nodes with input 0. (6) If nodes were found in previous step, then go to (1). (7) Use R5 to remove edges. (8) If there are nodes with a single incoming edge only, then use R6, R7 to remove them and go to (1). (9) Find nodes with positive outgoing edges only. (10) If nodes were found in previous step, then use R8 to remove them and go to (1). (11) Find circuits of length greater than 1 with positive edges only. Only use this step once. (12) If circuits were found in previous step, then reduce them using R9 and go to (1). (13) Compute the steady states of the reduced AND-NOT network. (14) Use the bijections given by the reduction steps to find the steady states of the original system. The algorithm has 3 main parts. In (1)-(12) we reduce the AND-NOT network; in (13) we compute the steady states of the reduced AND-NOT network; and in (14) we use these steady states to find the steady states of the initial AND-NOT network. Note that step R9 is used only once in the algorithm because none of the other steps create extra circuits.

Implementation and Computational Complexity.
We preliminarily implemented our algorithm in C++ and used the Boost Graph Library to manipulate graphs (code available upon request). We stored the one-to-one correspondence as an acyclic graph so that once the steady states of the reduced network are computed, one simply uses backward substitution to recover the steady states of the original network. The steady states of the reduced AND-NOT network are computed by exhaustive search.
Example 3.1. Consider the AND-NOT network given by: The wiring diagrams of this AND-NOT network, the reduced AND-NOT network, and the acyclic graph are in Figure 3 (see the Appendix for details about the format we use in our implementation). The reduced network is h(x 4 ) = x 4 , from which we easily obtain the steady states x 4 = 0, 1. The acyclic graph encodes the following substitution: For x 4 = 0 we obtain x 3 = 0, x 6 = 1, Since it is not known the average number of times each pattern in Figure 2 appears in a random AND-NOT network, it is difficult to predict the exact computational complexity of our algorithm. However, we present a heuristic estimation as follows. Let n be the number of nodes and e the number of edges; we denote with T the computational complexity. We first focus on steps (1)-(12) of the algorithm. In the worst case scenario steps (11) and (12) will have to be done at the beginning, this contributes O(n + e) to T [23]. Each detection of an individual pattern and reduction step from (1) to (10) takes constant time for each node; then, each one of these steps of the algorithm (not counting the "go to" statements) contributes O(n) to T ; and, since we have to repeat this at most n times (counting the "go to" statements), we obtain that steps (1)-(10) contribute O(n 2 ) to T . Thus, steps (1)- (12) contribute O(n + e + n 2 ) = O(n 2 ) to T . Assuming that one can check in constant time if a state in the reduced network is a steady state, step (13) contributes O(2 m ) where m is the size of the reduced AND-NOT network.
Note that the reduction part and backward substitution part contribute O(n 2 ) to T ; that is, they run in polynomial time.
Although T = O(n 2 + 2 m ) is of little improvement if m ≈ n, Boolean models of biological systems are not arbitrary and have especial properties. For example, they are sparse and have motifs such as feedforward loops (which we considered in R3−R5); also, they have few steady states when compared to random networks. Hence, one can argue that for Boolean models of biological systems, the reduced network are likely to be very small. Indeed, for the three Boolean models that we consider in the next section, the size of the reduced networks, m, did not exceed ln(n). The fact that m did not exceed ln(n) is important, because for family of networks where m = O(ln(n)), we obtain that T = O(n 2 + n k ) for some k > 0. Thus, one can conjecture that under some conditions our algorithm (including steady state computation) runs in polynomial time, but a formal statement and proof of this conjecture is outside the scope of this manuscript. However, our statistical analysis in Section 4.4 supports this.

Applications
In this section we apply our reduction algorithm to three published networks and random networks, and demonstrate that it can result in a significant reduction of the network's dimension. We denote two negative (positive) edges between i and j by a bidirectional negative (positive) edge, •-• (◭-◮); if the edges have different signs we denote them by •-◮.

Th-lymphocyte Differentiation.
Here we consider an AND-NOT model for Th-cell differentiation [4,14], The wiring diagram is shown in Figure 4 (left). The state space of this model has 2 26 ≈ 6.7 × 10 7 states.
By using our algorithm we reduce this AND-NOT network to the AND-NOT network shown in Figure 4 (right), h : {0, 1} 2 → {0, 1} 2 , given by h(x 1 , x 22 ) = (x 1 ∧¬x 22 , ¬x 1 ∧x 22 ). Notice that its state space has only 4 states, which is about 7 orders of magnitude smaller than the original state space. Since this is a small network, it is easy to find its steady states: 00, 01, and 10. Therefore, our results guarantee that the original AND-NOT model has 3 steady states which can be recovered from the steady states of the reduced network. The timing of our implementation was .00273562s (average of 1000000 repetitions).

ERBB2 Activation.
Here we consider an AND-NOT network model of ERBB2 activation based on the Boolean model in [24] (left). The wiring diagram of the equivalent AND-NOT model is shown in Figure 5. The state space of this model has 2 24 ≈ 1.6 × 10 7 states.  15 , f 24 = x 13 ∧ x 14 ∧ x 15 By using our algorithm we reduce this AND-NOT network to the AND-NOT network shown in Figure 5 (right), h : {0, 1} 1 → {0, 1} 1 given by h(x 20 ) = 1. Notice that its state space has only 2 states, which is about 7 orders of magnitude smaller than the original state space. It is easy to see that this network has a unique steady state (x 20 = 1). Therefore, the original AND-NOT model also has a unique steady state, which can be recovered from the steady state of the reduced network. The timing of the reduction is .00266462s (average of 1000000 repetitions).

T-cell receptor.
Here we consider an AND-NOT network model of the T-cell receptor based on the Boolean model in [25] (left). The wiring diagram is shown in Figure 6. The state space of this model has 2 43 ≈ 8.8 × 10 12 states. exhaustive search, one can use additional tools such as polynomial algebra [29], but as mentioned before, this was not necessary for our simulations. Table 1. Timimg, t, for the reduction algorithm (in seconds). The mean and standard deviation of t are denoted by µ t and σ t , respectively. Last row: best fit polynomial t = cn k .  We can see in Figure 7 that our reduction algorithm scales well with the number of nodes. Furthermore, our algorithm can reduce networks with 1000000 nodes. We also see that for very sparse networks (i.e. large values of γ) our algorithm scales very well. As sparsity is lost (i.e. as γ decreases), our algorithm becomes less and less scalable; however, as mentioned before, the value of γ for biochemical networks is usually between 2 and 3 for which our algorithm performs well. Also, the timings look polynomial (linear on the log-log scale), especially for large γ and n. The best fit polynomial of the form t = cn k for γ = 2 was given by c ≈ 10 −8 and k ≈ 2.025; and for γ = 3 was given by c ≈ 2 × 10 −6 and k ≈ 1.197. Note: Although using the timings at n = 10 2 for the estimation of c and k would give a smaller value of k, we did not use them because the linear relationship between log(t) and log(n) seemed to start at n = 10 3 .

Discussion
Since the problem of analyzing BNs is hard for large networks, many reduction algorithms have been proposed [11][12][13]. However, it is not clear if such algorithms scale well with the size of the network. In order to optimize reduction algorithms, it is necessary to focus on specific families of BNs.
The family of AND-NOT networks has been proposed as a special family simple enough for theoretical analysis, but general enough for modeling [14][15][16]. Thus, we propose an algorithm for network reduction for the family of AND-NOT networks. A key property of our algorithm is that it preserves steady states, so it can be very useful in steady state analysis. We applied our algorithm to three AND-NOT network models, namely, Th-cell differentiation, ERBB2 activation, and T-cell receptor. Our reduction algorithm performed very well with these models; the state space of the reduced networks were several orders of magnitude smaller than the original state space. This greatly simplified steady state computation. Using random AND-NOT networks, we showed that our algorithm scales well with the number of nodes and can handle large sparse AND-NOT networks with up to 1000000 nodes. To the best of our knowledge, no other algorithm can handle AND-NOT networks or any other class of (nonlinear) BNs of this size.
It is important to mention that since our reduction algorithm is defined using the wiring diagram only, it has the special property that it runs in polynomial time. That is, we have developed a polynomial-time algorithm that reduces the problem of finding steady states of an AND-NOT network, f : {0, 1} n → {0, 1} n , into the problem of finding the steady states of a smaller AND-NOT network, g : {0, 1} m → {0, 1} m , where m ≤ n. Also, the steady states of the reduced AND-NOT network can be used to compute the steady states of the original AND-NOT network in polynomial time. Thus, our algorithm transforms an NP-complete problem of input size n into an NP-complete problem of input size m, where m ≤ n. While this represents no theoretical improvement in the complexity of the problem, it does represent a significant improvement in the practical ability to analyze actual models that arise in molecular systems biology because they are very sparse, and in that case we typically have m << n. Also, this could provide a novel way to solve NP-complete problems by first using our polynomial-time algorithm as a pre-processing step.