Analysis of gene network robustness based on saturated fixed point attractors

The analysis of gene network robustness to noise and mutation is important for fundamental and practical reasons. Robustness refers to the stability of the equilibrium expression state of a gene network to variations of the initial expression state and network topology. Numerical simulation of these variations is commonly used for the assessment of robustness. Since there exists a great number of possible gene network topologies and initial states, even millions of simulations may be still too small to give reliable results. When the initial and equilibrium expression states are restricted to being saturated (i.e., their elements can only take values 1 or −1 corresponding to maximum activation and maximum repression of genes), an analytical gene network robustness assessment is possible. We present this analytical treatment based on determination of the saturated fixed point attractors for sigmoidal function models. The analysis can determine (a) for a given network, which and how many saturated equilibrium states exist and which and how many saturated initial states converge to each of these saturated equilibrium states and (b) for a given saturated equilibrium state or a given pair of saturated equilibrium and initial states, which and how many gene networks, referred to as viable, share this saturated equilibrium state or the pair of saturated equilibrium and initial states. We also show that the viable networks sharing a given saturated equilibrium state must follow certain patterns. These capabilities of the analytical treatment make it possible to properly define and accurately determine robustness to noise and mutation for gene networks. Previous network research conclusions drawn from performing millions of simulations follow directly from the results of our analytical treatment. Furthermore, the analytical results provide criteria for the identification of model validity and suggest modified models of gene network dynamics. The yeast cell-cycle network is used as an illustration of the practical application of this analytical treatment.

: The cell-cycle network of the budding yeast (A) and simplified cell-cycle network with only one checkpoint "cell size" (B) (reproduced from the above paper).
The dynamics of the network was defined by Li et al. as where 1 and 0 correspond to active and inactive states of the gene, i.e., 0 instead of -1 is used to represent inactive state, and a ij is similar to w ij in dynamics − 1, (i = 1, 2, . . . , n), a ij takes values 1 (green arrow) and -1 (red arrow), respectively. Moreover, the yellow loop is used to represent "self-degeneration" (value -1).
Using model (1) Li, et al. found that there exit 7 saturated fixed point attractors (upon considering 0 as -1) and all of the 2 11 = 2, 048 possible saturated initial expression states converge to one of the seven fixed point attractors (see Table 1).  Note that dynamics (1) is different from dynamics (2). For j a ij S j (t) = 0, dynamics (2) gives S i (t + 1) = 0, not S i (t). Dynamics (2) with the W constructed directly from the connectivities in Figure 1(B) will not give the same result as that from dynamics (1). Hereafter, 0 representing an inactive state by Li et al. will be replaced by -1.
All the information given by the simplified model for the yeast cell-cycle network (Figure 1(B)) will be considered as "experimental information" for the budding yeast. The procedure shown below is only an illustration for network construction and robustness analysis by our analytical treatment combining the "experimental information" of budding yeast. A biologist familiar with the yeast cell-cycle network possibly may construct better networks W .

Network construction
Define the node order from 1 to 11 as that given in Table 1, i.e., Cln3 is node 1, MBF is node 2, etc. We first construct all viable networks sharing the most stable saturated equilibrium expression state, the first fixed point attractor in Table 1 S for dynamics (2). According to Theorem 2, these viable networks will also share the other saturated equilibrium expression state Theorem 6 gives the criterion to construct all of such networks W , i.e., for the given saturated equilibrium state S 1 the following inequalities must be satisfied by W : under the condition w ij ∈ [−a, a]. Here, 2, 3, 4, 6, 7, 8, 10, 11}. (8) When w ij can take any value within [−a, a], there is an infinite number of solutions for (5,6). As mentioned by Li, et al., "the overall dynamic properties of the network are not very sensitive to the choice of these parameters" (w ij ), but the connectivity patterns of the network, i.e., the regulatory influence between genes (activation, repression and absence) is important for gene network robustness analysis. Therefore, we restrict that w ij can only take the discrete values 1 (activation), -1 (repression) and 0 (absence).
When w ij only takes values [−1, 0, 1], to satisfy (5,6), then each row of W must have 5 or more nonzero elements due to β ≥ 5. Otherwise, the network would not have a saturated equilibrium state. This problem occurs not only for networks with less than 5 genes, but also for larger networks with sparse connectivities between genes. For example, Node 1 (Cln3) in Figure 1(B) is a pure "parent" node, which does not have any regulation coming from all other "children" nodes, i.e., all w 1j = 0 for j = 1, and for S 1 the condition (6) does not hold: To avoid this problem, the factor β is introduced such that so to satisfy the condition (5,6),ŵ ij can only take values [−1, 0, 1] without any restriction on the number of nonzero elements in each row ofŴ . For the sake of notational simplicity, in the sequel we still use W instead ofŴ , but write dynamics (2) as The necessary and sufficient condition for S to be a saturated equilibrium state for dynamics (11) then become For saturated equilibrium S 1 , (12,13) may be written as − 11 j=1,j =5,9 or 11 j=1,j =5,9 11 j=1,j =5,9 Condition (16,17) will be used to construct all viable networks W for dynamics (11) sharing saturated equilibrium state S 1 . When w ij takes values [-1,0,1], for all 9 possible combinations of (w i5 , w i9 ) we can determine all permitted row patterns satisfying (16,17) for the other nine w ij 's in each combination of (w i5 , w i9 ).
To avoid −S and 0 as an equilibrium state and fixed point, similarly, the modified dynamics is proposed. The necessary and sufficient condition to have S as a saturated equilibrium state for dynamics (18) is and for S 1 (19,20) becomes or 11 j=1,j =5,9 Condition (23,24) can be used to construct all viable networks W for dynamics (18) with a given set of θ i 's sharing saturated equilibrium state S 1 . Since there is no unambiguous biological interpretation of the values of θ i as [−1, 0, 1] for w ij , we will not construct all such viable networks here.
The summation on the lefthand side is a negative number ≤ −3. To satisfy this condition, three to nine w ij (j = 5, 9) can be -1. All permitted row patterns are counted in Table 2. Table 2. Permitted row patterns of w ij (j = 5, 9) for w i5 + w i9 = −2 Number of -1 Number of 0 Number of 1 Number of patterns Formula Number To satisfy this condition, two to nine w ij can be -1. All permitted row patterns are given in Table 3. Table 3. Permitted row patterns of w ij (j = 5, 9) for w i5 + w i9 = −1 Number of -1 Number of 0 Number of 1 Number of patterns Formula Number To satisfy this condition, there may be one to nine w ij being -1. All permitted row patterns are given in Table 4. Table 4. Permitted row patterns of w ij (j = 5, 9) for w i5 + w i9 = 0 Number of -1 Number of 0 Number of 1 Number of patterns Formula Number According to (16) we have 11 j=1,j =5,9 w ij ≤ 0.
To satisfy this condition, nine w ij may be all 0, or one to nine w ij may be -1. All permitted row patterns are given in Table 5. Table 5. Permitted row patterns of w ij (j = 5, 9) for w i5 + w i9 = 1 Number of -1 Number of 0 Number of 1 Number of patterns Formula Number 18 9 0 0 C 9 9 1 sum 10,907 1.1.5 Determination of w ij (j = 5, 9) with w i5 + w i9 = 2 According to (16) we have 11 j=1,j =5,9 w ij ≤ 1.
1.2 Determination of w ij for i = 5, 9 in dynamics (11) We then determine all permitted row patterns of w ij for i ∈ J − (S 1 ), i.e., i = 5, 9. For i = 5, 9, inequality (15) − j=1,j =5,9 needs to be satisfied. Condition (15) can be rewritten as which is the same as the inequality (16) used to determine w ij for i = 5, 9, except that w ij is replaces by −w ij . Therefore, all the permitted row patterns of w ij for i = 5, 9 multiplied by -1 are the permitted row patterns for i = 5, 9. Hence, the total number of permitted row patterns for each i(i = 1 − 11) is the same, 72, 219. Then, the total number of viable networks for dynamics (11) with n = 11 and sharing the saturated equilibrium state S 1 (and S 2 ) is 72, 219 11 ≈ 2.7872 × 10 53 .
As shown below, according to Figure 1 There are many choices of networks consistent with the experimental observations.

Construction of the yeast cell-cycle network W with dynamics (11)
According to the definition for the green and red arrows as well as the yellow loop given by Li et al., the network W 0 directly constructed from Figure 1(B) is W 0 does not satisfy condition (16,17) for any saturated state and does not have a saturated equilibrium state for dynamics (11). However, W 0 will be used as a starting point (considered as an experimental observation of the connections of the network) to construct networks with the saturated equilibrium expression state S 1 . The construction of networks reduces to satisfying condition (16, 17) or (23, 24) consistent as much as possible with the connectivities in the experimental observation.
There are three choices: 1. Suppose that the information (w 1j (j = 1) are all zero) given in Figure  1(B) is correct. Then the only choice is w 11 = 1.
2. If a biologist has observation that there exist activating regulations from other nodes to Clin3, then some appropriate w 1j (j = 5, 9) can be set to 1.

Row 2 (i = 2 ∈ J − (S 1 ))
For i = 2 ∈ J − (S 1 ), we also have (w 25 , w 29 ) = (0, 0) in W 0 and 11 j=1,j =5,9 w 2j ≥ 1 should be satisfied. The summation on the lefthand side is zero for row 2 of W 0 . So we need to add at least one 1 in any position j = 5, 9. Without any additional information about activation for node 2, we add 1 in j = 2, i.e., a "self activation regulation". The other way is to set w 21 to be 2, but this is beyond our restriction w ij ∈ [−1, 0, 1]. If a biologist has sufficient information to prove such a choice, it may be adopted. Therefore, we choose row 2 as ( 1 1 0 0 0 0 0 0 0 −1 0 ).

Row 4 (i
Similar to rows 2 and 3, we need to add 1 in some j = 5, 9 in row 4 of W 0 . As we do not have any additional information, we still add 1 in j = 4 which cancels -1 in W 0 to give row 4 as w 6j ≥ 1 should be satisfied, but it is 0 for row 6 in W 0 . We need to add at least 1 in j = 5, 9. Without any additional information, we add 1 in j = 6 which cancels -1 there to give w 8j ≥ 1 + w 85 + w 89 = 0, but it is -1 for row 8 in W 0 . We need to add at least one 1 in j = 5, 9. Comparing Figure 1(B) with 1(A), we see that node 8 is a combination of node Clb5,6 and node DNA Replication. There is a green arrow between the two nodes which was ignored after combination. To represent this green arrow we can add 1 at j = 8 as "auto activation regulation" to give one choice as ( 0 1 0 0 0 0 −1 1 −1 0 0 ).
Moreover, there are bi-activating regulations between nodes Mcm1/SFF and Clb1,2, but only one green arrow from node Clb5,6 to node Clb1,2. There is a possibility to have a green arrow back from node Clb1,2 to node Clb5,6. So we may have another choice for row 8 as w 10,j ≥ 1 + w 10,5 + w 10,9 = −1.
Combining all choices of the 11 rows, we finally obtain two yeast cell-cycle networks for dynamics (11) and W 1 and W 2 differ only for w 8,10 .
For n = 11, there are 2 11 = 2, 048 saturated states. All of the 2,048 states were tested by condition (40) for W 1 and W 2 , respectively, to determine which of them are saturated equilibrium states of W 1 and W 2 . The test for 2,048 states took only 0.01 seconds by Matlab on a Dell Precision Workstation T3400. The saturated equilibrium expression states for a given network W with threshold vector Θ in dynamics (18) can be determined by using the condition S i n j=1 w ij S j − θ i ≥ 1, (i = 1, 2, . . . , n).
The resultant saturated equilibrium expression states for W 1 , W 2 and W 0 with Θ have been obtained as shown below.
1. W 1 W 1 has two saturated equilibrium expression states for dynamics (11)

W 2
W 2 has four saturated equilibrium expression states for dynamics (11) The S 1 and S 3 are just the 1st and 3rd fixed point attractors in Table 1.
3. W 0 with Θ For W 0 with parameters Θ, there is only a single saturated equilibrium state for dynamics (18)

Robustness to noise
First, the numbers of saturated initial expression states converging to each equilibrium expression state for W 1 , W 2 are determined by either directly solving the dynamics (11) or using the modified condition of Theorem 5 For W 1 , the CPU times are 0.8 and 0.3 seconds, respectively to check all 2,048 saturated states, i.e., using Theorem 5 the CPU time is ∼ 41% of that for direct solving the sigmoidal function. The results are given in Table 7. The robustness to noise R nt is defined as where m is the total number of fixed points. The robustness to noise R n i of a given saturated equilibrium expression state S i for a gene network W is specified by the ratio of the number N i of saturated initial expression states converging to S i , with respect to the total number 2 n of possible saturated initial states R n i = N i /2 n .
Note that for W 1 , W 2 , no saturated initial state converges to the unstable fixed point 0. Therefore, in the calculation of R nt , we ignore 0 and only consider the saturated equilibrium states. The robustness to noise R nt and R n i are given in Table 8. There are significant differences between R n i (i = 1, 2, 3, 4) for W 2 . Obviously, the saturated equilibrium states S 1 , S 2 are much more stable than S 3 , S 4 . The robustness to noise R n ij of a given viable pair of saturated equilibrium and initial expression states S i and S j (0) for a gene network W is specified by the ratio of the number N ij of neighboring saturated initial expression states differing from S j (0) by only one element and still converging to S i , with respect to the total number n of possible one element differing saturated initial states R n ij = N ij /n.
The robustness to noise R n ij for each viable pair of saturated equilibrium and initial expression states was calculated. The distribution of R n ij , i.e., how many pairs with the same value of R n ij , is given in Tables 9 and 10. Table 9. The distribution of R n ij for the networks W 0 with Θ and W 1 Final state R n ij (for W 0 with Θ) R n ij (for W 1 ) 11/11 10/11 S 1 2,048 1,024 S 2 1,024 The results show that W 0 with Θ is completely stable for any viable pair; for W 1 , there is one neighbour of S j (0) differing at the first element, which causes a change in the saturated equilibrium state S i ; for W 2 , the distribution of R n ij is divergent, and S 1 and S 2 are much more stable than S 3 and S 4 .

Robustness to mutation
The robustness to mutation R m i for saturated equilibrium state S i is defined as where N v W and N v w i respectively are the total numbers of viable (keeping the saturated equilibrium state unchanged) single w ij changes of W and its ith row; N W (= 242) is the total number of possible single w ij changes of W . R m i 's have been calculated for W 1 , W 2 and W 0 with the Θ given above as shown in Table 11. Note that for S 1 the robustness to mutation R m i is almost the same for W 1 , W 2 and W 0 with Θ. Robustness to mutation R m ij for a specified pair of saturated equilibrium and initial states S i and S j (0) of a viable network W is defined as where N v m ij is the total number of viable single w ij changes of W with respect to a specified states S i and S j (0). N v m ij can be obtained by determining how many networks in N v W , which share the saturated equilibrium state S i , also share the saturated initial expression state S j (0). The resultant distribution, i.e., how many viable pairs having the same R m ij for W 1 , W 2 and W 0 with the Θ given above is shown in Table 12