Attractor detection and enumeration algorithms for Boolean networks

The Boolean network (BN) is a mathematical model used to represent various biological processes such as gene regulatory networks. The state of a BN is determined from the previous state and eventually reaches a stable state called an attractor. Due to its significance for elucidating the whole system, extensive studies have been conducted on analysis of attractors. However, the problem of detecting an attractor from a given BN has been shown to be NP-hard, and for general BNs, the time complexity of most existing algorithms is not guaranteed to be less than O(2n). Therefore, the computational difficulty of attractor detection has been a big obstacle for analysis of BNs. This review highlights singleton/periodic attractor detection algorithms that have guaranteed computational complexities less than O(2n) time for particular classes of BNs under synchronous update in which the maximum indegree is limited to a constant, each Boolean function is AND or OR of literals, or each Boolean function is given as a nested canalyzing function. We also briefly review practically efficient algorithms for the problem.


Introduction
The Boolean network (BN) [1][2][3] is a mathematical model that is used to represent various biological processes such as gene regulatory networks [4][5][6], neural networks [7], cell cycle control networks [8], and signal transduction networks [9][10][11][12]. For example, when modeling a gene regulatory network, each node corresponds to a gene and is assigned a Boolean value of 0 (FALSE) or 1 (TRUE), which means the gene is inactive or active, respectively. The state of each node is updated according to a Boolean function assigned to it, i.e., the edges between nodes correspond to the regulation relationships between the nodes. In addition, the state of the entire system is called a global state, which eventually reaches one of two types of stable states: singleton attractor (point attractor) or periodic attractor (cyclic attractor) by repeating the update of the state. Once the global state reaches a singleton attractor, it will always remain in the same state, whereas a periodic attractor consists of multiple global states that are cyclically visited. There is a close relationship between attractors and biological processes. For example, in the case of a gene regulatory network, an attractor is a stable state of the entire system and is considered to be a stable state of a cell. Therefore, one attractor is often regarded as one cell type [3,13,14], and extensive studies have been conducted on attractor detection, due to its significance for elucidating the system. Although it might be pointed out that a BN is a fairly simplified model and it deviates from the actual biological systems, it is meaningful to start the analysis with a simple model in order to decipher a complex system.
A BN of n nodes has 2 n global states. Thus, if we want to detect or enumerate the attractors of a given BN, we can simply update all global states as initial states and find out which initial state reaches which stable state. If n is sufficiently small, this approach is not problematic. When n is large, however, a large amount of computational resources is required. Therefore, algorithms for biological BNs have so far been limited to small networks in most cases [15]. Furthermore, from a theoretical view point, detection of a singleton attractor and enumeration of singleton attractors have been shown to be NP-hard [16,17] and #P-hard [16], respectively. In fact, there is no existing method that works in o 2 n À Á time (i.e., provably faster than in O 2 n À Á time) for general BNs. Therefore, algorithms that work efficiently for limited classes of BNs have been proposed [18,19]. For instance, there are some algorithms that detect a singleton attractor in polynomial time for very restricted classes of BNs [20,21] and several o 2 n À Á time algorithms for reasonably wide classes of BNs [6]. On the other hand, no o 2 n À Á time algorithm is known for detection of a periodic attractor of period 3 or more for a reasonably wide class of BNs. This may be because detection of an attractor with a long period has been suggested to be PSPACE-hard [6], which is more difficult than problems belonging to class NP unless NP = PSPACE. This theoretical observation also suggests a practical difficulty of detection of a long periodic attractor because PSPACE-hard problems may not be efficiently solved by using practical solvers for NP-hard problems (e.g., SAT-solver, ILP-solver).
In this paper, we review o 2 n À Á time algorithms for singleton/periodic attractor detection for particular classes of BNs under synchronous update in which the maximum indegree is limited to a constant, each Boolean function is AND or OR of literals, or each Boolean function is given as a nested canalyzing function. The reasons why the time complexities of these special classes of BNs have been well-discussed and we focus on them in this paper are (i) each Boolean function can be represented and evaluated in constant space and constant time by limiting the maximum indegree, (ii) it is possible to develop o 2 n À Á time algorithms by introducing some constraints on types of Boolean functions even when there is no constraint on the maximum indgree, (iii) biologically important functions often have a nested canalyzing form [22][23][24][25][26], and (iv) there are some other BN classes (e.g., monotonic BNs [27][28][29], monomials BNs [30,31], semilattice BNs [32], BNs whose functions consist of only XOR operations [33], BNs whose functions consist of only AND and NOT operations [34]) but algorithms guaranteed to work in o 2 n À Á time for the classes have not been proposed as far as we know. We also briefly review practical algorithms for detection and enumeration of attractors for BNs.

Boolean network and attractors
. . . ; x n f g and a set of Boolean functions F ¼ f 1 ; f 2 ; . . . ; f n f g . Each x i takes a value of 0 (FALSE) or 1 (TRUE), and in the case of a gene regulatory network, x i corresponds to the gene i, and 0 and 1 indicate that the gene is inactive or active, respectively. In particular, when writing x i t ð Þ, it indicates the state of the gene i at time t, and the global state of . . . ; x i k À Á is assigned to each x i and it indicates a rule that x i is controlled by specific input nodes be a set of input nodes for x i . The number of input nodes is denoted by jIN x i ð Þj and is called indegree of x i . We use K to denote the maximum indegree for a BN. The state of x i at time t þ 1 is determined from the state of IN x i ð Þ at time t according to the corresponding Boolean function . This may be simplified and written as Þ if there is no confusion. Furthermore, the global state at time t þ 1 can be written as The network structure of BN is represented by As an example, Fig. 1(a) shows the following BN and the corresponding graph: Fig. 1. Example of (a) Boolean network and (b) its state transition diagram.
T. Mori and T. Akutsu Computational and Structural Biotechnology Journal 20 (2022) 2512-2520 where x^y and x denote conjunction (AND) of x and y, and negation (NOT) of x, respectively. Since this BN has three nodes, the number of global states is 2 3 in total. At each time step, the transition of each global state is computed by Þ and a directed graph representing such transition of the states is called a state transition diagram ( Fig. 1(b)). As we can see from this state transition diagram, the global state at time t þ 1 is uniquely determined by the global state at time t. For example, if x t ð Þ ¼ 1; 0; 1 ½ , the global state will be 0; 0; 1 ½ at time t þ 1 and then 0; 0; 0 ½ at time t þ 2. However, once it reaches 0; 0; 0 ½ , it remains 0; 0; 0 ½ no matter how many updates are repeated. This state is called a singleton attractor or point attractor. Meanwhile, the states 1; 0; 0 ½ and 0; 1; 1 ½ repeat the transition between these states and are called a periodic attractor or cyclic attractor. More generally, an attractor is represented by a set of states Þ and a 0 ¼ f a pÀ1 À Á hold. If p ¼ 1, it corresponds to a singleton attractor. If p > 1, it corresponds to a periodic attractor and p is called the period of the attractor (the more detailed reviews on the attractor can be found in [35]). Since an attractor can be regarded as a stable state of the entire system, detection and enumeration of attractors are important analysis for detailed understanding of the system. Therefore, extensive studies have been done on the distribution and length of attractors. In particular, the NK model of BN is often focused on, which consists of N nodes and each node has K input nodes randomly selected and is assigned a Boolean function randomly selected from all possible 2 2 K functions. In the analysis of this NK model, it is well known that the expected number of singleton attractors is 1 for every K > 0 and n P K [36,37]. As for the length and distribution of periodic attractors in the NK model, several theoretical studies have been done [38,39].
There exist mainly two types of BNs, synchronous BNs whose states of all nodes are updated all at once and asynchronous BNs whose states of nodes are updated asynchronously such as random order asynchronous, general asynchronous, and deterministic asynchronous update models [11]. Asynchronous BNs are often employed since the update strategy is reasonable to take into account biological processes which occur in different time scales [6,40]. However, each state of the state transition diagram of an asynchronous BN may have more than one outgoing edges. Therefore, an attractor in an asynchronous BN is often defined as a strong connected component (SCC) without outgoing edges (referred to as bottom SCC or terminal SCC), and the attractor detection is considered as the bottom SCC detection [41]. Such a complex structure of attractor prevents us from applying SAT-based algorithm (described in Section 3.1.1) efficiently since each SAT formula becomes large [42]. Meanwhile, the fact that the synchronous model and the general asynchronous model have the same set of singleton attractors can be shown [6]. Although we deal with synchronous BNs in this paper unless otherwise stated since extensive studies on attractor detection have been done due to their simplicity, some practical attractor detection algorithms for asynchronous BNs are described in Section 5. In addition to synchronous and asynchronous BNs, another semantics on BNs called most permissive has been proposed and studied [43]. It was shown in [43] under this semantics that singleton attractor detection (i.e., detection of a fixed-point) remains NP-complete but reachability and attractor membership problems can be solved better than in PSPACE (under reasonable assumptions on the complexity classes).

Detection and enumeration of singleton attractors
One of the simplest methods for detecting or enumerating attractors is to repeat the state transitions with all possible global states as initial states for a given BN and to find out which state each initial state finally reaches. This task can be done by constructing and analyzing the state transition diagram. However, the state transition diagram for a BN with n nodes consists of 2 n nodes and 2 n edges. Therefore, if n is large, it is impossible to construct and analyze the state transition diagram. In fact, the singleton attractor detection problem has been shown to be an NP-hard problem by using a reduction from 3-SAT (Boolean satisfiability problem) [6,17]. Thus, various algorithms that work efficiently for specific classes of BNs have been proposed. In this section, we introduce singleton attractor detection/enumeration algorithms for BNs with maximum indegree K, BNs in which each Boolean function is AND or OR of literals, and BNs configured by nested canalyzing functions.

Singleton attractor detection for BNs with maximum indegree K
This section deals with BNs with maximum indegree K, where K is any constant larger than 1. Since the number of states of K input nodes is 2 K and both 0 and 1 can be assigned as an output value for each of such states, there exist 2 2 K Boolean functions with K inputs. That is, if K is a constant, each Boolean function is specified by a constant space and can be evaluated in constant time.

Reduction to SAT
We have already mentioned that SAT is used to show the NPhardness of the attractor detection problem, but it can also be applied to the singleton attractor detection itself [44][45][46]. More specifically, the singleton attractor detection for BNs with maximum indegree K can be replaced by the K þ 1 ð Þ -SAT of n variables. SAT is a problem of determining whether there is a 0 À 1 vector (assignment) a that satisfies f a ð Þ ¼ 1 when a Boolean function f is given by the following conjunctive normal form (CNF), where ' i;j is the jth literal of the ith clause. Note that ' is called a literal if it is a Boolean variable or its negation, and OR of literals is called a clause. Especially when each clause consists of at most k literals, it is called k-SAT. For example, for a BN with K ¼ 2, we consider the following transformations for the two-variable Boolean operations,^(conjunction, AND), _ (disjunction, OR), and È (exclusive OR, XOR) (the cases of constant and unary functions are omitted since they are trivial).
where, in a singleton attractor, x i t þ 1 ð Þ¼x i t ð Þ holds for all i ¼ 1; . . . ; n, so t is omitted to express the Boolean variables. This means that each Boolean function can be converted to a 3-SAT consisting of at most four clauses. Furthermore, by applying this transformation to all Boolean functions assigned to x i i ¼ 1; . . . :; n ð Þand then combining them with^, the singleton attractor detection problem of a BN with n nodes and K ¼ 2 can be converted into a 3-SAT problem consisting of n variables and at most 4n clauses. In general, the singleton attractor detection for a BN with n variables and the maximum indegree K can be converted in polynomial time to K þ 1 ð Þ -SAT with at most 2 Kþ1 Á n clauses and n variables. Although readers may think it is strange to use reduction to K þ 1 ð Þ -SAT (not to K-SAT), it is reasonable because 2-SAT is solvable in polynomial time whereas singleton attractor detection for K ¼ 2 is NP-hard [6]. To actually find a singleton attractor, we need to apply an efficient SAT solver after converting to SAT by the above method. Although SAT is an NP-hard problem, various practical solvers were developed which can solve large-scale SAT instances [47]. Furthermore, various theoretically efficient algorithms have also been developed. For example, an O 1:3303 n À Á time algorithm is known for 3-SAT [48], which implies that the singleton attractor detection problem can be solved in O 1:3303 n Á poly n ð Þ À Á time for BNs with K ¼ 2, where poly n ð Þ means some polynomial function of n.

Using feedback vertex set
The singleton attractor detection can also be done by using the feedback vertex set (FVS) in graph theory [16]. The FVS is a subset of nodes U # V in a directed graph G V; E ð Þ such that removal of all incoming edges to U eliminates all directed cycles (Fig. 2). In particular, the FVS with the smallest number of nodes is called the minimum FVS. If the state of each node v 2 U included in FVS is fixed, the state of IN v ð Þ is irrelevant to the update of v, so that, the state of FVS is propagated to all nodes except FVS in at most n À 1 steps in a graph G with n nodes. In other words, by fixing the state of FVS of a directed graph G V; E ð Þcorresponding to a given BN N V; F ð Þ, the BN reaches a stable state in at most n À 1 steps. However, it should be noted that not all stable states correspond to singleton attractors.
For sible states of FVS is 2 jUj . It is known that the problem of finding an FVS is NP-hard, but some practical algorithms have been proposed [49][50][51], and they can be used for finding singleton attractors. The FVS-based algorithm has been extended for enumeration of periodic attractors in BNs [6] and for more general non-linear models of biological networks [52]. In addition, Mori and Mochizuki also reported the relationship between the expected number of singleton attractors and the feedback arc sets [53].

Recursive algorithms
Here, we introduce a simple recursive algorithm having a guaranteed average-case complexity for singleton attractor enumeration [54]. This recursive algorithm considers a partial global state x 1 ; . . . ; x m ½ for m < n, and iteratively examines whether there is a contradiction between the partial global states and f i x ð Þ assigned to x i while incrementing m. If there is no contradiction between all x i and f i x ð Þ, output it as a singleton attractor, otherwise stop examining the current assignment and proceed to the next assignment. For example, consider the BN shown in Fig. 1(a) as an example and give x ¼ ; ; ½ as an initial state, where represents a state that has not yet been determined. Then, we examine the par- . At this time, x 2 and x 3 are undetermined, thus f x ð Þ is consistent. Therefore, we expand the subset of nodes to be assigned and examine x ¼ 0; 0; ½ . Since there is no contradiction in f x ð Þ as well, we further expand the subset and examine x ¼ 0; 0; 0 ½ . Here, f x ð Þ is consistent and the states of all nodes are determined. Therefore, x ¼ 0; 0; 0 ½ is reported as a singleton attractor. In the next step, x is not a singleton attractor and nothing is output. Then, we go back to the previous recursion level and examine . In this case, although x 3 is not yet determined, for K ¼ 2 and K ¼ 10, respectively. Therefore, when K is small, it is sufficiently faster than O 2 n À Á . In this algorithm, the partial global states are determined in a given order of x 1 ; . . . ; x n . However, if the nodes are sorted in advance in the descending order of their outdegrees, the number of nodes whose states are determined at each recursive step will increase, so that it is expected that the singleton attractor can be detected with a smaller number of trials. Actually, the algorithm based on this idea has been shown to work in O 1:19 n À Á time and O 1:57 n À Á time for K ¼ 2 and K ¼ 10, respectively. Furthermore, this algorithm can be extended for enumeration of periodic attractors with a small fixed period p [54].

Singleton attractor detection for AND/OR BNs
This section considers AND/OR BNs, in which each Boolean function is AND or OR of literals. For singleton attractor detection for AND/OR BNs, there is an algorithm that works in O 1:792 n À Á time, which is a combination of recursive calls and existing methods for solving SAT [55]. To explain the algorithm, first we assume that the following BN function is assigned for x i .
In this case, the consistent 0 À 1 assignments for x 1 ; x i ½ are 0; 0 ½ ; 1; 0 ½ , and 1; 1 ½ , whereas 0; 1 ½ is inconsistent since x i must be 0 if x 1 is 0, so that the assignment for x 1 ; x i ½ ¼ 0; 1 ½ does not need to be considered anymore. Thus a singleton attractor may be obtained by examining the other three assignments and repeating the same procedure for the other nodes. Since only 3 among 2 2 ¼ 4 assignments are examined per two nodes, this might lead to an O 3 n=2 Á poly n ð Þ ¼ O 1:733 n Á poly n ð Þ À Á time algorithm. However, this procedure cannot be continued if there are no remaining edges or only nodes with self-loops. Let U be the set of nodes whose 0 À 1 assignments are already determined by this procedure. In such case, the assignments of the set of remaining nodes W (i.e., W ¼ V n U) are determined by the following procedure.
if jUj > an, examine all possible assignments on W, otherwise, compute consistent assignments on W using a SAT algorithm.  [57]. It is to be noted that these results might be slightly improved by using a recent O 1:223 m Á poly n ð Þ À Á time algorithm for SAT with m clauses [58].

Singleton attractor detection for nested canalyzing BNs
A Boolean function f v (assigned to a node v) is called nested canalyzing if it has the following form: . . . ; x n f g and 1 6 k 1 < k 2 < Á Á Á. A BN is called an nc-BN if nested canalyzing functions are assigned to all nodes. It should be noted that both AND and OR functions of literals are special cases of nested canalyzing functions. Since biologically important functions often have a nested canalyzing form [22][23][24][25][26], it is meaningful to consider attractor detection algorithms for nc-BNs.

Recursive algorithm for singleton attractor detection for nc-BNs
Here, we briefly introduce the basic idea of a recursive algorithm SattNC for singleton attractor detection for nc-BNs [19].
The first OR part of f v (i.e., ' 1 _ Á Á Á _ ' k 1 À1 ) is called the initial clause. Suppose that u appears positively in the initial clause of some f v where u -v. If 1 is assigned to u, then we have f v ¼ 1 and thus the state of v becomes 1 regardless of the state of other input nodes to v. This means that assigning 1 to v determines the states of two nodes (u and v). On the other hand, if 0 is assigned to u, the state of v may not be determined. Conversely, suppose that u appears negatively in f v (i.e., f v ¼ u _ f v 0). In this case, u ¼ 0 determines the states of two nodes (u and v), whereas u ¼ 1 determines the state of one node (u). This procedure may be applied recursively. Let G q ð Þ be the number of recursive calls with q unassigned variables in such a recursive procedure. Then, we may have G q ð Þ 6 G q À 1 ð Þþ G q À 2 ð Þ from the above discussion. If we can repeat this procedure until the states of all nodes are determined, G q ð Þ will be O 1:619 q À Á (the Fibonnaci number). However, in many cases, this procedure cannot be repeated so many times. Furthermore, assigning u ¼ 1 (resp., u ¼ 0) gives a constraint f u ¼ 1 (resp., f u ¼ 0). In order to cope with these issues, we need to solve SAT for a set of nested canalyzing functions, which is more general than the standard SAT for clauses (i.e., a set of OR of literals). SattNC was obtained by combining the above recursive procedure with a newly developed SAT algorithm for nested canalyzing functions, and was shown to work in O 1:871 n À Á time [19].

Singleton attractor detection for nc-BNs with bounded treewidth
It is known that many NP-hard problems can be solved in polynomial time using dynamic programming if an input graph has a fixed treewidth, where the treewidth is an integer value measuring how close a graph is to a tree [59,60]. To formally define the treewidth, we consider the tree decomposition that transforms called a bag, each node v i 2 V must appear in at least one B t , nodes in each edge e 2 E must be included in at least one B t , and B t 's containing each node must be connected in T (see Fig. 3). Note that the tree decomposition is not uniquely determined. The width of the decomposition is defined by max t2VT jB t j À 1 ð Þ and the treewidth of G is the smallest width among all tree decompositions of G.
It is known that the singleton attractor detection problem can be solved in O n 2 wþ1 ð Þ Á poly n ð Þ À Á if a given BN N V; F ð Þ is composed of nested canalyzing functions and the treewidth of G V; E ð Þ is bounded by w [61]. This algorithm applies dynamic programming to T V T ; E T ð Þand computes partial assignments from leaves to the root in a bottom-up manner. Chang et al. have improved this algorithm for the special cases of AND/OR BNs and nc-BNs [62]. There is also an algorithm that reduces singleton attractor detection for nc-BNs with bounded treewidth to a constraint satisfaction problem for bounded treewidth [61,63].

Detection of periodic attractors
As mentioned in Section 1, detection of a periodic attractor seems much harder than detection of a singleton attractor. Nevertheless, some o 2 n À Á time algorithms have been developed for detection of periodic attractors with a short period. This section briefly introduces such algorithms.

Reduction to the singleton attractor detection
A simple strategy for detecting periodic attractors is to reduce it to the singleton attractor detection problem. Here, we introduce a simple reduction algorithm that constructs p-multiplied network N p V p ; F p ð Þ from an input BN N V; F ð Þ [61]. N p V p ; F p ð Þ is defined as follows:  Fig. 1  (a).
where each x t ð Þ i is regarded as a node and f t ð Þ i ffi f i means that f t ð Þ i is the same Boolean function as f i except that input variables are different (see Fig. 4). A singleton attractor at N p V p ; F p ð Þclearly corresponds to an attractor of N V; F ð Þ with period q that divides p. However, in order to guarantee that the detected attractor has exactly period p, for all q ¼ 2; . . . ; p, x This means that the state of N V; F ð Þat time t ¼ 1 must be different from that of at time t ¼ 2; . . . ; p, so that N V; F ð Þdoes not take the same global state at t ¼ q 1 and t ¼ q 2 1 6 q 1 < q 2 6 p ð Þ . Here, let / t ð Þ be a function from 2; . . . ; p f gto 1; . . . ; n f gand w t ð Þ be a function from 2; . . . ; p f gto 0; 1 f g. To check if there is a periodic attractor with period p of N V; F ð Þ, it is enough to check whether there exist an singleton attractor of N p V p ; F p ð Þand functions / t ð Þ and w t ð Þ such that

2-periodic attractor detection for AND/OR BNs
This section describes a 2-periodic attractor detection algorithm for AND/OR BNs using N 2 V 2 ; F 2 [61]. First, we transform an AND/ OR BN to an OR BN. Let x i be an AND node assigned the following function: This Boolean function can be transformed to an OR function by replacing it to and negating all x i t þ 1 ð Þin the function f j for all x j 2 V. Therefore, in this section, we can assume that N 2 V 2 ; F 2 is an OR BN.
The basic strategy of the 2-periodic attractor detection algorithm is similar to that of the O 1:587 n À Á time algorithm for an AND/OR BN, that is, we recursively examine 0 À 1 assignments and finally apply an SAT algorithm. However, here we use the special property on N 2 V 2 ; F 2 . The procedure of 2-periodic attractor detection is to first construct N 2 V 2 ; F 2 and set all nodes to be unassigned. Then, we recursively examine 0 À 1 assignments for unassigned nodes x such that U x ð Þ P 3 until no such nodes exist or the number of assigned nodes is more than H, where U x ð Þ is the number of unassigned neighboring nodes of x and H is a parameter. Here, let A be the set of assigned nodes, and let A 1 ¼ A \ V 2 1 and A 2 ¼ A \ V 2 2 (w.l.o.g., jA 1 j P jA 2 j), where V 2 1 and V 2 2 are the sets of nodes V 2 corresponding to t ¼ 1 and t ¼ 2, respectively. If jAj > H, we examine all possible assignments for V 1 n A 1 , otherwise recursively examine assignments for paths and cycles (because all nodes with U x ð Þ P 3 have already been assigned), and then finally solve SAT. The resulting algorithm has been shown to work in O 1:985 n À Á time by letting H ¼ 0:3196n [61]. It is known that detection of 2-periodic attractor can be done in linear time for a positive OR BN in which all Boolean functions are OR functions and all variables appear positively [21,61]. However, periodic attractor detection remains NP-hard for a positive BN in which each function is an AND or OR function [64].

Periodic attractor detection for nc-BNs with bounded treewidth
The algorithm presented in Section 3.3.2 can be extended for detection of a p-periodic attractor for an nc-BN with bounded treewidth, using the p-multiplied network N p V p ; F p ð Þ [61]. The extended algorithm is based on the following proposition: if the graph G V; E ð Þ associated with N V; F ð Þ has the treewidth w, the graph G p V p ; F p ð Þassociated with N p V p ; F p ð Þhas a tree decomposition with the width less than p w þ 1 ð Þ , and for each x 2 V, x 1 ð Þ ; . . . ; x p ð Þ are included in the same B t [61]. Then, we can apply a dynamic programming procedure to the tree decomposition of G p V p ; F p ð Þ associated with N p V p ; F p ð Þ , where some modifications are needed to detect an attractor with exactly period p. The resulting algorithm works in O n 2p wþ1 ð Þ Á poly n ð Þ À Á time [61]. Furthermore, the algorithm works in O g d; p; w ð ÞÁpoly n ð Þ ð Þ time if the maximum indegree is bounded by a constant d where g is a non-polynomial function depending on d; p; w but not depending on n. It gives a fixed-parameter algorithm when p; w, and d are parameters (i.e., the exponential factor of the time complexity depends only on some parameters for the algorithm) [59,60].

Practical algorithms
In the previous sections, we have introduced singleton/periodic attractor detection and enumeration algorithms that are guaranteed to work in less than O 2 n À Á time. However, these algorithms are limited to particular subclasses of BNs and no such an algorithm has been known for general BNs. On the other hand, there are many algorithms that do not have such a theoretical guarantee but work efficiently in practice.

SAT and logic-based approaches
One major practical approach is to use SAT and related logicbased methods. As mentioned in Section 3.1.1, attractor detection problems can be reduced to SAT. Furthermore, many practically efficient SAT solvers have been developed [44,45]. Therefore, it is reasonable to develop practically efficient attractor detection/enumeration methods for BNs and related models using such SAT solvers. For example, Dubrova and Teslenko used a SAT solver to search for a path of length p on the state transition diagram and detected an attractor by checking if it contains a loop [45], and de Jong and Page transformed the problem of searching for a stable state of a network described by the pairwise-linear differentiation equation model into SAT [44]. For other logic-based approaches, Devloo et al. developed an algorithm applying constraint programming to the singleton attractor detection and enumeration [65]. Inoue provided an algorithm that directly encodes a BN into a logic program and computes a singleton attractor based on that logic program [66], and Abdallah et al. proposed an algorithm based on answer set programming (ASP) that enumerates all attractors without creating an entire state transition diagram [67]. In addition to SAT, Binary decision diagrams (BDDs) have also been utilized to solve large scale logic-based problems in various fields. Therefore, it is reasonable to apply BDDs in place of SAT solvers. Indeed, some methods have been developed based on BDDs that can detect attractors for large networks [68][69][70]. Integer linear programming (ILP) is another useful method to efficiently solve Boolean constraints and thus has been applied to attractor detection and related problems on BNs [71][72][73].

Network reduction-based approaches
Another major practical approach is to use network reduction.
For example, suppose that v has an input node u and u has only one input node w. Then, for both detection and enumeration of singleton attractors, we can remove u by letting w as an input to v.
Although this example is a very simple one, various reduction methods were developed where logic-based or other methods were used to solve reduced instances at the final stage. For example, Veliz-Cuba et al. proposed various reduction rules including the above mentioned one, elimination of redundant edges, replacement of functions, and simplification of Boolean functions [74]. Veliz-Cuba et al. developed a singleton attractor enumeration algorithm that transforms input networks to AND-NOT networks, then applies network reduction, and finally applies an algebraic method (compute the Gröbner basis to perform a generalized version of Gaussian elimination) [75]. He et al. developed a reduction method based on removal of unstable partial states and identification of constant nodes [76]. Saadatpour et al. provided rigorous proofs of conservation of attractors in synchronous and asynchronous BNs for various reduction methods [77]. Beneš et al. proposed another network reduction method called interleaved transition guided reduction (ITGR) for detecting of bottom SCCs of asynchronous BNs [78]. They succeeded in handling large Boolean networks of real data (up to 350 variables) as well as of synthetic data (up to 1100 variables) based on the reduction technique. Gao et al. characterized periodic attractors of a conjunctive BN (i.e., each Boolean function is AND of literals) whose underlying directed graph is strongly connected by establishing bijection between the set of periodic attractors and the set of binary necklaces (i.e., character strings over the binary set 0; 1 f g, where their all rotations are dealt with as equivalent) of a certain length [79]. Chen et al. extended this study to BNs over weakly connected directed graph by applying network reduction [80]. It is to be noted that network reduction methods were utilized also in some other methods mentioned in this section.

Network decomposition-based approaches
Divide-and-conquer is one of the major general techniques to efficiently solve various combinatorial problems in computer science. Therefore, it is reasonable to apply divide-and-conquer to the attractor detection/enumeration problems on BNs. Indeed, various methods have been developed for the problem by decomposing the original network into subnetworks and then reconstructing global attractors by integrating local attractors for the subnetworks, where the techniques introduced in Sections 5.1 and 5.2 are also utilized in many of them. As a pioneering work in this direction, Irons proposed an algorithm that combines partial state and predecessors [81]. Mizera et al. developed another algorithm by using decomposition of the original BN into stronglyconnected components (SCCs) [41], where an SCC is a well-known concept in graph theory and is a maximal subnetwork in which all node pairs are connected by directed paths. Su et al. significantly improved this algorithm by further partitioning SCCs [82]. Zañudo and Albert introduced the concept of a stable motif, which is an SCC that stabilizes in attractors and made use of stable motifs to efficiently find attractors in the whole BN [83]. Klarner et al. introduced a similar concept, trap space, and combined it with ILP [72] and model checking [84]. Choo and Cho developed a method based on hierarchically partitioning with focusing on attractors corresponding to a particular phenotype of interest instead of considering all attractors [85]. Tamaki developed an algorithm combining path-decomposition and partial states [86]. It should be noted that the tree decomposition-based methods mentioned in Sections 3.3.2 and 4.3 can also be regarded as network decomposition-based ones.

FVS-based approaches
The FVS-based singleton attractor detections for synchronous BNs are described in Section 3.1.2. Meanwhile, practical FVSbased approaches for asynchronous random Boolean networks (ARBNs) have also been proposed [42,87]. Skodawessely and Klemm focused on a reduced dynamics to make it easier to find attractors of ARBNs [87]. A reduced dynamics is obtained by retaining a Boolean state of a node in FVS and then removing state transitions from the state transition diagram. Finally, the attractors of an ARBN can be found from the sets of fixed points of the reduced dynamics. Although this approach seems to be widely applicable, it is often difficult to handle large networks (n P 30) since it requires large memory space to store all the state vectors of an attractor during the traverse of the original state transition diagram. Van Giang et al. were inspired by the concept of reduced dynamics and presented the relations between FVSs and dynamics of BNs with formal proofs [42]. Furthermore, they proposed another FVS-based method for detecting attractors in ARBNs based on the relations. Its main idea is to systematically remove edges from the transition diagram and then filter out the set of fixed points of them to get a candidate set of states which one-to-one corresponds to the set of attractors of the given ARBN. The computational experiments using real biological networks and random N-K networks showed that the method succeeded in handling large networks whose sizes are up to 101 without any network reductions.

Other approaches
It is well-known that starting from an arbitrary state and repeating update of states, the trajectory will finally fall into a singleton or periodic attractor. Therefore, starting from many initial states, we may be able to get directed attractors or some statistics on attractor distributions. Although such sampling-based studies have been conducted on small-size BNs, it is difficult to apply such methods to analysis of large-scale BNs.
Another major approach is to use semi-tensor product (STP) [5], which is an extension of matrix product. In this approach, many problems on BNs are defined as matrix-based problems. Since matrices play central roles in control theory, various concepts and methodologies in control theory have been applied to BNs and thus many studies have been done using STP [88]. For example, a state transition diagram can be represented as a binary matrix A of size 2 n Â 2 n in which A ij ¼ 1if and only if there exists a transition from state i to state j [5]. Then, singleton attractors correspond to elements of A such that A ii ¼ 1. Therefore, once A is constructed, detection and enumeration of singleton attractors become trivial [5]. Although most of existing STP-based methods need to handle 2 n Â 2 n or larger size matrices and thus can only handle small-size BNs, some efforts have been done to address this complexity issue [89,90].

Conclusion
In this article, we briefly reviewed algorithms for singleton/periodic attractor detection/enumeration with focusing on synchronous BNs. Enumeration of all singleton and periodic attractors can trivially be done by constructing and analyzing a state transition diagram. However, such an approach needs O 2 n Á poly n ð Þ À Á or more computation time because there exist 2 n possible states. Since 2 n is often too large, it is important to develop theoretically and/or practically efficient algorithms. However, it is not an easy task because it is known that detection of a singleton attractor is NP-hard, enumeration of a singleton attractor is #Phard, and detection of a long periodic attractor is suggested to be PSPACE-hard [6].
We have seen that singleton attractor detection can be done in provably faster than in O 2 n À Á time (i.e., still exponential, but o 2 n À Á time) for some special but reasonably wide classes of BNs such as AND-OR BNs and BNs consisting of nested canalyzing functions, by combining recursive assignment techniques with SAT algorithms. However, these presented algorithms are not necessarily optimal. Therefore, improvement of the time complexity of these algorithms is left as an open problem. We have also seen that singleton attractor detection can be done in polynomial time for BNs with bounded treewidth and bounded maximum degree, using tree decomposition and dynamic programming as in many other polynomial-time algorithms for graphs with bounded treewidth. As for enumeration of singleton attractors, we have reviewed a simple recursive algorithm that works in o 2 n À Á time in the average case.
Detection of a periodic attractor is much more difficult. We have seen that detection of an attractor with period two can be done in o 2 n À Á time for AND-OR BNs. However, as far as we know, there is no o 2 n À Á time algorithm that can detect an attractor with period three (or more) for reasonably wide classes of BNs. Therefore, development of such an algorithm is left as an open problem. Recently, in order to cope with this difficulty, an algorithm has been proposed which detects a long periodic attractor in o 2 n À Á expected time, assuming that partial 0 À 1 assignments in a desired attractor is probabilistically given in advance as a priori [15]. This method succeeded in detecting attractors of large networks with long periods that were difficult to find by existing methods. From a practical viewpoint, many methods have been developed for detection/enumeration of singleton/periodic attractors for synchronous and asynchronous BNs. However, it is unclear which methods are the most useful in practice, where it may depend on structures of the target BNs. Therefore, rigorous comparison of state-of-the-art methods should be done. Furthermore, it seems that detection/enumeration of long attractors remains quite difficult also in practice. Therefore, practically efficient methods should be developed for detection/enumeration of long attractors.