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 ( 2 n ) . 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 ( 2 n ) 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 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 time (i.e., provably faster than in 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 time algorithms for reasonably wide classes of BNs [6]. On the other hand, no 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 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 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 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

A BN consists of a set of nodes and a set of Boolean functions . Each takes a value of 0 (FALSE) or 1 (TRUE), and in the case of a gene regulatory network, corresponds to the gene i, and 0 and 1 indicate that the gene is inactive or active, respectively. In particular, when writing , it indicates the state of the gene i at time t, and the global state of the entire system at time t is denoted by . A Boolean function is assigned to each and it indicates a rule that is controlled by specific input nodes . Here, let be a set of input nodes for . The number of input nodes is denoted by and is called indegree of . We use K to denote the maximum indegree for a BN. The state of at time is determined from the state of 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 can be written as . The network structure of BN is represented by a directed graph , where . As an example, Fig. 1(a) shows the following BN and the corresponding graph:where and 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 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 is uniquely determined by the global state at time t. For example, if , the global state will be at time and then at time . However, once it reaches , it remains no matter how many updates are repeated. This state is called a singleton attractor or point attractor. Meanwhile, the states and 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 , where and hold. If , it corresponds to a singleton attractor. If , 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 functions. In the analysis of this NK model, it is well known that the expected number of singleton attractors is 1 for every and [36], [37]. As for the length and distribution of periodic attractors in the NK model, several theoretical studies have been done [38], [39].

Fig. 1

Example of (a) Boolean network and (b) its state transition diagram.

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 nodes and 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 and both 0 and 1 can be assigned as an output value for each of such states, there exist 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 NP-hardness 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 -SAT of n variables. SAT is a problem of determining whether there is a vector (assignment) that satisfies when a Boolean function f is given by the following conjunctive normal form (CNF),where 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 , 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, holds for all , 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 and then combining them with , the singleton attractor detection problem of a BN with n nodes and can be converted into a 3-SAT problem consisting of n variables and at most 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 -SAT with at most clauses and n variables. Although readers may think it is strange to use reduction to -SAT (not to K-SAT), it is reasonable because 2-SAT is solvable in polynomial time whereas singleton attractor detection for 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 time algorithm is known for 3-SAT [48], which implies that the singleton attractor detection problem can be solved in time for BNs with , where 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 in a directed graph 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 included in FVS is fixed, the state of is irrelevant to the update of v, so that, the state of FVS is propagated to all nodes except FVS in at most steps in a graph G with n nodes. In other words, by fixing the state of FVS of a directed graph corresponding to a given BN , the BN reaches a stable state in at most steps. However, it should be noted that not all stable states correspond to singleton attractors. For example, if FVS = and , BN reaches a specific stable state by fixing , but this state is not a singleton attractor for a given BN since must be 1. Therefore, in order to enumerate singleton attractors, it is enough to assign all possible states to the minimum FVS and to update the states until the BN reaches stable states, and then to check their consistencies. Assuming that the size of minimum FVS is and each Boolean function of can be evaluated in polynomial time, enumeration of all singleton attractors can be done in time since the number of possible states of FVS is . 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].

Fig. 2

Example of a BN for illustration of FVS. is the minimum FVS. is a FVS but is not the minimum one.

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 for , and iteratively examines whether there is a contradiction between the partial global states and assigned to while incrementing m. If there is no contradiction between all and , 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 as an initial state, where represents a state that has not yet been determined. Then, we examine the partial global state of . At this time, and are undetermined, thus is consistent. Therefore, we expand the subset of nodes to be assigned and examine . Since there is no contradiction in as well, we further expand the subset and examine . Here, is consistent and the states of all nodes are determined. Therefore, is reported as a singleton attractor. In the next step, is examined. Since , is not a singleton attractor and nothing is output. Then, we go back to the previous recursion level and examine . In this case, although is not yet determined, is already established, so that the examination of is stopped and then is examined. By repeating this procedure, it is possible to list all singleton attractors. Letting , theoretical analysis shows that this algorithm works in time on the average, where the average is taken over all NK models, and . For example, the resulting time complexity is and for and , respectively. Therefore, when K is small, it is sufficiently faster than . In this algorithm, the partial global states are determined in a given order of . 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 time and time for and , 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 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 .In this case, the consistent assignments for are , and , whereas is inconsistent since must be 0 if is 0, so that the assignment for 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 assignments are examined per two nodes, this might lead to an 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 assignments are already determined by this procedure. In such case, the assignments of the set of remaining nodes W (i.e., ) are determined by the following procedure.

if , examine all possible assignments on W,

otherwise, compute consistent assignments on W using a SAT algorithm.

It is shown that by letting and utilizing an time algorithm for SAT with m clauses [56], the resulting algorithm works in time. By further omitting unnecessary examinations of partial assignments, an time algorithm was developed [57]. It is to be noted that these results might be slightly improved by using a recent time algorithm for SAT with m clauses [58].

Singleton attractor detection for nested canalyzing BNs

A Boolean function (assigned to a node v) is called nested canalyzing if it has the following form:where and . 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 for singleton attractor detection for nc-BNs [19].

The first OR part of (i.e., ) is called the initial clause. Suppose that u appears positively in the initial clause of some where . If 1 is assigned to u, then we have 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 (i.e., ). In this case, determines the states of two nodes (u and v), whereas determines the state of one node (u). This procedure may be applied recursively. Let be the number of recursive calls with q unassigned variables in such a recursive procedure. Then, we may havefrom the above discussion. If we can repeat this procedure until the states of all nodes are determined, will be (the Fibonnaci number). However, in many cases, this procedure cannot be repeated so many times. Furthermore, assigning (resp., ) gives a constraint (resp., ). 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). was obtained by combining the above recursive procedure with a newly developed SAT algorithm for nested canalyzing functions, and was shown to work in 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 to a pair , where is a rooted tree, is a family of subsets of nodes , each is called a bag, each node must appear in at least one , nodes in each edge must be included in at least one , and ’s containing each node must be connected in (see Fig. 3). Note that the tree decomposition is not uniquely determined. The width of the decomposition is defined by and the treewidth of G is the smallest width among all tree decompositions of G.

Fig. 3

Example of graph and its tree decomposition with width two, where are bags.

It is known that the singleton attractor detection problem can be solved in if a given BN is composed of nested canalyzing functions and the treewidth of is bounded by w [61]. This algorithm applies dynamic programming to 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 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 from an input BN [61]. is defined as follows:where each is regarded as a node and means that is the same Boolean function as except that input variables are different (see Fig. 4). A singleton attractor at clearly corresponds to an attractor of with period q that divides p. However, in order to guarantee that the detected attractor has exactly period p, for all , must hold for some . This means that the state of at time must be different from that of at time , so that does not take the same global state at and . Here, let be a function from to and be a function from to . To check if there is a periodic attractor with period p of , it is enough to check whether there exist an singleton attractor of and functions and such that and holds for all . The possible number of s and s are and , respectively, which are in polynomial of n for a constant p. Therefore, if singleton attractor detection can be done in time, p-periodic attractor detection canbe done in time. However, even if the time algorithm for singleton attractor detection for AND/OR BNs is applied, it does not yield an time algorithm because it would take time () to detect a periodic attractor with period 2.

Fig. 4

p-multiplied network corresponding to given in Fig. 1(a).

2-periodic attractor detection for AND/OR BNs

This section describes a 2-periodic attractor detection algorithm for AND/OR BNs using [61]. First, we transform an AND/OR BN to an OR BN. Let be an AND node assigned the following function:This Boolean function can be transformed to an OR function by replacing it toand negating all in the function for all . Therefore, in this section, we can assume that is an OR BN.

The basic strategy of the 2-periodic attractor detection algorithm is similar to that of the time algorithm for an AND/OR BN, that is, we recursively examine assignments and finally apply an SAT algorithm. However, here we use the special property on . The procedure of 2-periodic attractor detection is to first construct and set all nodes to be unassigned. Then, we recursively examine assignments for unassigned nodes x such that until no such nodes exist or the number of assigned nodes is more than H, where 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 and (w.l.o.g., ), where and are the sets of nodes corresponding to and , respectively. If , we examine all possible assignments for , otherwise recursively examine assignments for paths and cycles (because all nodes with have already been assigned), and then finally solve SAT. The resulting algorithm has been shown to work in time by letting [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 [61]. The extended algorithm is based on the following proposition: if the graph associated with has the treewidth w, the graph associated with has a tree decomposition with the width less than , and for each , are included in the same [61]. Then, we can apply a dynamic programming procedure to the tree decomposition of associated with , where some modifications are needed to detect an attractor with exactly period p. The resulting algorithm works in time [61]. Furthermore, the algorithm works in time if the maximum indegree is bounded by a constant d where g is a non-polynomial function depending on but not depending on n. It gives a fixed-parameter algorithm when , 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 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 logic-based 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 , 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, 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 strongly-connected 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, 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 FVS-based 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 () 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 in which if and only if there exists a transition from state i to state j [5]. Then, singleton attractors correspond to elements of A such that . Therefore, once A is constructed, detection and enumeration of singleton attractors become trivial [5]. Although most of existing STP-based methods need to handle 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 or more computation time because there exist possible states. Since 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 #P-hard, 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 time (i.e., still exponential, but 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 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 time for AND-OR BNs. However, as far as we know, there is no 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 expected time, assuming that partial 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.

CRediT authorship contribution statement

Tomoya Mori: Investigation, Writing – original draft. Tatsuya Akutsu: Conceptualization, Investigation, Writing – original draft.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.