Boolean network (BN) models are widely used in molecular and systems biology to capture coarse-grained dynamics of a variety of regulatory networks, with a particular focus on features such as steady state behavior [1–22]. One advantage of discrete models of this type is that, for small models, the entire dynamics can be explored by exhaustive enumeration of all state transitions. Since the size of the state space of a Boolean model with *n* nodes is 2^{
n
}, this approach becomes unfeasible for larger models, those with more than approximately 30 variables, depending on the computational resources available. Also, for larger models, finding steady states (fixed points in this manuscript) through sampling is not effective anymore either, since even large attractors can be missed entirely by this approach. On the theoretical side, it has been shown that the problem of finding, or even counting, steady states of Boolean networks is NP-hard [23, 24], so that any algorithm for this problem will eventually encounter serious limitations. Since the size of published models has increased in recent years, including models with 100 or more nodes [15, 17, 21, 22], it is important to develop more efficient methods that find all steady states of a Boolean model.

Several methods have been proposed in the literature for dealing with this problem, including exact as well as heuristic methods. We provide a brief review of the different types here. For this purpose, we represent a Boolean network as follows. Let

*K*={0,1}, and assume that the network has

*n* nodes

*x*
_{1},…,

*x*
_{
n
}. Each node

*x*
_{
i
} has associated to it a Boolean function

*f*
_{
i
}:

*K*
^{
n
}→

*K*. Thus, we can represent the Boolean network as a function

One can represent the variable dependencies through the *dependency graph* of the network, whose nodes are the variables *x*
_{1},…,*x*
_{
n
}. There is an edge *x*
_{
i
}→*x*
_{
j
} if *x*
_{
i
} appears in the function *f*
_{
j
}, that is, the state of *x*
_{
j
} depends on the state of *x*
_{
i
}. The problem of finding steady states is then formulated as finding all states *x*∈*K*
^{
n
} such that *f*(*x*)=*x*.

One approach to the problem is model reduction. Some existing *reduction methods* use a “steady-state approximation” [25–28] to reduce the number of variables. Intuitively, if a function depends on a variable, e.g., *f*
_{
i
}=*f*
_{
i
}(*x*
_{
j
},*x*
_{
k
},*x*
_{
l
}), then we can remove variable *x*
_{
j
} from the network by replacing *f*
_{
i
}(*x*
_{
j
},*x*
_{
k
},*x*
_{
l
}) with the new function *f*
_{
i
}(*f*
_{
j
}(*x*
_{1},…,*x*
_{
n
}),*x*
_{
k
},*x*
_{
l
}). By repeating this process, one obtains a reduced network that in practice is much smaller than the original network. The stopping criteria for reduction methods is that variables can be removed only if the steady state information is preserved. The steady states of the reduced network are in algorithmic one-to-one correspondence with the steady states of the original network. More precisely, the reduction algorithm decomposes a large system into a smaller system and a set of equations in triangular form, so that when the steady states of the reduced system are found, the steady states of the original systems can be found simply by backwards substitution. That is, the existence of the one-to-one correspondence is not just theoretical.

Another method uses the fact that one can represent a Boolean function as a polynomial function in the variables *x*
_{1},…,*x*
_{
n
}, with coefficients in the finite number system *K*={0,1} (with integer addition and multiplication modulo 2). The problem of finding the steady states of a Boolean network in *n* variables, as above, can then be reformulated as the problem of finding the solutions to a system of polynomial equations *p*
_{
i
}:=*f*
_{
i
}(*x*
_{1},…,*x*
_{
n
})−*x*
_{
i
}=0;*i*=1,…,*n*
[29–31]. Then, the roots of the system of polynomial equations is encoded by the set {*p*
_{1},…,*p*
_{
n
}}. Using tools from computational algebra it is possible to find another set that has the same roots (a Gröbner basis), such that it is possible to do a generalized version of Gaussian elimination. These computations can be done using several different software packages developed for this purpose.

A graph-theoretic method, *Minimal Feedback Vertex Sets*, consists of finding a set of vertices in the dependency graph of the network that “generate” all steady states. More precisely, one finds a set *S*⊂{1,…,*n*} such that knowing the coordinates *x*
_{
i
}, for all *i*∈*S*, of a steady state completely and algorithmically determines the other coordinates of the steady state. It turns out that so-called feedback vertex sets have this property. In practice, by finding a minimal feedback vertex set, one reduces the problem from checking 2^{
n
} states to the problem of checking 2^{|S|} states, where |*S*| is typically much smaller than *n*
[23]. A feedback vertex set can be found by removing vertices from the graph until the graph has no directed cycles. A minimal feedback vertex set can be found by finding the smallest number of vertices that we need to remove from the graph so that it does not have directed cycles.

SAT methods have also been used for the purpose of finding steady states of Boolean networks, which are used to determine whether a Boolean expression in several variables has a variable assignment that makes the expression true; see [32–35]. In essence, the system of Boolean equations, *f*
_{
i
}=*x*
_{
i
}, is rewritten as a single equation *G*(*x*)=1, and then the problem of finding the steady states becomes the problem of finding when the equation *G*(*x*)=1 is satisfied. For example, Melkman, Tamura, and Akutsu [33, 35] used SAT algorithms to find steady states of AND/OR Boolean networks, i.e., Boolean networks in which the *f*
_{
i
} contain only the AND and OR operators, with a time complexity of *O*(1.587^{
n
}) (where *n* is the number of nodes). Dubrova and Teslenko [34] also developed a SAT-based algorithm to find all attractors of a Boolean network with very good performance characteristics. The methodology was tested on Boolean networks with sizes ranging from 12 to 52. It was also tested using random networks with up to 7000 nodes and average in-degree less than 2. For a fixed in-degree of 2 the maximum size networks tested have 2000 nodes.

Integer programming-based method have also been used to find the steady states of Boolean networks, Tamura, Hayashida, and Akutsu [36]. In essence, the system of Boolean equations is rewritten as a set of inequalities *A*
*x*≤*b*,*x*≥0 and the goal is to maximize a linear function of the form *c*
^{
T
}
*x*.

*Strategic Sampling*, (Zhang, Hayashida, Akutsu, Ching, and Ng, [37]) is a recursive search approach to identify all steady states of a random Boolean network with maximum in-degree 2, with an average time complexity of *O*(1.19^{
n
}) (where *n* is the number of nodes). The idea is that the equations are solved recursively: First one considers the solutions of the equation *f*
_{1}=*x*
_{1}. Since the *f*
_{
i
}’s depend on few variables in practice, one only has to keep track of the variables that appear in *f*
_{1}. Then, one finds the solutions of *f*
_{2}=*x*
_{2} that are compatible with the solutions previously found. The process continues until one finds solutions of all equations. In the worst case, however, algorithm complexity can be *O*(*n*2^{
n
}) [31].

Finally, the problem of finding attractors has also been studied by using Binary Decision Diagrams (BDD) [38–41]. The idea is to represent the Boolean functions as a directed graph that efficiently encodes the functions by allowing fast evaluation. Then, by combining the BDD representation of all the Boolean functions, the problem of finding steady states becomes a search problem in the larger BDD. Many of these methods were tested on some biologically relevant networks with fewer than 100 nodes.

In this paper, we present a new method for computing steady states of a Boolean network, combining a graph theoretic reduction/transformation method with an approach using computational algebra. We show that the method performs favorably on some types of networks in comparison with other methods on a collection of benchmark networks, consisting of both published models and random networks with certain properties, namely Kauffman networks and networks whose in-degree distribution satisfies a power law.