An Efficient Algorithm for Computing Attractors of Synchronous And Asynchronous Boolean Networks

Biological networks, such as genetic regulatory networks, often contain positive and negative feedback loops that settle down to dynamically stable patterns. Identifying these patterns, the so-called attractors, can provide important insights for biologists to understand the molecular mechanisms underlying many coordinated cellular processes such as cellular division, differentiation, and homeostasis. Both synchronous and asynchronous Boolean networks have been used to simulate genetic regulatory networks and identify their attractors. The common methods of computing attractors are that start with a randomly selected initial state and finish with exhaustive search of the state space of a network. However, the time complexity of these methods grows exponentially with respect to the number and length of attractors. Here, we build two algorithms to achieve the computation of attractors in synchronous and asynchronous Boolean networks. For the synchronous scenario, combing with iterative methods and reduced order binary decision diagrams (ROBDD), we propose an improved algorithm to compute attractors. For another algorithm, the attractors of synchronous Boolean networks are utilized in asynchronous Boolean translation functions to derive attractors of asynchronous scenario. The proposed algorithms are implemented in a procedure called geneFAtt. Compared to existing tools such as genYsis, geneFAtt is significantly faster in computing attractors for empirical experimental systems. Availability The software package is available at https://sites.google.com/site/desheng619/download.


Introduction
Biological networks contribute a mathematical analysis of connections found in ecological, evolutionary, and physiological studies, such as genetic regulatory networks [1]. Pursuit for the nature of these networks is the central problem for biologists [2][3][4]. In the past decades, a wide variety of research has focused on modeling genetic regulatory networks using Boolean networks and search for their attractors [5][6][7][8][9][10]. Computing the attractors in the Boolean networks is critical in understanding corresponding genetic regulatory networks and coordinated cellular processes such as cell cycle and cell differentiation in living organisms [6,11]. In classical Boolean networks (CBNs), all nodes update their values at the same time called as synchronous Boolean network (SBNs). However, a criticism of CBNs as models of genetic regulatory networks is that genes do not update their values all simultaneously. To reflect this property of gene regulatory networks, Harvey and Bossomaier defined asynchronous Boolean networks (ABNs) where the random nodes were selected at each time and updated [12]. Since that, depending on the different update schemes, Boolean networks can be generally categorized into synchronous Boolean networks [7][8][9] and asynchronous Boolean networks [7,10,13,14]. For the same update schemes with different priority of activator or inhibitor in genetic regulatory networks, classical equations [15], prior inhibitor equations [10] and a combination of these two [16] are three types of Boolean translation functions. Therefore, given a Boolean network, there will be 2|3 different methods to represent its Boolean translation function.
In a synchronous Boolean network, all genes update their values simultaneously at consecutive time points. Heidel et al. [17] and Farrow et al. [18] have proposed a scalar equation approach to compute attractors in SBNs. Based on the former, Zhao [19] has proven that the way of computing attractors in SBNs is a NPcomplete problem. Dubrova et al. [20] have presented two tools -BooleNet and Bns -to compute attractors of SBNs. By contrast, in an asynchronous Boolean network, all genes update their values at different time points. Because each interaction between two nodes of a biological network follows distinct kinetics, it is generally thought that ABNs more realistically represent biological networks. However, due to the complexity of ABNs, the algorithms for computing network attractors are still mostly based on SBNs.
Previously, Garg et al. proposed a solution to compute the attractors in both SBNs and one class of ABNs [10]. First, they demonstrated that there were four types of attractors in a Boolean network: self loop, simple loop, syn-complex loop [or simple loop (type2)], and asyn-complex loop, shown as Fig. 1. The first two types (i.e. self loop and simple loop) were shared by SBNs and ABNs. But the latter two types, the syn-complex loop and the asyn-complex loop, were respectively found in SBNs and ABNs. Subsequently, they developed a series of algorithms which could be applied to compute the four types of attractors in a given Boolean network. Based on Garg's contribution, Ay F et al. [16] gave a faster method to list the states of self loops and one outgoing edge. Both Garg et al. and Ay et al. used the ROBDD data structure to support their algorithms.
Here, we developed two algorithms to further improve the computation of complex attractors in both SBNs and ABNs. First, based on the works of Garg et al., and Ay et al., we show that iterative computing can be used to accelerate the identification of the attractors of SBNs. Second, we develop a method to compute the asyn-complex loop (complex loop) using syn-complex loop, which allows us to simplify the computation of attractors of complex loops in ABNs. We have a software package to accomplish our two algorithms which are used to locate attractors of Boolean dynamic networks (for both SBNs and ABNs), with significantly reduced time. The structure of this paper is organized as follows: Section 2 gives the methods to compute attractors and splits them into two subsections. Section 2.1 presents iterative computing attractors' theory and its algorithms for SBNs. Section 2.2 proves a novel algorithm to locate attractors of ABNs from attractors of SBNs by asynchronous Boolean translation functions (ABTF). Section 3 tests feasibility and efficiency of our algorithm by several classical experimental benchmarks. Section 4 gives a conclusion and description of the future work.

Methods
This section gives two methods to compute attractors in both SBNs and ABNs. The first subsection presents iterative computing attractors' theory and its algorithms for SBNs. The second subsection provides a novel algorithm to locate attractors of ABNs from attractors of SBNs by asynchronous Boolean translation functions.

Computing Attractors in Synchronous Boolean Networks
In a synchronous Boolean network, all nodes update their values simultaneously at consecutive time points [18]. In another word, at a given time t[N, each node has only one Boolean value: 1 (Active) or 0 (Inhibit) [19]. Then, the equation of a synchronous Boolean network with n nodes is shown as Eq. 1 [19].
x n, tz1~fn (x 1, t , x 2, t , x 3, t , :::, x n, t ); ð1Þ where x i, t is a node in SBNs, f i represents the Boolean function of node x i, t , (x 1, t , x 2, t , x 3, t , :::, x n, t ) is a state in S, S denotes the universal set with 2 n different states, x i, t [f0,1g, 1ƒiƒn. It can be simplified as follows.
where X t~( x 1, t , x 2, t , x 3, t , :::, x n, t ), x i, t [f0,1g, 1ƒiƒn, F syn~( f 1 , f 2 , f 3 , :::, f n ) is the synchronous Boolean translation function (SBTF) from f0,1g n to f0,1g n . S i is a subset of S. FR(S i ,F syn ) is a set of forward reachable states, which are all the states that can be reached from the states set S i by F syn . BR(S i ,F syn ) is a set of backward reachable states, which are all the states that can reach the states set S i by F syn . All the states in FR(S i ,F syn ) or BR(S i ,F syn ) are different.   Remark 1. Definition 1 defines an attractor of a synchronous Boolean network. So similarly, we can define an attractor of an asynchronous Boolean network, when using F asyn instead of F syn , shown as Fig. 1(a),(b) and (d) in solid line boxes. F asyn represents an asynchronous Boolean translation function which will be introduced in section. We also use Atts syn and Atts asyn to represent all the attractors of a synchronous Boolean network and its asynchronous Boolean network, respectively. If a state X i is in an attractor, X j is one of its transient states, where X j [fBR(X i ,F ){FR(X i ,F)g, F [fF syn ,F asyn g, shown as Fig. 1 in dotted line boxes.
Because an attractor of a Boolean network is not known in advance, a common way to address this problem is setting a randomly chosen state as the initial state and exhaustively searching the entire state space. This approach has been successfully applied in several studies to compute the network attractors using empirically derived biological networks. However, the computational burden of this approach increases exponentially with respect to the number and length of attractors. Thus, it limits the application of this method for large biological networks.
Due to the recurrent nature of attractors, we reason that iterative computing algorithms can be applied on the Boolean translation functions of SBNs, like X tzj~Xt . An important implication is that identifying all attractors (Definition 1) does not require the computation of the entire state space. This suggests that we can use iterative computing to accelerate the identification of attractors in a given Boolean network. In the following, we present three theorems and their proof for iterative computation. Incorporating these theorems, an algorithm is demonstrated to compute attractors of SBNs.
Theorems of computing attractors using iterative computing in synchronous boolean networks. According to Eq. 2, it is easily inferred that X tz3~Fsyn (F syn (F syn (X t )))F 3 syn (X t ), where F 3 syn is synchronous Boolean translation function F syn after iteratively computing X t three times. Therefore, a simplified form of iterative computational equations is described as below.
where F j syn is F syn after j times iterative computing. F j syn can compute the state X t to state X tzj directly instead of j iterative computing steps by F syn . If state X tzj is same with state X t , that means state X t is in an attractor, which can be located as much as j steps iterative computation.
Definition 2. D j : it represents the states that will return to themselves after j iterations, where j[N z . This can be described by Eq. 4.
Definition 2 gives a simplified description of attractors, whose states could return to themselves after finite iterations. A shallow example can be supposed that, in a synchronous Boolean network, there are two attractors with length of 1 and 3 respectively. D 3 is the sum of the two attractors. Because the attractor whose length is 1 could also return to itself after 3 iterations. This feature can be proved by Theorem 1.
Theorem 1. For all S Att (Atts syn , if Length(S Att )~m, mDj (m is a factor of j), then, D j )S Att .
Proof. Let j~m|n, m,n[N z ; VX [S Att , Length(S Att )~m; According to Theorem 1, a set of attractors, whose length is j, can be located after j steps of iterative computing, shown as Eq. 5.
If a synchronous Boolean translation function F syn is same with F j syn , that means all the states can return to themselves by less than jz1 iterations. This is an important character to identify the numerous attractors in the SBNs, which has been proved by Theorem 2 and 3.
Theorem 2. Given a synchronous Boolean translation function F syn with n nodes, after d iterations, if D d~S , the period of this synchronous D d~S and Definition 2; The period of this synchronous Boolean network is d. Atts syn~Dd~S : An improved algorithm to compute attractors in synchronous boolean networks. Combining iterative computing (Theorem 1, 2, 3 and Eq. 5) and the ROBDD data structure, we have developed Algorithm 1 to compute attractors in SBNs. The input of Algorithm 1 is the synchronous Boolean translation function F syn ; its output is states of attractors (Atts syn ) and number of all attractors (Atts num). Specifically, Algorithm 1 starts with an initializing part, which initializes all the necessary variables. This is followed by a resolving part, which computes and deletes redundant attractors in a network. The resolving part further contains four components. The first component (Lines 12-13) will continue the next iterative computing and delete the visited attractors based on Theorem 1 and Eq. 5. The second component (Lines 15-20) judges whether synchronous Boolean translation functions are periodic or not, which has been proved by Theorem 2 and 3. The third component (Lines 22-27) verifies whether there will be a new attractor generated after one iteration. That is, if existing a new attractor, Algorithm 1 will add it into attractors (Atts syn ). Meanwhile, it will also update the attractors' number (Att num) and continue the next iteration. If not, Algorithm 1 will go to the next iterative computing. The last component (Lines 11, 29) contains the fix-point condition. When satisfying this condition, it will output attractors and number of attractors. For more detailed information, please read the Algorithms 1.
In the initializing part, j is the times of iterations. w is an empty set. Atts syn and Att num are the attractors and number of attractors, respectively. total set is the fix-point condition to judge whether the algorithm can be terminated or not. D j represents the states that will return to themselves after j iterations by synchronous Boolean translation function F syn . Atts syn , total set and D j are ROBDD data structures. In the main resolving part, D j :size() represents state number of D j . BR(D j ,F syn ) are all the states that can reach to D j by synchronous Boolean translation function F syn .
Algorithm 1 is different with Garg et al., which randomly picks up a state from state space and computes its forward reachable states to get an attractor. If you want to find out the attractors whose length is j, it needs to exhaustively search the state space. However, our algorithm can easily compute the same attractors in j times iterative computing.

Computing Attractors in Asynchronous Boolean Networks
As mentioned earlier, SBNs and ABNs differ in nodes updating schemes of Boolean translation functions. Instead of updating values of all the nodes simultaneously, ABNs only allow some of the nodes to update their values at a time point. For this reason, the computing of attractors in ABNs is more time consuming. Especially, it needs more intermediate steps when there are more than one bit different between two states.

Analysis of attractors in asynchronous boolean
networks. It is essential to give a simple description of types of attractors. As represented in Fig. 1, there are four types of attractors, self loop, simple loop, syn-complex loop and asyn-complex loop in both SBNs and ABNs. A self loop is a single state attractor, shown in Fig.1 (a). A simple loop includes two or more states, where every state is connected with only another state, and any two adjacent states differ from each other by only one bit, shown in Fig.1 (b). A syn-complex loop is similar to simple loop, but any two adjacent states differ from each other by more than one bit, shown in Fig.1 (c). A asyn-complex loop includes multiple interlinked states: every state is connected with more than one states, and there is only one bit different between any two adjacent states, shown in Fig.1 (d). Fig. 1 (a)(b)(c) and Fig. 1 (a)(b)(d) stand for the different types of attractors in SBNs and ABNs, respectively. According to the properties of self loops and simple loops, they can easily be identified in SBNs, which also are same in ABNs. Interestingly, a closer examination of the structure of the syncomplex loops and asyn-complex loops suggests that every asyn-complex loop contains one syn-complex loop or some transient states. This suggests that it is possible to use syn-complex loop to easily locate the states in asyn-complex_loop by asynchronous Boolean translation functions.

Algorithm 1: Iterative Computing Attractors on Synchronous Boolean Translation Function.
Function: Iterate Compute Attractors SBTF (F syn ) will compute all the attractors of SBTF iteratively; Input: The synchronous Boolean translation function F syn ; Output: All attractors (Atts syn ) and number of attractors (Att num); 1 begin 2//Initializing part 3 begin 4 int j/1; j is times of iteration 5 int Att num/0; //Att num is the number of attractors 6 set Atts syn /w; Atts syn is a set of attractors 7 set total set/S; total set is a set of unvisited states 8 set D 0 /w; Initialize D 0 as empty set 9 end 10 //Main resolving part 11 while total set=w do 12 report attractors D j and its number D j :size()=j 24 Att num/Att numzD j :size()= j; //Update the value of Att num 25 Atts syn /Atts syn zD j ; //Add the new attractors into Atts syn 26 total set/total set{BR(D j ,F syn ); //total set deletes the states can reach D j 27 end. 28 end. 29 return Atts syn , Att num; 30 end.
One example is shown in Fig. 2, here i=m, m=n, i=n, i,m,n[N z . Fig. 2(a) is an asynchronous attractor, where the current state and its next states are different by one bit. Suppose that the i th bit of the state s 1 and s 2 is different. If i~m, it means that state s 1 and s 3 are the same. The situation is also true for i and n. If state s 3 and s 4 are different at the n th bit, then state s 2 and s 5 must differ at the n th bit. Otherwise, state s 4 cannot reach state s 5 by changing one bit. When state s 5 and s 6 differ at the n th bit, state s 6 and s 1 will be different at the i th bit, and vice versa. Fig. 2(b) shows the corresponding synchronous attractor to Fig. 2(a). The difference is that state s 2 and s 4 differ in the m th and bits simultaneously. Other relations are the same except for state s 3 . Therefore, an asyn-complex_loop contains one syn-complex loop or some transient states. That means we can use syn-complex loop to easily locate the states in asyn-complex_loop by asynchronous Boolean translation function F asyn .
An algorithm to compute attractors in asynchronous boolean networks. An important implication of the above analysis is that the attractors of the ABNs can be derived from the attractors of the SBNs using synchronous and asynchronous Boolean translation functions. Therefore, we use attractors computed by Algorithm 1 in SBNs as the basis input for the new algorithm (Algorithm 2) to compute attractors of ABNs.
Specifically, Algorithm 2 can be divided into three parts: initializing part (Lines 3-6), main resolving part (Lines 8-41) and checking unvisited states part (Lines 44-52). The initializing part also initializes all necessary variables in Algorithm 2. The main resolving part also can be split into five components. The first component (Lines 11-15) means sk is a self loop. The second component (Lines 18-21) means FR(s,F syn ) is a simple loop. The third component (Lines 24-29) means FR(s,F asyn ) is an unvisited asyn-complex loop and FR(s,F syn ) is an unvisited syn-complex loop. The forth component (Lines 32-34) means FR(s,F asyn ) is a visited asyncomplex loop and FR(s,F syn ) is an unvisited syn-complex loop. The fifth component (Lines 37-40) means FR(s,F asyn ) is the set of transient states and FR(s,F syn ) is an unvisited syn-complex loop. After the main resolving part, if there exists unvisited states, Algorithm 2 will go to the checking unvisited states part. This part will check the unvisited states and report the left asyn-complex loops. For more detailed information, please read the Algorithms 2.
Furthermore, in the initializing part, s is a state, w is an empty set, asyn state set is a set of recording unvisited states by F asyn , S is the universal set. In the main resolving part and checking unvisited part, pick up one state(Atts syn ) represents picking up anyone state from attractors Atts syn of SBNs. FR(s,F syn ) and FR(s,F asyn ) are the reachable states from s by synchronous Boolean translation function F syn and asynchronous Boolean translation function F asyn , respectively. BR(s,F asyn ) is the reachable states to s by asynchronous Boolean translation function Algorithm 2: Compute and classify Four Types Attractors.
Function. Compute Attractors() computes and classfies attractors as four types, shown in Fig. 1; Input. Attractors of SBN(Atts syn ), SBTF(F syn ) and ABTF(F asyn ); Output. Four types attractors (a)(b)(c)(d) of SBTF and ABTF shown as Fig. 1; 1 begin 2 // Initializing part 3 begin 4 set s/w; // s is a state 5 set asyn state set/S; / asyn state set is a set of unvisited states by F asyn 6 end 7 // Main resolving part 8 while (Atts syn =w) do 9 s/pick up one state(Atts syn ); // s is any one state in Atts syn 10 // s is a self loop, shown as Fig. 1

Results and Discussion
We have implemented our methodology in a software package called geneFAtt, which is based on a ROBDD data structure named BuDDy [21]. In this package, there are two parts, source code and benchmark. The source code part implements functions of Algorithm 1 (iterative computing attractors in synchronous Boolean networks) and Algorithm 2 (computing and classifying attractors as four types of synchronous Boolean networks and asynchronous Boolean networks). Algorithm 1 computes the attractors of SBNs, which are then used as the input of Algorithm 2 to classify the attractors of ABNs. Because the self loop ( Fig. 1(a)) and simple loop ( Fig. 1(b)) of SBNs are same with its corresponding type of attractors of ABNs, respectively. We can easily use syncomplex loops (Fig. 1(c)) to locate the asyn-complex loops ( Fig. 1(d)) by the asynchronous Boolean translation functions. The benchmark part includes five biological networks, Mammalian Cell [7], T-helper [22], Dendritic Cell [10], T-cell Receptor [23], and Protein-ex [17]. We ran geneFAtt and genYsis [10] with the five synchronous/asynchronous biological models and have compared their running results.
As shown in Table 1, the first four biological networks, Mammalian Cell [7], T-helper [22], Dendritic Cell [10] and T-cell Receptor [23] have been studied in [10]. Protein-ex is an extended case from Heidel et al. [17]. It also represents a kind of biological networks that all the states are in attractors. Because we adopt to the same methods of modeling the SBNs and ABNs to setup the synchronous and asynchronous Boolean translation functions refereed to Garg et al. Meanwhile, using the same inputs, geneFAtt can obtain the same attractors as those by genYsis [10] shown in Table 1, which suggests that our algorithms are valid for the analysis of biological networks.
With the same validity, geneFAtt shows more efficient and feasible compared to genYsis. Table 2 gives the running time and running time efficiency ratio (RTER, Eq. 6) of the five biological networks which have been produced by the two procedures genYsis and geneFAtt. T 1 and T 2 represent the running time on every biological network of the two procedures genYsis and geneFAtt, respectively. All experiments are performed on an IntelH Core TM CPU 4300 1.80 GHz with 2GB memory and a Ubuntu 9.04 Linux server. Importantly, we note that the running time of geneFAtt is much more shorter than genYsis [10]. Specifically, compared with genYsis, geneFAtt improves the running time of Mammalian Cell [7], T-helper [22], Dendritic Cell [10], T-cell Receptor [23], and Protein-ex [17] Figure 2. Diagrams of an attractor in asynchronous (a) and synchronous (b) Boolean networks. Each state is represented by a circle, and is designated as S n . The variable i represents that the i th bit of the state s 1 and s 2 is different, which is also same as m and n. The numbers ½m,n indicate that state s 2 and s 4 differ by the m th and n th bits respectively. The (i,n) and (n,i) represents when state s 5 and s 6 o differ at the n th bit, state s 6 and s 1 will be different at the i th bit, and vice versa. The difference between the two representations (i.e. synchronous versus asynchronous) of the attractor is that s 2 and s 4 differ in the m th and n th bits, m=n, m,n[N z . That means we can use syn-complex loop to easily locate the states in asyn-complex_loop by asynchronous Boolean translation function F asyn . doi:10.1371/journal.pone.0060593.g002

Conclusions
This paper has addressed a method to compute attractors in SBNs and ABNs. We have developed a new iterative computing algorithm to identify the attractors Atts syn of SBNs. Meanwhile, another computing and classifying algorithm has been proposed to locate the attractors of ABNs rapidly. Our approaches give a significant acceleration of computing and locating attractors in biological networks. It is a big challenge that how to identify attractors in the large biological networks. Based on above research, it is expected that we can find a pathway to resolve this problem from the modeling methods of biological networks.