### Algorithm

We begin with simple examples to explain the basic idea of our proposed algorithm using *a priori* information.

First, we consider the simplest case of *n* = 1 and the target attractor is a point attractor such that *x*_{1} = *b*_{1} (we omit time step *t* because we consider a point attractor), where *b*_{1} corresponds to a point attractor and takes 0 (resp., 1) with probability 0.5. Suppose that we examine *x*_{1} = 1 first, and then *x*_{1} = 0. Since the first trial succeeds with probability 0.5, the number of expected trials is . On the other hand, suppose that we know that *b*_{1} = 0 holds with probability . In this case, we examine *x*_{1} = 0 first, and then *x*_{1} = 1. Then, the expected number of trials is , which is less than .

Next, we extend this example for the case of *n* = 2. Suppose that the target attractor is a point attractor such that *x*_{i} = *b*_{i} (*i* = 1, 2), where each *b*_{i} takes 0 (resp., 1) with probability 0.5. If there is no prior information, we examine all 2-bit patterns in any order, for example, 11,10,01,00. Then, the expected number of trials is . On the other hand, suppose that we know that *b*_{i} = 0 holds with probability for each *i* ∈ {1, 2}. In this case, we examine 2-bit patterns in the following order: 00,01,10,11. Then, the expected number of trials is , which is smaller than .

We can extend this idea to general *n*. The algorithm is quite simple (although its analysis is involved). It starts testing for attractors from the most plausible vector **x**_{0} ∈ {0,1}^{n} with following the trajectory starting from it. If the search from this vector fails, the algorithm examines trajectories from the vectors each of which has one bit different from the first one. It further repeats the same procedure by increasing the number of bits different from the first vector. The following gives the formal description of the algorithm, where we assume without loss of generality (w.l.o.g.) that the prior probability of 0 is larger than that of 1 for all bits (otherwise, it is enough to exchange the roles of 0 and 1 for the corresponding bits).

- For
*k*= 0 to*n*, do STEP 2 and 3. - For all
**x**_{0}∈ {0,1}^{n}that contains*k*bits with value 1, do STEP 3. - Let
**x**(0) =**x**_{0}. Compute**x**(*t*+ 1) =**f**(**x**(*t*)) repeatedly until**x**(*t*+ 1) =**x**(*t*′) holds for some*t*′ ≤*t*(which means 〈**x**(*t*′),…,**x**(*t*)〉 is an attractor). If the attractor is a desired one, output it and exit.

In the following, we refer to this algorithm as **ATTapriori**. Note that how to decide whether the current attractor is a desired one is not a trivial task. If we know some exact criteria (e.g., allowable range of the attractor period and/or states of specific nodes), we can add a subroutine to check it. Otherwise, the algorithm can be terminated if the number of trials exceeds a specified number or CPU time exceeds the time limit. Then, we may manually check (maybe with some user-customized computer program to rule out obviously non-desired ones) whether or not there exists a desired one among the attractors reported so far.

It is to be noted that the algorithm can be modified so that it can enumerate all attractors by removing “If the attractor is a desired one, output it and exit.” of STEP 3 and merging the identical attractors. However, it would take more than *O*(2^{n}) time (because it will examine all 2^{n} starting vectors) and thus the resulting algorithm is meaningless. Since the purpose of this algorithm is to make use of *a priori* information on some global state in a specific attractor, it is reasonable that the algorithm is not useful for enumerating all attractors.

### Time complexity analysis

Here we analyze the expected time complexity of **ATTapriori**.

We consider a BN with *n* nodes. Suppose that *x*_{i} = 0 holds with probability *p* for all *i* = 1, …, *n* in some specific global state **x**_{t} of the target attractor, where *p* > 0.5. Note that if *x*_{i} = 1 has the probability *p* (*p* > 0.5) it is enough to exchange the roles of 0 and 1 for such nodes. In theoretical analyses, the objective is to give an upper bound of the number of trials until the algorithm examines **x**_{t} as a starting vector. Therefore, we need to modify STEP 3 as follows:

If **x**_{0} = **x**_{t}, output it and exit.

Of course, this modified algorithm is meaningless in practice because it is impossible to know **x**_{t} in advance. However, the expected number of trials (i.e., the expected number of tested **x**_{0}s) for this modified algorithm gives an upper bound of that for the original algorithm because it examines the whole trajectory beginning from **x**_{0} (much more than one global state **x**_{0} per trial).

**ATTapriori** examines bit vectors from those with a smaller number of 1s to a larger number of 1s (e,g., 000, 001, 010, 100, 011, 101, 110, 111) where ties can be broken in an arbitrary manner. We can see that the expected number of trials *E*(*n*, *p*) until reaching a desired attractor is given by

where

Note that *p*^{n−k} ⋅ (1 − *p*)^{k} is the success probability per vector **x**_{0} containing *k* bits with value 1, and *S*(*n*, *k*) + *j* denotes the number of trials (i.e., the number of different staring vectors) until the current vector (i.e., *S*(*n*, *k*) + *j*th vector) is examined.

**Lemma 1**. *The expected number of trials E*(*n*, *p*) *is O*(*f*(*α**)^{n} ⋅ *n*^{2}), *where* , *f*(*α*) = *p*^{1−α}(1 − *p*)^{α} *β*^{2}, *and* .

*Proof*. In order to analyze the order of *E*(*n*, *p*), we divide this expectation due to vectors **x**_{0} containing at most ⌈*n*/2⌉ bits with value 1 (i.e., at most half of the bits have value 1) and those containing at least ⌈*n*/2⌉ + 1 bits with value 1. First, we evaluate the former part, that is, the partial sum *F*(*n*, *p*) of *E*(*n*, *p*) defined by

Since holds for *k* ≤ ⌈*n*/2⌉, we have

Next, we need to find an upper bound of ‘max’ in the last part of the above inequality. To this end, we let *k* = *αn*. By using upper and lower bounds of Stirling’s approximation , we have

where . Here we define *β* by

Then, the exponential factor of *F*(*n*, *p*) is bounded by max_{α≤0.5} *f*(*α*)^{n} where

Since it is difficult to directly maximize *f*(*α*), we derive *α* maximizing *g*(*α*) where

By differentiating *g*(*α*) with respect to *α*, we have

By solving the last equality, we have

Furthermore, we can verify that holds. In order to evaluate the partial sum of *E*(*n*, *p*) due to vectors containing at least ⌈*n*/2⌉ + 1 bits with value 1,

holds for some constant *γ* > 0. Since *S*(*n*, *k*) ≤ 2^{n} holds for any *k* ≤ *n*,

holds for some constant *d*. Thus, by letting , the partial sum of *E*(*n*, *p*) due to the latter half of vectors is bounded as

Therefore, the expected number of trials *E*(*n*, *p*) is bounded as follows:

In Lemma 1, we assumed that *Prob*(*x*_{i} = 0) = *p* holds for all *x*_{i}. However, the probability is not usually the same. Therefore, we consider the case where *Prob*(*x*_{i} = 0) = *p* + *a*_{i} holds for all *x*_{i}, where *a*_{i} ≥ 0 (*i* = 1, …, *n*) and . That is, we let , where we can exchange the roles of *x*_{i} = 0 and *x*_{i} = 1 if . Note that we are considering a situation that we have *a priori* information of some specific state in the target attractor and thus this assumption is reasonable.

**Theorem 1**. *Suppose that x*_{i} = *b*_{i} (*b*_{i} ∈ {0, 1}) *holds with probability greater than or equal to p* (*p* > 0.5) *in one global state in the target attractor. Then*, **ATTapriori** *finds the target attractor using an O*([*p*^{1−α}(1 − *p*)^{α} *β*^{2}]^{n} ⋅ *n*^{2}) *expected number of trials, where* , .

*Proof*. The proof strategy is to show that the expected number of trials in this setting is upper bounded by that in the case of *Prob*(*x*_{i} = 0) = *p* for all *x*_{i}.

Let *p*_{i} = *Prob*(*x*_{i} = 0) = *p* + *a*_{i} (*a*_{i} ≥ 0) and *D*_{P} denote the corresponding probability distribution. Let **b** denote the step when **b** is examined, where **b** is a 0–1 vector of length *n*. Let **b**_{i} denote the *i*th element (*i.e*., *i*th bit) of **b**. Let *δ*(*x*, *y*) = 0 be the delta function: *δ*(*x*, *y*) = 1 if *x* = *y*, otherwise *δ*(*x*, *y*) = 0. Recall that **ATTapriori** examines bit vectors from those with a smaller number of 1s to a larger number of 1s, where ties can be broken in an arbitrary manner. For example, in the case of *n* = 3, we can examine in the order of 000,001,010,100,011,101,110,111, where we have *I*_{000} = 1, *I*_{001} = 2, *I*_{010} = 3, ⋯.

Then, the expected number of trials *E*(*n*, *D*_{P}) is given by

because the success probability for vector **b** is given by

and the number of trials done until the examination of **b** is **b**. We prove the theorem (*i.e*., *E*(*n*, *D*_{p}) ≤ *E*(*n*, *p*)) by mathematical induction on the number of bits with *a*_{i} > 0.

In the base case, *p*_{i} = *p* (i.e., *a*_{i} = 0) holds for all *i* = 1, …, *n*. Therefore, the theorem follows from Lemma 1.

In the inductive step, we show that the expected number of trials does not decrease if we change one *p*_{i} = *p* + *a*_{i} with *a*_{i} > 0 to *p*_{i} = *p*. Suppose that the theorem holds for any in which *h* nodes satisfy *p*_{i} > *p*. Let *D*_{P} be a distribution in which *h* + 1 nodes satisfy *p*_{i} > *p*. We assume w.l.o.g. that *p*_{1} = *p* + *a*_{1} > *p*. Let be the distribution that is obtained from *D*_{P} by changing the value of *p*_{1} to *p*. From the induction hypothesis, we have

Note that *p*_{1} = *p* is used in calculation of *p*_{b}. Let *B*_{0} (resp., *B*_{1}) be a set of 0–1 vectors **b** such that **b**_{1} = 0 (resp., **b**_{1} = 1).

Then, *E*(*n*, *D*_{P}) is written as

because prior probabilities *p* and 1 − *p* for the first bit (i.e., *b*_{1}) in are replaced by *p* + *a*_{1} and 1 − (*p* + *a*_{1}) in *D*_{P}. Let **b**^{(0)} (resp., **b**^{(1)}) be a bit vector obtained from **b** by letting **b**_{1} = 0 (resp., **b**_{1} = 1). Define *C*_{0} and *C*_{1} by

Here we note that holds because the number of bits with 1 of **b**^{(0)} is smaller than that of **b**^{(1)}. Since holds from the definition of *p*_{b}, we have . Therefore, we have

which completes the mathematical induction.

It is to be noted that the expected time complexity of **ATTapriori** is *O*([*p*^{1−α}(1 − *p*)^{α} *β*^{2}]^{n}*poly*(*n*)) if each Boolean function can be evaluated in polynomial time and the length of each trajectory (including an attractor cycle) is polynomially bounded.

Since most practical Boolean functions can be evaluated in polynomial time and the lengths of trajectories in most practical BNs are considered to be not very long, this is a reasonable assumption.

Even if the trajectory is not bounded by a polynomial, we can modify the algorithm so that each trial is finished if the length of the trajectory exceeds some given steps (some polynomial steps). Since the proofs of Lemma 1 and Theorem 1 only discuss whether a given initial state is the same as the target state, this modification does not affect these theoretical results. It may affect the correctness of the original algorithm (i.e, the original algorithm may miss the desired attractor) because the period of the desired attractor may not be bounded by a specified polynomial, where the period of an attractor corresponds to a cyclic part of the trajectory (not the whole part of the trajectory). However, it is reasonable not to consider very long (e.g., exponentially long) attractors because any usual organism cannot live for exponentially long periods. For example, suppose that 1 step corresponds to 1 micro second (i.e., 10^{−6} second). Then, even for *n* = 100, 2^{n} is greater than 4 × 10^{16} years. Limiting the length of the trajectory is also useful to reduce the space complexity. The algorithm need to memorize all states in the trajectory in order to find a cycle, which implies that the worst-case space complexity would be *O*(2^{n} ⋅ *n*) because the length of the trajectory can be *O*(2^{n}) in the worst case. However, if we limit the length of the trajectory by some constant and it is enough to find one desired attractor, the space complexity would be *O*(*n*).

In some cases, satisfy that *Prob*(*x*_{i} = 0) ≥ *p* but there is no information on . In this case, for *k* = 0, …, *n*_{1}, we examine assignments on , each of which has *k* 1s. Furthermore, we examine assignments on for each of assignments in the random order, where *n*_{2} = *n* − *n*_{1}. Then, the resulting expected number of trials is

Since is from Lemma 1, this number is

Suppose that *n*_{1} = *cn* holds for some constant *c*. Then, this number will be

Since *f*(*α**) < 2, this number is *O*((2 − *δ*)^{n}*n*^{2}) for some constant *δ*, where *δ* depends on *p* and *c*.

The algorithm can be extended for the case in which there exist several levels of uncertainty (*i.e*., we can use several *p*_{1}, *p*_{2}, …, *p*_{h} depending on node types). Note that the algorithm is deterministic and terminates once the desired attractor is found. Of course, how to define the desired attractor needs some discussions as discussed before. If we know some criteria to automatically decide whether or not a given attractor is a desired one, we can terminate the algorithm as soon as such an attractor is detected. Otherwise, the algorithm can be terminated if the number of trials exceeds a specified number or CPU time exceeds the time limit, where the number may be determined based on computational experiments as done below.

The **ATTapriori** algorithm is implemented in C, and the source code is available at https://github.com/takutsu5/AttPrior. In this version, at most 3 different values can be specified as probabilities (including probability 1.0) and bit vectors are examined from higher probabilities one to lower ones under a constraint that the number of changed bits for each probability class must not exceed the specified threshold.

### Analysis of synthetic random *N*–*K* networks

In order to evaluate the efficiency of the proposed algorithm, we performed computational experiments to detect attractors using synthetic random *N*–*K* networks, where *N* and *K* indicate the number of nodes and the indegree, respectively, and the update function of each node is controlled by exactly *K* nodes which are uniformly randomly selected. In this section, random *N*–*K* networks were generated using the *generateRandomNKnetwork* function provided by the R-package *BoolNet* [27].

First, we measured the average number of trials of 1000 experiments until a bit vector of global state examined by **ATTapriori** algorithm corresponds with the target attractor. In this experiment, random *N*–*K* networks of *N* = 5, 10, 15, 20, 25, 30, 35, 40, and *K* = 2 were generated, and one of the attractors detected by the *getAttractor* function of *BoolNet* was set as the target attractor for each network. In order to avoid matching with the target attractors in the first trials, the bit vectors generated by changing each bit of the target attractors according to the prior probability are given as the initial states.

We assigned the same prior probability *p* to all nodes, and changed *p* from 0.7 to 0.9 in steps of 0.05, and then the upper bound of the expected number of trials given by Theorem 1 was calculated for each *p*.

The results of comparing the average number of trials with the theoretical value using random *N*–*K* networks are shown in Fig 2A. The average number of trials increases exponentially with the increase of the network size *N*, and is upper-bounded by the expected number of trials obtained theoretically. On the other hand, the average number of trials decreases as the prior probability increases, and 3.50 × 10^{10} and 9.62 × 10^{6} trials are required when the prior probabilities are 0.7 and 0.9 in the case of *N* = 40, which indicates that the higher prior probability enables us to detect attractors much faster than without *a priori* information.

Fig 2. Analysis of synthetic random networks and the angiogenesis network.

(A) The number of expected trials in synthetic random networks. The average number of trials and the expected number of trials for attractor detection of synthetic random *N*–*K* networks with the prior probability from 0.7 to 0.9. The horizontal axis shows the size of input random *N*–*K* networks, where *K* is fixed at 2. The dotted lines indicate the regression lines for the experimental results. (B) The lengths of trajectories of synthetic random networks (RBNs) and the angiogenesis network. The average lengths of trajectories of 100 trials from random initial states to attractors for random *N*–*K* networks (*N* = 5, 10, 20, 50, 100 and *K* = 2) and the angiogenesis network (*N* = 142). The symbol ‘*’ and the bars indicate the averages of period lengths and the standard errors, respectively. (C) The total lengths of trajectories required to reach the target attractors. The **ATTapriori** and **Random** are the cases with and without *a priori* information, respectively. The result shows that the average total lengths of trajectories of 100 experiments not including period lengths until reaching the target attractors for the BN sizes of *N* = 10, 20, 30, and 40, where the prior probability in **ATTapriori** was set to 0.90.

Second, we performed another experiment to show the distribution of the lengths of trajectories until reaching attractors. Fig 2B indicates the average length of trajectories of 100 experiments from random initial states to attractors without period lengths for random *N*–*K* networks whose sizes are *N* = 5, 10, 20, 50, 100 and indegrees are *K* = 2. As a real network, the lengths of trajectories of the angiogenesis network of size *N* = 142 was also measured (its details are described in the next section) [9]. The symbol ‘*’ indicates the averages of period lengths. Although the average lengths of trajectories becomes longer as the size gets larger, it is less than 300 even when the size of network is 100. In addition, the average length of trajectories of the angiogenesis network is much shorter than that of random networks of size *N* = 100.

Finally, we compared the total lengths of trajectories until reaching a target attractor with (**ATTapriori**) and without (**Random**) *a priori* information. In this experiment, a target attractor was set in advance, which was detected by *BoolNet* package of R. The **Random** method starts from a random initial state and repeats state transitions until an attractor is found. If the found attractor is the target attractor, it outputs the total lengths of trajectories and terminates. Otherwise, it starts from another random initial state and repeats the above procedure. On the other hand, in **ATTapriori**, the same procedure was done except that the initial states are decided by **ATTapriori** algorithm. Fig 2C shows that the average total lengths of trajectories of 100 experiments not including period lengths until reaching the target attractors for the sizes of *N* = 10, 20, 30, and 40, where the prior probability in **ATTapriori** was set to 0.90. From the results, the average total length of trajectories of **ATTapriori** was shortened to about 1/40 at the maximum. Additional information and files can be found in S2 File.

### Effect of microenvironment on EC behavior revisited (BN with *n* = 142 nodes)

We next applied the proposed algorithm to the analysis of a BN describing the effect of several microenvironments on EC behavior during angiogenesis [9]. The authors of the original study developed two BN versions: a full BN with *n* = 142 nodes which we analyze here, and a reduced BN with *n* = 64 nodes, which was analyzed in the original study. The two BNs are formulated in the classical way, where each node corresponds to a molecular entity or a conceptual placeholder (such as *shear stress*). Using the reduced model, the authors identified in total 67 microenvironments of which 35 induce typical, and 32 atypical EC behavior. The authors used three signatures of four molecular markers (see Introduction; Materials and methods) to decide which EC behavior (tip, stalk, or phalanx) was triggered by each microenvironment. A microenvironment induces typical behavior if all detected attractors show an EC marker configuration corresponding to a single and stable EC signature. Otherwise, the microenvironment is regarded as inducing atypical behavior. For example, if the EC marker signature does not correspond to any of the three behaviors, the induced behavior is regarded as atypical. Here, we investigated the attractor landscape in the full model with *n* = 142 nodes as provided by the authors [9]. We use the same definition as the authors to interpret EC behavior based on the signature of the four EC markers, and use the same microenvironments inducing a specific behavior. (*I.e*., microenvironment numbers 1–2, 3–16, and 17–35 induce phalanx, stalk, and tip, respectively; numbers 36–67 induce atypical behavior). Furthermore, our analysis requires initial guesses for each microenvironment, and *a priori* probabilities. We base our initial guesses on the in total 67 microenvironment configurations by the authors of the original study [9], which each comprise 16 nodes. The authors also include 10 nodes with constitutive activity, which we also incorporate in our initial guess. Hence, in total, we can assign 26 node values to our initial guess based on the original study, and assume that their *a priori* probability is 1 (Materials and methods). Since the initial values of all of the remaining 116 nodes are not reported in the original study [9], we test two scenarios: one in which they are set to 0 (OFF scenario), and one in which they are set to 1 (ON scenario). We call this type of guess with partially known, and partially unknown initial values a semi-informed guess. For these remaining nodes, we test for each of these scenarios (ON and OFF) four arbitrarily chosen *a priori* probabilities: 0.7, 0.8, 0.9, and 1 (Materials and methods).

First, our attractor analysis shows that for the microenvironments inducing typical behavior (numbers 1 to 35), point attractors or periodic attractors with period 2 can be identified (Fig 3 and Fig A in S1 File). The detected periodic attractors for the microenvironments (numbers 36 to 67) inducing atypical behavior show higher variation in their lengths, with *per* = 2, 4, 6, 8, 18 (Fig A in S1 File).

Fig 3. Typical endothelial cell (EC) behaviors induced by 35 microenvironments.

Comparison of the induced EC behavior from 35 microenvironments (rows) to the results from the original study [9], *Original* column, with our analysis, *ON* and *OFF* columns. In the underlying BN three types of EC behavior, phalanx, stalk, and tip, were predicted to be induced dependent on the microenvironment. In the original study [9] based on a BN with *n* = 64 nodes (*Original*), microenvironments 1–2 induce phalanx behavior (blue), 3–16 induce stalk behavior (yellow), and 17–35 induce tip behavior (grey). We used the full network with *n* = 142 nodes for our analysis requiring an initial guess and *a priori* probabilities. We could assign 26 node values based on the original study, and tested two scenarios for the remaining 116 nodes as our initial guess: the ON scenario, in which the remaining nodes are set to 1 (*ON* column), and the OFF scenario (*OFF* column), in which the remaining nodes are set to 0. Results using *a priori* probability 0.7 shown. For both types of guesses (ON and OFF scenario), both point, and periodic attractors with period 2 were detected. Per.: periodic. *No periodic attractors detected. (The same microenvironment numbering is applied as in the original study).

Second, we confirm that the 33 microenvironments (numbers 3 to 35) previously found to induce either stalk or tip behavior induce the expected behavior in the full network (Fig 3). Similarly, all 32 microenvironments (numbers 36 to 67) identified to induce atypical behavior also induce atypical behavior in the full network (Fig B in S1 File).

Third, the two microenvironments (numbers 1, 2) previously described to induce only phalanx behavior were found to induce two types of behavior: phalanx and stalk behavior (Fig 3). This finding indicates that the conclusion about the two previously defined microenvironments (numbers 1,2) in the reduced network does not directly apply to the full network.

#### Analysis of a cell cycle control network (bBM with *n* = 3158 nodes).

We next considered a cell cycle control network of *S. cerevisiae* [28] formulated as a bBM with 3158 nodes. Due to the bBM formalism, a few considerations need to be taken for the attractor analysis (see Introduction). The **ATTapriori** algorithm may find attractors in which two or more mutually exclusive states (*e.g*., DNA both replicating and replicated at the same time) are simultaneously true, or a necessary component is not true. While these attractors are technically possible, they do not convey any biological meaning. We regard these attractors as technical artefacts and refer to them as *invalid*. In addition, the algorithm may remove a gene or the last true state of a specific site, effectively creating a deletion mutant of the corresponding component(s) and hence altering the genotype of the model. Next, we regard all other attractors to correspond to a biologically relevant genotype, with either of two possible corresponding phenotypes: viable or lethal. We regard an attractor to correspond to a viable phenotype if all macroscopic processes are traversed, as this would be required for a cell to successfully duplicate. This is only possible for periodic attractors. Similarly, we regard an attractor to correspond to a lethal phenotype if the attractor is a point attractor, or a periodic attractor which does not traverse all macroscopic states. Such a lethal phenotype may correspond to a mutant in which an essential gene is missing.

First, we studied the attractors of the genotype corresponding to the wildtype phenotype where we used the previously known periodic attractor as initial guess. We regard this an informed guess, where we also set two necessary network inputs *Nutrients, Histones* to 1, and four other inputs, corresponding to chemical inhibitors of the cell cycle, to 0. We tested four different, arbitrarily chosen *a priori* probabilities, 0.7, 0.8, 0.9, and 1 in different combinations, except for the in total six inputs, where the *a priori* probabilities are all set to 1 (see Materials and methods). We discovered 66 unique periodic attractors with period 186, and 125 unique point attractors (Fig 4). We then analyzed these 66 unique periodic attractors (Table A in S1 File) by manually scanning for nodes that are constitutively 0, and compare this to the wildtype attractor. If the constitutively inactive states all belong to one gene, we consider this gene to be absent in the attractor. We found that one attractor corresponds to the wildtype, 39 correspond to viable deletion mutants (gene-specific sites at The Saccharomyces Genome Database (SGD) [31]), two correspond to deletion mutants with conflicting reported phenotypes (PP2A, Net1), and 16 to deletion of hypothetical kinases and phosphatases introduced in the model construction process (not testable). Only one corresponds to a lethal deletion mutant; *sec2*, where the essential function is not part of the model. The remaining 7 are technical artefacts, where an infeasible pattern of macroscopic states was active in the initial perturbation. Hence, all the newly discovered periodic attractors correspond either to deletion mutants or to technical artefacts, and no new basins of attractions in the wildtype were discovered. The 125 unique point attractors (Table B in S1 File) can by definition only correspond to lethal phenotypes. We manually examined the attractors for inactive nodes, and assigned lethal mutations if the inactive nodes corresponded to one or several genes. We found that 71 of the 125 detected point attractors indeed correspond to known lethal mutants. 28 mutants have a viable phenotype *in vivo*, but redundancy mechanisms that compensate for loss of their functions are not included in the model. 18 mutants showed a missing structural component, and hence, the phenotype is lethal (model phenomenon, technical artefact). The remaining eight attractors constitute technical and other mutants which could not be clearly categorized, or where the reported *in vivo* phenotype is conflicting.

Fig 4. Attractor analysis results of cell cycle control network.

Attractor analysis result of a bBM cell cycle control network with *n* = 3158 nodes. Results from in total eight initial guesses are shown: the wildtype genotype, the *clb1* genotype with a corresponding viable phenotype, the *cln3* genotype with a corresponding lethal phenotype (due to absence of Bcks2 in the underlying model), the *clb1clb2* genotype with a lethal corresponding phenotype, and four initial guesses where one of the following cell cycle inhibitors is active, hydroxyurea, latrunculin, nocodazole, and pheromone. Several combiantions of *a priori* probabilities were used, representative results shown. (A) Barplot of the found types of attractors (point or periodic) for the eight initial guesses. (B) Data table. WT: wildtype, HU: hydroxyurea, Lata: Latrunculin A, Noco: nocodazole, Pher: Pheromone, per: period.

Fig 4 also shows the attractors obtained using other initial guesses: the viable *clb1* mutant phenotype, and the two lethal mutant phenotypes *cln3* (due to the absence of Bck2) and *clb1clb2*.

Second, we perturbed the network with the four cell cycle inhibitors hydroxyurea (HU), Latrunculin A (LatA), nocodazole (Noco), and pheromone (Pher), and analyzed the attractors. Perturbing the cell cycle network with either of the four inhibitors arrests cell cycle progress and hence, we expect point attractors. Fig 4 shows that the inhibitory effect is reflected in the identified attractors, as only point attractors were observed, corresponding to lethal phenotypes.

Third, we introduced two types of functional mutations: AND to OR mutations (functional promiscuity), and OR to AND mutations (functional restriction), with mutation rates of 0.1%, 1%, and 10%, and calculated the attractors. The two types of functional mutations promote different kinds of attractors with increasing mutation rates: Functional promiscuity mutations promote periodic attractors (Fig C in S1 File), whereas functional restriction mutations promote point attractors (Fig D in S1 File). Overall, the periodic attractors had a high variation, with lengths of 2, 3, 4, and 32, and none of them was found to correspond to a viable phenotype. Furthermore, with increasing mutation rate, the detected attractors did no longer correspond to a biologically meaningful pattern and are regarded as invalid. In comparison, in the wildtype and *clb1* mutant (Fig 4), periodic attractors corresponding to viable phenotypes were found. The majority of the point attractors detected in the *cln3* and *clb1clb2* mutants have a biologically valid signature, and correspond to lethal phenotypes. We observe that mutations in the Boolean functions with rates as low as 0.1% induce lethal behavior.

Finally, we tested random initial states as initial guesses. Our previously known periodic attractor corresponding to the wildtype is an informed guess about an existing periodic and biologically relevant attractor. Here, we used 10 random initial states to test if our algorithm can detect periodic and/or biologically relevant attractors with random, and hence, uninformed guesses. We tested three scenarios, one with completely random initial guesses, one in which the input *Nutrients* was set to on, and one in which *Nutrients* was set to off. In all three scenarios, only point attractors were detected (Fig E in S1 File).