Analysis of gene network robustness based on saturated fixed point attractors

EURASIP Journal on Bioinformatics and Systems Biology Analysis of gene network robustness based on saturated fixed point attractors Genyuan Li 0 Herschel Rabitz 0 0 Department of Chemistry, Princeton University , Princeton, NJ 08544 , USA 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. Robustness of gene networks; Fixed point attractor; Sigmoidal function; Yeast cell-cycle network 1 Introduction A subset of genes in a cell whose protein products mutually regulate one another’s expression at the transcriptional level will be referred to as a ‘gene network’ in this paper. The concentration of proteins encoded by the genes changes in time due to auto- and cross-regulation by the gene products. Each of such network is considered as a dynamical system. Gene networks must be robust with respect to ever-changing environments. Robustness here refers to the ability of a gene network to respond to shortterm changes in the environment and quickly return to its functional steady state. Moreover, a gene network itself may endure small structural changes and mutations, while still retaining its desired steady state. The robustness of gene networks depends on their topology with some networks being more stable than others. The analysis of the relationship between the topology of a gene network and its robustness to environmental and structural perturbations is important both theoretically and practically [1-10]. Recently, Wagner and coworkers considered the robustness of gene networks [11-14], and a similar assessment was given by Cho et al. [15] In these works, a simplified model proposed by Wagner was used to describe the dynamics of the gene expression states [11,12]. Let G = (G1, . . . , Gn) represent the n genes in a network. The concentration of proteins encoded by the genes (G1, . . . , Gn) is denoted by P = (P1, . . . , Pn). For computational convenience, the admissible concentration range for each Pi is normalized and restricted to the interval [ 0, 1], where Pi = 1 corresponds to the maximum possible concentration, i.e., the corresponding gene Gi is in a state of maximum transcriptional activation. It is also assumed that Pi = 0.5 means that gene i is 50% ‘on’. The dynamics of the expression states of the genes in a network is often described by some sigmoidal function gc(x) Pi(t + τ ) = gc ⎝ ⎛ n j=1 ⎞ wijPj(t)⎠ , (i = 1, 2, . . . , n) ( 1 ) where τ is a time constant characteristic of the process under consideration. In some work, τ was set to be 1. The constant wij ∈ describes the strength of interaction (i.e., transcriptional regulation) of the product of gene j with gene i, i.e., the degree of transcriptional activation (wij > 0), repression (wij < 0), or absence (wij = 0). These constants define a matrix of connectivities W = (wij) within the network. To facilitate the analytical treatment, the variable transformation S = 2P − (1 . . . 1)T is employed. Using the sigmoidal function σ proposed by Siegal [16] and Cho [15] for gc, ( 1 ) becomes ⎛ n Si(t + τ ) = σ ⎝ j=1 = wijSj(t)⎠ ⎞ 2 n j=1 to describe the dynamics for each element of primitive objects vi (cells, nuclei, fibers, and synapses), where ga is a sigmoidal threshold function, T ab is similar to wij in ( 3 ); via, vib denote the elements of vector vi; ha determines the threshold of ga. The long-time behavior of this system has been studied, and, in some cases, is controlled by a simple limit set with Si ∈[ −1, 1], and Si = 0 corresponding to 50% of gene i being ‘on’. Notwithstanding the simplicity of ( 3 ), variants of this model have been successfully used to study (a) the robustness of gene regulatory networks [12,16,17], (b) the role of robustness in evolutionary innovation [18,19], and (c) how recombination can produce negative epistasis [20]. Mjolsness at al. [21] proposed a model τa ddvtia = ga n b=1 Tabvib + ha − λavia, (a = 1, 2, . . . , n) ( 2 ) ( 3 ) ( 4 ) ( 5 ) which is similar to ( 3 ) with the additional parameters λa and ha. This model has been successfully applied to treat the blastoderm of Drosophila melanogaster [21]. Similarly, Mendoza and Alvarez-Buylla [22] used a model ⎞ xi(t + 1) = H ⎝ wijxj(t) − θi⎠ , ⎛ n j=1 where H is the Heaviside step function H(x) = 1, if x > 0, 0, if x ≤ 0 ( 6 ) ( 7 ) to describe the dynamics of a genetic regulatory network for Arabidopsis thaliana flower morphogenesis. This model is also similar to ( 3 ) except that the sigmoidal function σ is replaced by the Heaviside step function and a threshold parameter θi is included. All these models present simplified descriptions of gene network dynamics. Nevertheless, the models are still useful for obtaining insights into the dynamics of gene networks. In the following analysis for gene network robustness, we will employ the sigmoidal function model in ( 3 ), and its modification with threshold parameters. The robustness of a gene network specified by W to noise (environmental) and mutation (structural) perturbations may be expressed as the stability of the final equilibrium (or steady) expression state S(∞) obtained from the solution of ( 3 ), respectively, to changes of initial expression state S(0) and to changes of W. A complete and reliable robustness analysis would seem to call for an exhaustive sampling over the space of possible initial states for a given network and then repeating the same simulations for all possible networks. This task is infeasible as there are many possible networks W, and each W may have many initial/final equilibrium expression states. Consider a simple case where Si(0) and Si(∞) (i = 1, 2, . . . , n) can only take values −1,1 and wij can only take values −1, 0, 1. In this case, there are 2n possible initial/final states and 3n2 possible gene networks. For a modest network of size n = 20, there are 220 = 1, 048, 576 initial/final expression states and 3400 ≈ 7 × 10190 gene networks. Even if one arbitrarily makes the restriction that 75% of the wij interactions are zero, and that the remaining 25% of the wij interactions can only take nonzero values −1,1 to reduce the possible number of networks 100 [13], for n = 20 there are still C400 ×2100 possible gene networks. Further restriction may also be applied to reduce the possible number of initial expression states [13]. Even under these restrictions, solving ( 3 ) for all possible initial expression states and gene networks is still infeasible, and only a small fraction of them can be randomly sampled and simulated. Such limitations leave open the reliability of the conclusions obtained from the simulations. Previous work [13] concerned networks with connectivity W whose expression dynamics start from a prespecified initial state S(0) at some time t = 0 and arrive at a prespecified stable equilibrium or ‘target’ expression state S(∞); these networks are referred to as ‘viable’. Then, the values of some elements of S(0) or W were changed for each viable network to check whether S(∞) is reached. These studies entailed performing millions of simulations with different network topologies and initial expression states. Although the number of simulations is much smaller than 2n and 3n2 for n = 20, some specific conclusions were obtained. First, the fraction of viable networks, that is, networks that arrive at a prespecified target expression state S(∞) given an initial gene expression state S(0) to the total number of possible networks, is generally very small. For moderately sized networks of n = 20 genes (with the number M of nonzero wij set to be 200, and the fraction d of elements different between S(0) and S(∞) set to 0.5), the fraction of viable networks was found to be vf = 5.1 × 10−9 ± 1.7 × 10−10. Due to the large numbers 2n and 3n2 , the qualitative correctness of this conclusion is clear. Since there are 2n possible equilibrium states, even if we only consider the factor of expression states, the probability that a network W arrives at a prespecified S(∞) is expected on the order of 1/2n. The viable networks in this prior work could be organized as a graph with each node corresponding to a network of a given topology, and two nodes are connected by an edge if they differ by a single regulatory interaction (i.e., they differ in one element wij). Remarkably, this graph is connected and can be easily traversed by gradual changes of the network topology. Thus, highly robust topologies can evolve from topologies with low robustness through gradual Darwinian topological changes. These results are claimed to be valid for discrete and continuous wij taking values [ −1, 0, 1] and over the interval [ −a, a], respectively. While simulations are valuable, they do not provide a complete picture. Ciliberti et al. [13] considered the case where each element of the initial and equilibrium expression states, S(0) and S(∞), can only take the values 1 or −1 corresponding to maximum and minimum possible protein concentrations. We call them saturated expression states. Under this condition, the present paper provides an analytical robustness assessment of gene networks whose dynamics can be described by ( 3 ) and its modification with threshold parameters. This 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; (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 identification of model validity and suggest modified models of gene network dynamics. The paper is organized as follows: Section 2 first defines the saturated state and saturated fixed point attractor for dynamics ( 3 ), and then gives the necessary and sufficient condition for a gene network to have a given saturated equilibrium state. Sections 3 and 4 analyze the robustness to noise and mutation, respectively. Section 5 proposes a modification of dynamics ( 3 ) with threshold parameters. Section 6 gives an illustration of the practical application of this analytical treatment: the model construction of the yeast cell-cycle network and its robustness assessment. The details of the treatment of the yeast cell-cycle network are given in an Additional file 1: Supplementary information. Finally, Section 7 presents conclusions. Mathematical proofs of the theorems in the main text are given in the Appendix. 2 Saturated states and fixed point attractors In this work, the initial and equilibrium expression states Si(0) and Si(∞) for gene i can only be either active (Si(0) and Si(∞) = 1) or inactive (Si(0) and Si(∞) = −1) [11,12,14]. The initial and equilibrium expression states with Si(0) and Si(∞) = ±1 are referred to as saturated initial and equilibrium expression states. Under the condition that the initial and equilibrium expression states S(0) and S(∞) are saturated, we may analyze the robustness and evolvability of gene networks analytically, as explained below. 2.1 Saturated sigmoidal function A continuous function f (x) defined on satisfying (f1) f (x) = 1 if x ≥ 1, f (x) = −1 if x ≤ −1, (f2) f (x) is a strictly increasing and continuous function for x ∈[ −1, 1] and f (0) = 0, is called a saturated sigmoidal function. Furthermore, if (f3) f (x) ≥ x for x ∈ (0, 1] and f (x) ≤ x for x ∈[ −1, 0), we call f (x)] a dissipative saturated sigmoidal function. Note that for the particular sigmoidal function σβ (x) in domain [-1, 1] 2 x(t + τ ) = σβ (x(t)) = 1 + e−βx(t) − 1, ( 8 ) the conditions (f1) to (f3) are approximately satisfied when β is sufficiently large. For example, when β = 5 and 10, we have σβ ( 1 ) = 0.9866 and 0.9999 (approximately 1); and σβ (−1) = −0.9866 and −0.9999 (approximately −1), respectively. Therefore, in numerical simulation, σβ (x) can be considered as a dissipative saturated sigmoidal function for a sufficiently large β. We refer to β = 5 and 10 as having 0.99 and 0.9999 confidence levels for the dissipative saturated sigmoidal function ( 8 ), because |σβ (±1)| is equal to 0.99 and 0.9999, respectively (see Figure 1). In the sequel, we set β ≥ 5. 2.2 Necessary and sufficient condition for a saturated state S to be an equilibrium state or a fixed point attractor of dynamics (3) with a given W When the equilibrium states S(∞) are saturated, the analysis of robustness and evolvability for gene networks can be readily performed by utilizing Feng and Tirozzi’s treatment for neural networks [23]. Definition 1. A saturated state S in [ −1, 1]n is called a saturated fixed point attractor (or saturated equilibrium state) of dynamics ( 3 ), if there exists a nonempty neighborhood B(S) of S such that lim S(t) = S t→∞ for S(0) ∈ B(S) and n j=1 wijSj = 0, for all i. ( 9 ) ( 10 ) ( 11 ) It is easy to see that a saturated state S is a saturated fixed point of dynamics ( 3 ) with 0.99 (β = 5) or 0.9999 (β = 10) confidence level when and n n j=1 j=1 wij sign(Sj) ≥ β, if Si = 1, wij sign(Sj) ≤ −β, if Si = −1, i.e., S is a saturated fixed point with 0.99 (β = 5) or 0.9999 (β = 10) confidence level. When n j=1 wij sign(Sj) < 5, for any i, then −0.99 < Si(t + τ ) < 0.99 for any t, and the gene network cannot have a saturated fixed point. Let S be a saturated state. We define J+(S) = {i, Si = 1}, J−(S) = {i, Si = −1} ( 12 ) which denote, respectively, the two sets of all integers in {1, 2, . . . , n} with Si = 1 and Si = −1. Figure 1 The function σ β (x). 2 1 −2 −2 y=x β=1 β=2 β=5 β=10 0 x −1 1 2 Theorem 1. The necessary and sufficient condition for a saturated state S to be an equilibrium expression state or a fixed point attractor of dynamics ( 3 ) with a given matrix W is or equivalently Proof. If a saturated state S is a fixed point of ( 3 ), it must satisfy ( 9,10 ), i.e., ( 13,14 ). If ( 13,14 ) are satisfied, so are ( 9,10 ), then ( 11 ) is satisfied and S is a fixed point. A saturated state S satisfying ( 13,14 ) (or equivalently (15)) is a fixed point attractor. The proof is given in Theorem A2 in the Appendix . There are a total 2n saturated states S. Since ( 15 ) only involves simple multiplication and summation, for a modest n one can test all 2n saturated vectors S for a given W and find all of its saturated fixed point attractors without iteratively solving the sigmoidal function ( 3 ). For n = 11 and 211 = 2, 048, the test takes 0.01 s by Matlab on a Dell Precision Workstation T3400. Example 1. Consider the network model proposed by Azevedo et al. [20] for the gap gene system of Drosophila melanogaster shown in Figure 2. The authors obtained the equilibrium state S1 = (−11 − 11) by iteratively solving a sigmoidal function similar to ( 3 ). Using ( 9,10 ), we can determine this solution by noting ⎡ 2 −1 −1 −5 ⎤ ⎡ −1 ⎤ ⎡ −7 ⎤ W · S1 = ⎢⎢⎣ −−18 −11 01 −14 ⎥⎦⎥ ⎢⎢⎣ −11 ⎥⎦⎥ = ⎢⎢⎣ −150 ⎥⎦⎥ 0 −1 −6 1 1 6 which shows that the necessary and sufficient condition is satisfied with 0.99 confidence level 4 j=1 4 j=1 wijSj ≥ 5, i = 2, 4 ∈ J+(S1), wijSj ≤ −5, i = 1, 3 ∈ J−(S1). There are 24 = 16 saturated states. Similar tests for the other 15 saturated states were performed. The case S2 = −S1 = (1 − 11 − 1) is the only other saturated equilibrium state for the network. 3 Robustness to noise Robustness to noise may be assessed for (a) each saturated equilibrium expression state or (b) a specified pair of saturated equilibrium and initial expression states. In case 1, we need to compare how many of 2n possible saturated initial expression states converge to each saturated equilibrium expression state; in case 2, we need to determine how many neighbours (differing only in one element) of the S(0) converge to the same saturated equilibrium expression state for a given W [13]. In either case, we need to establish the condition under which a saturated initial expression state converges to a given saturated equilibrium expression state of W. 3.1 Relationship between saturated initial expression states and saturated equilibrium expression states Theorem 2. If S is a saturated equilibrium expression state (or a fixed point attractor) of dynamics ( 3 ) with a given W, so is −S. Proof. Since S is a saturated fixed point attractor of ( 3 ) with a given W, then ( 15 ) is satisfied. Multiplying by −1 within and outside the parentheses in ( 15 ) will not change its right-hand side: ⎛ (−1)Si ⎝ (−1) ⎛ (−Si) ⎝ n j=1 n j=1 wijSj⎠ ⎞ ⎞ ≥ β, (i = 1, 2, . . . , n), wij(−Sj)⎠ ≥ β, (i = 1, 2, . . . , n) which proves that −S is also a saturated fixed point attractor, i.e., a saturated equilibrium expression state for W. Example 1 demonstrates its validity. Corollary 1. Dynamics ( 3 ) either does not have a saturated equilibrium expression state, or has an even number of saturated equilibrium expression states. This result can be obtained immediately from Theorem 2. Example 2. Consider a gene network given by ⎡ 0 1 1 −1 −1 −1 ⎤ ⎢⎢ 11 10 01 −−11 −−11 −−11 ⎥⎥⎥ W = ⎢⎣⎢⎢⎢ −−11 −−11 −−11 10 01 11 ⎦⎥⎥⎥ −1 −1 −1 1 1 0 ( 16 ) ( 17 ) with the discrete values [-1, 0, 1] for wij, and there is no auto-regulation (wii = 0). The six genes are separated into two groups {1, 2, 3} and {4, 5, 6}. The regulations are activating within each group, but repressing between the two groups. For such a simple system, it is easy to find by observation that amongst the 26 = 64 saturated states only two states S1 = ( 1 1 1 −1 −1 −1 ), S2 = ( −1 −1 −1 1 1 1 ) satisfy the necessary and sufficient condition ( 15 ) to be saturated equilibrium expression states with 0.99 confidence level. An examination of ( 15 ) for all 64 saturated states proved this to be the case. Moreover, S2 = −S1 satisfies Theorem 2 and Corollary 1. Theorem 3. If S(t) converges to a saturated equilibrium expression state S, then −S(t) converges to −S. Proof. See Theorem A3 in the Appendix. Example 3. The gene network given in ( 17 ) is used to show the validity of Theorem 3. The following two initial saturated states S1(0) = ( 1 1 −1 1 1 1 ), S2(0) = ( −1 −1 1 −1 −1 −1 ), with S2(0) = −S1(0), were used for dynamics ( 3 ). The two solution trajectories are found to satisfy Theorem 3. Corollary 2. In the hypercube [ −1, 1]n, the volume of the region where points converge to a saturated equilibrium state S is equal to the volume of the region where points converge to −S. Proof. Considering that S(t) and −S(t) are symmetric and have the same distance to the origin (see Figure 3), then in [ −1, 1]n (which is symmetric to the origin) the volume of the region where points converge to S is equal to the volume of the region where points converge to −S. However, if S1 and S2 are two saturated fixed point attractors, but S2 = −S1, for a given W, generally in [ −1, 1]n the volume of the region where points converge to S1 may not be equal to the volume of the region where points converge to S2. Theorem 4. S = 0 is a fixed point, but may not be a fixed point attractor for a W having a saturated fixed point attractor S. Proof. S = 0 is a fixed point because Initial state ( 1 1 −1 1 1 1 ) Let B = {Sˆ ∈ n | Sˆ < } be an open ball of radius centered at the origin in n where is chosen sufficiently small such that there is only a single fixed point 0 within B . The sufficient condition for 0 to be a unique fixed point attractor is (see (110) of Theorem A1 in the Appendix) |wij| ≤ 2, the following condition is satisfied after a finite number k of iteration steps of ( 3 ) starting from S(0) wijSj(t ≥ k) > −ln (αi −1)− (αi −1)2 − 1 , i ∈ J+(S), A saturated fixed point attractor S of W must satisfy ( 15 ). Then, we have wijSj(t ≥ k) < ln (αi − 1) − (αi − 1)2 − 1 , i ∈ J−(S), ≥ β > 2, (i = 1, 2, . . . , n). αi = |wij| n j=1 Therefore, 0 may not be a fixed point attractor, but an unstable fixed point. For an unstable fixed point, there exist divergent or both convergent and divergent neighborhoods of 0. In the convergent region, trajectories will be attracted to 0; In the divergent region, trajectories will leave from the neighbourhood of 0. Therefore, it is possible upon starting from a saturated initial expression state S(0) that the solution trajectory arising from the sigmoidal function ( 3 ) may converge to 0, or first approach 0, but then enter the divergent region, and move away from that neighborhood. under the constraint that αi ≥ 2. Proof. See Theorem A4 in the Appendix . Theorem 5 provides a way to determine all possible saturated initial expression states converging to a saturated equilibrium expression state of a W. The constraint αi ≥ 2 implies that the gene network must contain more than one gene if wij only takes values in [−1, 0, 1]. This is always true. Theorem 5. A saturated initial expression state S(0) converges to a saturated equilibrium expression state S if Example 4. For the gene network given in ( 17 ) all αi = 5, and the sufficient condition for a saturated initial ( 21 ) ( 22 ) ( 23 ) expression state converging to a given saturated equilibrium expression state S is then n where S(t ≥ k) is the solution of ( 3 ) after k iterations starting from a saturated state S(0). All 26 = 64 saturated states are used as initial expression states for dynamics ( 3 ). The results are given in Table 1. Using the condition in ( 24,25 ), we found that each saturated equilibrium expression state for the gene network given in (17) has 22 saturated initial expression states (including the saturated equilibrium expression state itself ). The remaining 20 saturated initial states, which do not satisfy ( 24,25 ), converge to 0. For saturated initial expression states converging to one of the two saturated equilibrium expression states, the condition in ( 24,25 ) is satisfied starting from the iteration number k as either 0 or 1. For saturated initial expression states converging to 0, the condition given in ( 24,25 ) is never satisfied. Therefore, in this gene network, the sigmoidal function ( 3 ) either does not need to be solved, or only needs to be solved just once, to determine which and how many saturated initial expression states converge to a particular saturated equilibrium expression state. The computational effort will be reduced. For a network with 11 genes, using Theorem 5 has approximately 60% CPU time saving compared to completely solving ( 3 ). Figure 4 gives the projection of a trajectory converging to the final state 0 onto the two-dimensional (S1, S5)subspace. Since 0 may be an unstable fixed point, and the values of the elements of S(t) are not exactly zero, but are within a small region around 0. Therefore, the trajectory is sensitive to computational precision, that determines which region the trajectory enters. Continued iteration may lead to the trajectory entering the diverging region of 0 and leaving away from 0. When there are no limit cycles for the W, the trajectory must converge to one of the two saturated equilibrium expression states. The outcome depends on the precision used in the computation. Table 1 Number of saturated initial states converging to different final states for the gene network given in ( 17 ) Final state (1 1 1 −1 −1 −1) (−1 −1 −1 1 1 1) (0 0 0 0 0 0) Number of saturated initial states 22 22 Figure 5 shows two trajectories that first approach to 0 and then leave the neighborhood of 0 and converge to one of the two saturated equilibrium expression states (1 1 1 −1 −1 −1) and (−1 −1 −1 1 1 1), respectively. Such a property of 0 for dynamics ( 3 ) may not be meaningful biologically and causes confusion. This problem will be discussed in Section 5. 3.2 Definition for robustness to noise Robustness to noise may be defined as follows for a network W, for each of its saturated equilibrium expression states, and for a specified pair of saturated equilibrium and initial expression states, respectively. Definition 2. The robustness to noise Rnt of a given gene network W is specified as inversely proportional to the total number m of fixed points. Rnt = 0, if m = 0, 1/m, otherwise. (26) The subscript nt denotes ‘noise’ and ‘total’. According to the definition, larger values of m correspond to worse robustness to noise. This definition is reasonable because more fixed points a W has, then less saturated initial states converging to each of its equilibrium states. Changing an initial state has more chances to cause change of the equilibrium state. In this definition, m should also include limit cycles. Since there is no simple way (like the Banach fixed point theorem for the existence of fixed point attractor) to determine the existence of limit cycle for sigmoidal functions, in the current work, we are unable to include it. The robustness Rnt takes on values in [ 0, 1]. For W given in Example 2, there are three fixed points: 0, (1 1 1 − 1 −1 −1) and (−1 −1 −1 1 1 1). Therefore, Rnt = 1/3, which implies that there is a 2/3 probability for having a saturated equilibrium expression state change caused by a saturated initial expression state change. Definition 3. The robustness to noise Rni of a given saturated equilibrium expression state Si for a gene network W is specified by the ratio of the number Ni of saturated initial expression states converging to Si, to the total number 2n of possible saturated initial states Rni = Ni/2n. (27) The subscripts n and i denote ‘noise’ and the ith saturated equilibrium expression state. For saturated equilibrium expression states (1 1 1 −1 −1 −1) and (−1 −1 −1 1 1 1) in Example 2, Rni = 22/64 ≈ 1/3 (i = 1, 2). Ciliberti et al. [13] argued that the pair of saturated equilibrium and initial expression states, S(∞) and S(0), play 0.8 0.6 0.4 0.2 0 −0.2 −0.4 −0.6 −0.8 −1 0 Initial state (−1 −1 1 −1 1 −1) S1 S5 20 40 60 80 100 120 Iteration number, t a central role for a viable network, but the variation of initial expression state in realistic cases is often mild. They define one measure of robustness to noise as the probability Rv,I that a change in one gene’s expression state in the saturated initial expression state S(0) leaves the unchanged network’s saturated equilibrium state S(∞). Following this pattern, we also define the measure of robustness to noise for a given pair of saturated equilibrium and initial expression states with a viable network as follows. Definition 4. The robustness to noise Rnij of a given pair of saturated equilibrium and initial expression states Si and Sj(0) for a viable gene network W is specified by the ratio of the number Nij of neighboring saturated initial expression states differing from Sj(0) by only one element and still converging to Si, with respect to the total number n of possible one element differing saturated initial states Rnij = Nij/n. (28) Initial state (1 −1 −1 1 −1 −1) Initial state (−1 1 1 −1 1 1) 1 Example 5. Determination of Rnij for the network in Example 2. The network ( 17 ) has two equilibrium states S1, S2 with 22 saturated initial states Sj(0) converging to each. The Rnij for each pair of Si and Sj(0) was determined from solving dynamics ( 3 ) for all n one-element differing saturated initial states. The measure Rnij s for all possible pairs take only two values 1, 1/3 (i.e., Nij only takes value 6 or 2). The distribution of Rnij (i.e., how may pairs have the same Rnij ) is given in Table 2. For Si , the seven Sj(0)s in the 22 saturated initial states with Rnij = 1 are Si itself and those differing from Si by only one element; the other fifteen Sj(0)s in the 22 saturated initial states with Rnij = 1/3 are those differing from Si by two elements. The Sj(0) differing from Si by more than two elements does not converge to the Si. 4 Robustness to mutations For robustness to noise, we need to find all possible saturated equilibrium expression states for a given gene network W. In contrast, for robustness to mutations, we have the opposite task: for a given saturated equilibrium expression state we need to determine which and how many W s share this saturated equilibrium expression state. 4.1 Conditions under which gene networks share the same saturated equilibrium expression state Theorem 6. For a given saturated expression state S, all possible networks specified by particular W s having it as a saturated equilibrium expression state can be completely constructed by solving the following inequalities: j∈J+(S) j∈J+(S) wij − wij − Proof. Note that the rows of W are independent. When S is the saturated equilibrium expression state of W, the elements of each row (for example, the ith row wij(j = 1, 2, . . . , n)) of W must satisfy either ( 13 ) or ( 14 ), i.e., (29) or (30) which is an inequality with n variables. The Table 2 Distribution of Rnij for the network W in Example 2 Final state (1 1 1 -1 -1 -1) (-1 -1 -1 1 1 1) 1 7 inequality is solvable and has an infinite number of solutions. Those are the desired solutions with each wij located within the required range [ −a, a]. For discrete wij (only taking values [ −1, 0, 1]), the number of solutions is finite. All the solutions can be completely counted and determined by solving (29) or (30). From (29) or (30), we may draw the following conclusions: 1. For i ∈ J+(S), increasing the value of the first term of (29) (if the increase does not make the value of wij larger than the upper bound a) will not violate the inequality, i.e., increasing the value of wij(j ∈ J+(S)) will keep the same saturated equilibrium state S. This behavior implies that either increasing the activation or decreasing the repression influence of active gene j on active gene i at equilibrium state S will not change the saturated equilibrium state. 2. For i ∈ J+(S), decreasing the value of the second term of (29) (if the decrease does not make the value of wij smaller than the lower bound −a) will not violate the inequality, i.e., decreasing the value of wij(j ∈ J−(S)) will keep the same saturated equilibrium state S. This behavior implies that either decreasing the activation or increasing the repression influence of inactive gene j on active gene i at equilibrium state S will not change the saturated equilibrium state. 3. Similarly, for i ∈ J−(S), either increasing the activation or decreasing the repression influence of inactive gene j on inactive gene i at equilibrium state S will not change the saturated equilibrium state. 4. For i ∈ J−(S), either decreasing the activation or increasing the repression influence of active gene j on inactive gene i at equilibrium state S will not change the saturated equilibrium state. Example 6. Consider the determination of all possible gene networks W with the given saturated equilibrium expression state S = ( 1 1 1 −1 −1 −1 ). Here, the discrete values [ −1, 0, 1] are required for the elements wij. Thus, we seek to determine all gene networks W sharing the same two saturated equilibrium expression states given in Example 2 (due to Theorem 2, −S is also a saturated equilibrium expression state). First, consider i ∈ J+(S). Set β = 5 for the 0.99 confidence level. In this case, (29) becomes wij − wij ≥ 5, (i = 1, 2, 3). (31) 1/3 15 wij = −3 j=4 In this case, (32) is 3 j=1 wij ≥ 2. Note that the second term on the right-hand side of (32) takes on integer values from −3 to 3 corresponding to all three wij having the value either −1 or 1, respectively. We treat each circumstance separately: It is easy to see that there is only one choice for (wi4 wi5 wi6 ) = (−1 −1 −1) and four choices for (wi1 wi2 wi3 ) = (1 1 1), (0 1 1), (1 0 1), (1 1 0). Thus, we have four choices for wij(j = 1, 2, . . . , 6) (36) (38) (39) This criterion is impossible because j3=1 wij cannot be larger than 3. Therefore, altogether, there are only seven choices or permitted rows for wij when i ∈ J+(S). wij ≥ −1 j=4 In this case, (32) is 3 wij = −2 j=4 In this case, (32) is 3 j=1 wij ≥ 3. It is easy to see that there is only one choice for (wi1 wi2 wi3 ) = (1 1 1) and three choices for (wi4 wi5 wi6 ) = (0 −1 −1), (−1 0 −1), (−1 −1 0). Thus, we have another three choices for wij(j = 1, 2, . . . , 6) w5+ = (1 1 1 w6+ = (1 1 1 −1 0 −1 −1), Note that (38) is the same as (32) except that (wi1wi2 wi3) and (wi4wi5wi6) interchange their positions. Therefore, there are seven choices or permitted rows for wij when i ∈ J−(S): In the construction of wk+, for a given (wi4, wi5, wi6), all possible patterns of (wi1, wi2, wi3) are considered. For a given (wi1, wi2, wi3), a similar treatment for (wi4, wi5, wi6) was performed. Thus, for any wk+ , we can always find a wl+ differing from it only by a single wij. This is also true for w−. In this example, wk+ and wk−(k = 2, 3, . . . , 7) differ k from w1+ and w1− only by a single wij, respectively. Since each row of W has seven choices, altogether, there are 76 = 117, 649 gene networks, each specified by a particular W, sharing the same saturated equilibrium expression states These gene networks sharing the same saturated equilibrium expression states are referred to as ‘viable’ networks (see [20]). This definition is different from that given by Ciliberti et al. [13], where viable networks were those sharing a prespecified pair of saturated initial and equilibrium expression states. The viable networks defined by sharing a pair of saturated initial and equilibrium expression states are a subset of the viable networks defined by sharing a saturated equilibrium expression state only. The analysis here about viable networks implies that for a given saturated state, all viable networks having it as an equilibrium state must follow certain patterns, i.e., its rows must be chosen from finite permitted rows. The permitted rows for a given saturated equilibrium state have specific biological meaning and reflect the required connectivity patterns of each gene to other genes. This restriction distinguishes viable networks for a given equilibrium state from other viable networks with distinct equilibrium states as well as inviable networks. 4.2 Definitions of robustness to mutation The number of viable networks in the example above 76 = 117, 649 itself is large, but the total number of possible gene networks with n = 6 is 362 ≈ 1.5 × 1017. The fraction of viable gene networks in the total number of possible gene networks is 76/336 ≈ 7.8383 × 10−13, even smaller than that obtained previously in numerical simulations [13] for n = 20. Based on the above analysis for viable gene networks, it seems plausible to define robustness to mutation for a given saturated equilibrium state as the ratio of the number of viable networks to the total number of possible networks for a given n. If so, it would appear that none of the viable gene networks is robust to mutation because the ratio is very small. The latter inference is misleading because most of the possible gene networks have no similarity in topology with the viable gene networks having a specific saturated equilibrium state, and there is rarely a chance that a viable gene network will suddenly change to one of them. In normal circumstances, the structure perturbations due to mutation are small and the topology can only change gradually. A viable network may experience topology changes stepby-step, and in each step, only one wij changes. Ciliberti et al. [13] defined mutational robustness for a viable gene network as the fraction of its one-mutant neighbors that are also viable, and we follow the same criterion. In the following discussion, wij is restricted only to take the discrete values [ −1, 0, 1]. We will use Example 5 as an illustration. A gene network W is viable if and only if its ith row belongs to one of the seven rows in W + (if i ∈ J+(S)) or W − (if i ∈ J−(S)), respectively. Since wij only takes on three values [ −1, 0, 1], each wij may have two possible changes from its original value, and there is a total 2n2 = 2 × 62 = 72 single wij changes (i.e., 2n2 = 72 one-mutant neighbors) for any W. Rmi = ⎧ 0, ⎪⎪⎨ if W is inviable, ⎪⎪ NNWWv = i=n1 N2nwv2i , if W is viable. ⎩ Here the subscripts m and i in Rmi denote ‘mutation’ and the ith saturated equilibrium expression state; N Wv and Nwvi respectively are the total numbers of viable single wij changes (which is also the number of one-mutant viable neighbors) of W and its ith row; NW is the total number of possible single wij changes of W. For a given W with a specified saturated equilibrium state, its robustness to mutation can be readily calculated from Nwvi of each row. In Example 5, since each row of a viable network must be one of the seven rows in W + or W −, and Nwv+ = Nwv− = 6, 1 1 Nwv+ = Nwv− = 1(k = 2, . . . , 7), the value of Rmi depends k k on how many w1+ and w1− are contained in W. Example 7. Some inviable and viable networks Wk(k = 1, . . . , 8) with respect to the saturated equilibrium expression state (1 1 1 −1 −1 −1) (and (−1 −1 −1 1 1 1) by Suppose that only one wij can change. From W +, we see that each element w1+j in w1+ changing to 0 yields one of wk+(k = 2, ..., 7), which is still viable. Other changes are not viable. Therefore, the total number of viable single w1+j changes in w1+ is However, for each wk+j in wk+(k = 2, . . . , 7) only 0 changing to 1 (if j ∈ J+(S)) or −1 (if j ∈ J−(S)) gives w1+ which is still viable, and other changes are not viable. Thus, the total number of viable single wk+j changes in wk+(k = 2, . . . , 7) is v Nw+ = 1, k (k = 2, . . . , 7). It can be proved that this is also true for w1− and wk−(k = 2, . . . , 7), i.e., v Nw− = 6, 1 v Nw− = 1, k (k = 2, . . . , 7). (42) We then define robustness to mutation as follows: Definition 5. Robustness to mutation for a viable gene network W with a specified saturated equilibrium state Si is (40) (41) (43) Theorem 2) are given below and their N Wv and Rmi are given in Table 3. W1 is inviable because its first row does not belong to any row of W + or W −. The other Wk s are viable. As mentioned above, the value of Rmi of a viable W depends on how many w1+ and w1− are contained in the network. From W2 to W8, more and more w1+ and w1− are included. Thus, the corresponding N Wv and Rmi become larger and larger, and W8 is the most stable one with N Wv = 36 (i.e., having 36 one-mutant viable neighbors) and Rmi = 1/2. For n = 6 and the saturated equilibrium expression states (1 1 1 −1 −1 −1) and (−1 −1 −1 1 1 1), the viable gene networks can only take seven distinct values of Rmi given in Table 3 corresponding to W containing 0,1,. . . ,6 of w1+ and w−. Suppose k rows of a viable network are w+ 1 1 and w1−, then the value of its mutational robustness Rmi (k) is given by 6 Rmi (k) = Vwvi /(2 × 62) =[ 6k + 1(6 − k)] /72 i=1 = (5k + 6)/72, (k = 0, 1, . . . , 6). It can be readily proved that the number N (k) of all possible viable networks with robustness to mutation equal to Rmi (k) is N (k) = C6 × 66−k, (k = 0, 1, . . . , 6), k Table 3 The robustness to mutation Rmi of W k W k W1 W2 W3 Rmi 0 6/72 11/72 where the first term C6k denotes the number of combinations of k elements taken from six elements, representing how many possible positions that k rows of w1+ or w1− can take in six rows of W. Each of the remaining 6 − k rows of W has six choices from wk+ and wk−(k = 2, . . . , 7), and the second term 66−k gives the total possible number of combinations for the 6 − k rows. Figure 6 gives the distribution of Rmi . Following [13], we also define robustness to mutation, Rmij , as follows: Definition 6. Robustness to mutation for a viable gene network W with specified saturated equilibrium and initial states Si and Sj(0) is (46) Rmij = N mvij NW N mvij . Here, the subscript m and i, j in Rmij respectively denote ‘mutation’ and the ith saturated equilibrium expression state Si and jth saturated initial expression state Sj(0); N mvij is the total numbers of viable single wij changes of W with respect to Si and Sj(0), which can be obtained by testing how many viable networks in N Wv share Sj(0). Example 8. Determination of Rmij for W2 and W3 in Example 7. Rmij values were determined for the networks W2 and W3 given in Example 7. Both W2 and W3 have the same saturated equilibrium states S1, S2. For W2 each Si has 22 saturated initial states converging to it. For W3 each Si has 32 saturated initial states converging to it. The numbers N Wv of one-mutant viable networks sharing the same saturated equilibrium state Si(i = 1, 2) for W2 and W3 are 6 and 11, respectively (see Table 3). For W2, all 6 one-mutant neighbours sharing Si also share all 22 Sj(0)(j = 1 − 22), i.e., N mvij = N Wv = 6, Rmij = 6/72 for all 22 pairs. For W3, N mvij is different not only for distinct Si but also for distinct initial states. For S1, N mvij (≤ N Wv ) takes values: 5,7,11; for S2, N mvij (≤ N Wv ) takes values: 5,6,9,11, respectively. Table 4 gives the distribution of Rmij values for W3, i.e., how many pairs of Si and Sj(0) take these values. 4.3 Topology evolution of gene networks From the procedure to construct viable networks given in Example 6, we know that for each permitted row there always exists another permitted row differing by only a single wij from it. Therefore, for a viable network Wi, we can always find one or more viable networks, Wj’s, differing by only one wij from it. The above eight networks Wk(k = 1, 2, . . . , 8) are an example. They only differ from one another as adjacent neighbors with a single changed wij. These changes in topology correspond to the loss of a 4.5 4 ) (k 3.5 N , ks 3 r o tew 2.5 n fro 2 e bm1.5 u N 1 0.5 regulatory interaction (wij → 0), or to the appearance of a new regulatory interaction that was previously absent. The changes can be represented as a reversible path W1 ⇔ W2 ⇔ W3 ⇔ W4 W8 ⇔ W7 ⇔ W6 ⇔ W5 In going from W2 to W1, the gene network no longer attains the saturated equilibrium expression state. Thus, we may consider W1 as ‘dead’. In going from W2 to W3, however, not only is the saturated equilibrium expression state retained, but also the robustness to mutation becomes higher. Suppose that all possible single wij changes have the same probability, then the gene network with higher Rmi has a greater chance to ‘survive’. This implies that highly robust topologies can evolve from topologies with low robustness through gradual Darwinian topological changes or ‘natural selection’. Ciliberti et al. [13] suggested that all viable networks attaining a given gene expression state can be organized into a graph whose nodes are networks that differ in their Table 4 The distribution of Rmij for the network W 3 Final state S1 6/72 0 topology. Two networks (nodes) in the graph are connected by an edge if they differ in the value of only one regulatory interaction (wij). As proved above, for a viable network Wi, we can always find one or more viable networks, Wjs differing by only one wij from it. The number of viable neighbors differing by a single wij for a viable network W is simply the value of its N Wv (see Table 2). Therefore, any two viable networks Wi and Wj with k different elements wij can be connected by a path with k edges and k − 1 viable networks between them. For example, W2 and W8 have six different diagonal wijs, and they are connected by a path with six edges and five viable networks between them. This circumstance implies that all viable networks can be organized to comprise a large graph which can be easily traversed by a sequence of single wij changes of network topology. Thus, robustness is an evolvable property. To draw this conclusion, a previous study performed millions of simulations [13], but the analytical treatment here directly leads to this result. 5 Modified sigmoidal function with threshold parameters All the results obtained here are based on the sigmoidal function model ( 3 ) for gene networks. This model is a simplified picture, and caution is called for so as to not over-interpret the conclusions obtained from our analytical treatment. For example, we proved that −S is also a saturated equilibrium expression state if S is one; this conclusion may not be biologically meaningful. Another conclusion from our analysis is that 0 may be an unstable fixed point. Following a previous definition, 0 corresponds to all genes being ‘half-on’. This definition may not be appropriate under some circumstances, and instability of 0 introduces difficulty for biological interpretation. However, these considerations provide criteria to modify the mathematical model, for example, by using the more general sigmoidal function proposed in [23] to describe network dynamics. To remove −S and 0 from being an equilibrium state or a fixed point, the complex sigmoidal function given in [23] is unnecessary, we only need to slightly modify the sigmoidal model ( 3 ) by introducing a threshold parameter θi [21,22]: i.e., (50) may not be satisfied, and −S may not be a saturated equilibrium state. If θi < 0 for all i, then −β − θi > −β + θi, and it is possible that n j=1 wij(−Sj) ≤ −β + θi, if i ∈ J−(−S). In this case, (51) may not be satisfied, and −S may not be a saturated equilibrium state. If θi < 0 for all i ∈ J+(S) and θi > 0 for all i ∈ J−(S), then both (50) and (51) may not be satisfied for −S, and −S may not be a saturated equilibrium state for W. In Example 2, when the model (48) is used with θi = −2 (i = 1, 2, 3) and θi = 2(i = 4, 5, 6), then the network given in ( 17 ) only has a single saturated equilibrium state (1 1 1 −1 −1 −1), and all saturated initial states converge to it. Thus, the problem reduces to choosing the parameter θi and giving it biological interpretation. Then, we have Theorem 7. The necessary and sufficient condition for a saturated state S to be an equilibrium expression state or a fixed point attractor of the dynamics (48) with a given matrix W is wij − wij − j∈J−(S) j∈J−(S) ⎞ (56) (57) − 1, − 1 (48) (49) (50) (51) (54) (55) i.e., S is a saturated fixed point with 0.99 (β = 5) or 0.9999 (β = 10) confidence level. A saturated state S satisfying (58,59, or 60) is a fixed point attractor. The proof is given in Theorem A5 in the Appendix. Similarly, we can have Theorem 8. For a given saturated expression state S, all possible networks specified by particular W s having it as a saturated equilibrium expression state for dynamics (48) can be completely constructed by solving the following system of inequalities wij − wij − under the condition wij ∈[ −a, a], and θi ∈[ −b, b]. Proof. The proof is the same as that for Theorem 6. 6 Application to a yeast cell-cycle network A simple yeast cell-cycle network shown in Figure 7b with 11 nodes was proposed by Li et al. [24]. The dynamics of the network was defined by Li et al. as aijSj(t) > 0, aijSj(t) < 0, aijSj(t) = 0, illustrate our analytical treatment for network construction and its robustness analysis. Hereafter, 0 representing the inactive state by Li et al. will be replaced by −1. Only the main results are presented here; see the online Supplementary information (Additional file 1) for more details. 6.1 Construction of viable networks Define the node order from 1 to 11 as specified in Table 5, 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 5 S1 = (−1 −1 −1 −1 1 −1 −1 −1 1 −1 −1) (65) (66) (64) for dynamics ( 3 ). According to Theorem 2, these viable networks will also share the other saturated equilibrium expression state where 1 and 0 correspond to active and inactive states of the gene, i.e., 0 instead of −1 is used to represent the inactive state, and aij is wij in ( 3 ). Using model (64), Li et al. found that there exist 7 saturated fixed point attractors (considering 0 as −1) and all of the 211 = 2, 048 possible saturated initial expression states converge to one of the seven fixed point attractors (see Table 5). Note that dynamics (64) is different from dynamics ( 3 ). For j aijSj(t) = 0, dynamics ( 3 ) gives Si(t + 1) = 0, not Si(t). Dynamics ( 3 ) with the W constructed directly from the connectivities in Figure 7b will not give the same result as that given by dynamics (64). All the information given by the simplified model for yeast cell-cycle network (Figure 7b) will be considered as ‘available experimental information’ for budding yeast, and used as an example to As mentioned by Li et al. [24], ‘the overall dynamic properties of the network are not very sensitive to the choice of these parameters’ (wij), but the connectivity patterns of the network, i.e., the regulatory influence between genes (activation, repression, and absence) is important for determining gene network robustness. Therefore, we restrict wij to only take the discrete values 1 (activation), −1 (repression), and 0 (absence). Theorem 6 gives the criterion to construct all of such networks W. When wij only takes values [ −1, 0, 1], to satisfy (29,30) each row of W must have five or more nonzero elements due to β ≥ 5. Otherwise, the network would not have any saturated equilibrium states. This problem occurs not only for networks with less than five genes, but also for larger networks with sparse connectivities 0 0 0 0 0 0 0 between genes. For example, Node 1 (Cln3) in Figure 7b is a pure ‘parent’ node, which does not have any regulation coming from all other ‘children’ nodes, i.e., all w1j = 0 for j = 1, and for S1 the condition (30) does not hold: w1j − j∈J+(S1) j∈J−(S1) w1j = −w11 ≤ −β. W = βWˆ , To avoid this problem, the factor β may be introduced such that so to satisfy condition (29,30), wˆij can only take values [ −1, 0, 1] without any restriction on the number of nonzero elements in each row of Wˆ . For the sake of notational simplicity, in the sequel, we still use W instead of Wˆ , but write dynamics ( 3 ) as Si(t + τ ) = σ ⎝ wijSj(t)⎠ = ⎛ n j=1 ⎞ 2 n j=1 1 + exp −β wij − wij − For saturated equilibrium state S1, (70,71) may be rewritten as − − 11 wij ≥ 1 − wi5 − wi9, if i = 5, 9, wij ≤ −1 − wi5 − wi9, if i = 5, 9, 1 0 1 0 0 0 1 (67) (68) (70) (71) (72) (73) 0 0 0 0 0 0 0 or 11 Using the condition (74,75), all permitted row patterns sharing saturated equilibrium state S1 for dynamics (69) have been completely counted and determined (see Additional file 1: Supplementary information). Each row has 72,219 permitted patterns. Thus, the total number of viable networks sharing saturated equilibrium state S1 is As shown below, for the yeast cell-cycle network, the first row of W is restricted to be then the total number of viable networks for dynamics (69) is There are many choices of practically relevant networks. Similarly, the dynamics with threshold parameters (48) is also modified as ⎛ n ⎞ Si(t + τ ) = σ ⎝ j=1 wijSj(t)⎠ 2 n j=1 1 + exp −β( wijSj(t) − θi) − 1. (76) The necessary and sufficient condition to have S as a saturated equilibrium state for dynamics (76) becomes and for S1 (77,78) become and Condition (79,80) can be used to construct all viable networks W for dynamics (76) with a given set of θis sharing the saturated equilibrium state S1. Since there is no unambiguous biological interpretation for the values of θi, as [ −1, 0, 1], to represent activation, repression, absence for wij, we will not construct all such viable networks here. 6.2 Construction of yeast cell-cycle networks According to the definition for the green and red arrows along with the yellow loop in [24], the network directly constructed from the connectivities of Figure 7b is ⎡ −1 0 0 0 0 0 0 0 0 0 0 ⎤ 1 0 0 0 0 0 0 0 0 −1 0 ⎥ 1 0 0 0 0 0 0 0 0 −1 0 ⎥⎥ 0 0 1 −1 0 0 0 0 0 0 0 ⎥⎥ 0 0 0 −1 0 0 1 −1 0 −1 0 ⎥⎥ 0 0 0 0 0 −1 1 0 0 −1 1 ⎥⎥ . 0000 0010 0000 −0010 −1000 1000 −−−1111 −0110 −−1010 −0011 1001 ⎥⎦⎥⎥⎥⎥⎥⎥ 0 0 0 0 0 0 0 1 0 1 −1 (81) W0 does not satisfy condition (70,71) for any saturated state and does not have a saturated equilibrium state for dynamics (69). However, W0 will be used as the basis for the connectivities of the network to construct networks for dynamics (69) with the saturated equilibrium expression state S1. The construction of networks reduces to satisfying condition (74, 75) or (79, 80) as much as possible consistent with experimental observation. Two yeast cell-cycle networks for dynamics (69) have been obtained (the detailed procedure for their construction can be found in the Additional file 1: Supplementary information). W1 and W2 differ only for w8,10. We can also use W0 without any change, but introduce the threshold parameters for dynamics (76) satisfying condition (79, 80). One choice with the smallest magnitudes of θis 6.3 Saturated equilibrium expression states for constructed networks The saturated equilibrium expression states for a given network W in dynamics (69) can be determined by using the modified condition of ( 15 ) in Theorem 1 11 ⎛ n j=1 ⎞ For n = 11, there are 211 = 2, 048 saturated states. All of the 2,048 states were tested by condition (88) for W1 and W2, respectively, to determine which of them are saturated equilibrium states for W1 and W2. The test for 2,048 states took only 0.01 s by Matlab on a Dell Precision Workstation T3400. The saturated equilibrium expression states for a given network W with threshold vector in dynamics (76) can be determined by using the condition ⎛ n j=1 ⎞ W 2 979 979 45 45 Table 6 The number of saturated initial states converging to different equilibrium states for different gene networks The saturated equilibrium expression states for W1, W2, and W0 with are shown below. 6.4 Robustness to noise First, the number of saturated initial expression states converging to each equilibrium expression state for W1, W2 is determined by either directly solving the dynamics (69) or using modified condition of Theorem 5 (90) (91) the saturated equilibrium states. The resultant robustness to noise measures Rnt and Rni are given in Table 7. There are significant differences between Rni (i = 1, 2, 3, 4) for W2. Obviously, the saturated equilibrium states S1, S2 are much more stable than S3, S4. The robustness to noise Rnij for each pair of saturated equilibrium and initial expression states was calculated. The distribution of Rnij , i.e., how many pairs with the same value of Rnij , is given in Tables 8 and 9. The results show that W0 with is completely stable for any viable pair; for W1, there is one neighbour of Sj(0) differing at the first element, which causes changes in the saturated equilibrium state Si; for W2, the distribution of Rnij is divergent, and S1 and S2 are much more stable than S3 and S4. 6.5 Robustness to mutation The Rmi values have been calculated for W1, W2, and W0 with the given in (85) as shown in Table 10. Note that for S1, Rmi is almost the same for W1, W2, and W0 with . Robustness to mutation Rmij for a viable pair of specified saturated equilibrium and initial states Si and Sj(0) has also been calculated for W1, W2, and W0 with . The resultant distribution, i.e., how many pairs having the same Rmij for W1, W2, and W0 with the given above is shown in Figure 8 and Table S12 in Additional file 1: Supplementary information. 7 Conclusion Based on the determination of saturated fixed point attractors for the sigmoidal function model in ( 3 ) with a given gene network, W, one can analytically determine Table 7 The Robustness to noise Rnt and Rni for different gene networks Rnt 1 1/2 For W1, the CPU times are 0.8 and 0.3 s, respectively to check all 2,048 saturated states, i.e., using Theorem 5 the CPU time is approximately 41% of that for the direct solving of the sigmoidal function. The results are given in Table 6. Note that for W1, W2, no saturated initial state converges to the unstable fixed point 0. Therefore, in the calculation of Rnt , we ignore 0 and only consider Network W0 with W1 Table 8 The distribution of Rnij for the networks W 0 with and W 1 Table 10 The robustness to mutation Rmi of W 0 with , W 1, and W 2 model construction and analysis of other properties for gene networks. Appendix The appendix proves several theorems in the main text. Lemma 1 (Banach Fixed Point Theorem [25]). Let (X, d) be a non-empty complete metric space with a contraction mapping g : X → X. Then g admits a unique fixed point x∗ in X (i.e. g(x∗) = x∗). Furthermore, x∗ can be found as follows: start with an arbitrary element x0 ∈ X and define a sequence {xn} by xn = g(xn−1), then xn → x∗. A map g : X → X is called a contraction mapping on X if there exists q ∈[ 0, 1) such that d(g(x), g(y)) ≤ q d(x,y) where d denotes the distance, for all x, y ∈ X. A continuous function g satisfies the Lipschitz condition g(x) − g(y) p ≤ sup J(t) p x − y p t∈X where J(x) is the Jacobian of g(x). Its (i, j)th entry is which and how many saturated equilibrium expression states exist. Furthermore, for each saturated equilibrium expression state of a W, which and how many saturated initial expression states converging to it can also be determined. These results make it possible to establish the robustness of a given gene network to noise without performing a large number of simulations. Based on the necessary and sufficient condition for gene networks to share the same saturated equilibrium expression state, one can determine all the viable gene networks for a specified saturated equilibrium state. This result also makes it possible to establish the robustness to mutation for a network with a specified saturated equilibrium expression state or a specified pair of saturated equilibrium and initial expression states. The analytical treatment presented here proved that for a given saturated state, all viable gene networks having it as an equilibrium state must follow certain patterns, i.e., the rows of the corresponding W must be chosen from a finite number of permitted rows. The permitted rows for a given saturated equilibrium state have specific biological meaning and reflect the required connectivity patterns of each gene to other genes. This restriction distinguishes the viable networks for a given saturated equilibrium state from other viable networks with distinct saturated equilibrium states as well as inviable networks. The analysis also proved, without performing a very large numbers of simulations, that all viable networks can be organized as a large graph which can be easily traversed by a sequence of single wij changes of network topology. Thus, robustness is an evolvable property. Highly robust topologies can evolve from topologies with low robustness through gradual Darwinian topological changes or natural selection. The analytical treatment presented in this paper may be employed not only for robustness analysis but also for the Table 9 The distribution of Rnij for the network W 2 Rnij 2/11 3/11 4/11 5/11 6/11 7/11 8/11 9/11 10/11 1 7 1 1 3 7 29 29 146 146 791 791 NvW 143 143 131 131 Rmsi 0.59 From (95 to 97) we see that J 1 is the largest sum of the absolute values of the elements in each column; J ∞ is the largest sum of the absolute values of elements in each row; and J 2 is the square root of the largest eigenvalue for matrix JT (t)J(t). Now define d of g(x) and g(y) as the Lp-norm of their difference. If g is a contraction W0 with Θ, S1 W1 and S1 (S2) 140 120 rs100 i a p lbe 80 a i v f ro 60 e b m uN40 20 8 7 6 s r i pa 5 e l b a ifv 4 o r be 3 m u N 2 mapping, (92) requires that at least one of the Lp-norms of its Jacobian satisfies J p ≤ q < 1. (98) This condition is sufficient, but not necessary. It is possible that one of (95 to 97) is satisfied, but the other two may not. Such examples can be constructed. Theorem A1. The sufficient condition for the sigmoidal function in ( 3 ) to have a unique fixed point attractor 0 is |wij| ≤ 2, (i = 1, 2, . . . , n). (99) ui(k) = wijSj(k). Proof. Since 0 is a fixed point for any W, to see whether it is a unique fixed point attractor, we need to determine under what condition ( 3 ) is a contraction mapping. The (i, j)th entry of the Jacobian J(k+τ )(S(k)) for ( 3 ) with τ = 1 is J(k+1)(S(k)) = ∂Si(k + 1) ij ∂Sj(k) 2e−ui(k) = 1 + e−ui(k) 2 wij 2e−ui(k) = 1 + 2e−ui(k) + e−2ui(k) wij 2 = eui(k) + 2 + e−ui(k) wij where n (100) (101) The matrix form of the Jacobian J(k+1)(S(k)) is J(k+1)(S(k)) = V (k)W where V (k) is a diagonal matrix with 2 Vii(k) = eui(k) + 2 + e−ui(k) . The ∞-norm is J(k+1)(S(k)) ∞ = 2 max S∈[−i1,1]n eui(k) + 2 + e−ui(k) j=1 (i = 1, 2, · · · , n). | wij |, (104) n Because eui(k), e−ui(k) > 0, the maximum of 2/(eui(k) + 2 + e−ui(k)) is given by the minimum of eui(k) + e−ui(k) which can be obtained from d eui(k) + e−ui(k) dui(k) This is true if and only if eui(k) = e−ui(k) = eui(k) − e−ui(k) = 0. which yields ui(k) = 0, eui(k) = e−ui(k) = 1. The minimum for eui(k) + e−ui(k) is 2, and then 2 max S∈[−i1,1]n eui(k) + 2 + e−ui(k) = 2 1 . Figure 9 gives the comparison of ex + e−x and 2/(ex + 2 + e−x). x − e+ 30 x e 60 50 40 20 10 0−4 (102) (103) (105) (106) (107) i.e., j=1 n j=1 0.5 0.45 Therefore, we have Notice that 1/2 is the superior value in (107). In many cases, not all Si = 0 (i.e., ui(k)’s are not zero), then the factor in (107) has values smaller than 1/2, which implies that the condition for a contraction mapping given in (110) may be softened as The condition in (111) is sufficient, but not necessary. Theorem A2. A saturated state S satisfying ( 13 ) and ( 14 ) is a fixed point attractor. Proof. Suppose S is a fixed point of ( 3 ). Let B = {Sˆ ∈ n | X = Sˆ − S < }, where is chosen sufficiently small such that there is only a single fixed point S within B . X = 0 is a fixed point in B with representation X because S is a fixed point for Sˆ in B . If we can prove that X = 0 is a fixed point attractor in B , then so is S. Subtracting Si on the both sides of ( 3 ) yields Xi(t + τ ) = Si(t + τ ) − 1 = j=1 where β ≥ β. The (i, j)th entry of the Jacobian is = = = = = n j=1 2 n j=1 n j=1 n j=1 2 2 n j=1 2 J(t+τ)(X(t)) = ij 2e−β exp − j=1 n j=1 (112) 2e−β exp − n j=1 wijXj(t) ≈ 1 were used. The L1-norm of Here, the condition X(t) < row vector Ji(t+τ) is 2 Ji(t+τ)(X(t)) 1 ≈ eβ + 2 + e−β n j=1 |wij|. When n j=1 |wij| < ≥ eβ + 2 + e−β eβ + 2 + e−β 2 2 > 75, if β = 5, 11, 014, if β = 10, which is often the case in practice for W satisfying ( 13 ) and ( 14 ), we have Ji(t+τ)(X(t)) 1 < 1, i ∈ J+(S). (113) (114) (115) (116) The same result can be obtained for i ∈ J−(S). Therefore, we have J(t+τ )(X(t)) ∞ = max Ji(t+τ )(X(t)) 1 < 1, i i.e., X = 0 is a fixed point attractor, and thus so is S. respectively. It can be proved that S˜i(t + τ ) = −Si(t + τ ). The proof is given below. − 1, (i = 1, 2, . . . , n), − 1, (i = 1, 2, . . . , n), Theorem A3. If an S(t) converges to a saturated equilibrium expression state S, then −S(t) converges to the saturated equilibrium expression state −S. Proof. For S(t) and −S(t) , the corresponding equations are 2 j=1 (117) (118) (119) (120) 1 + exp wijSj(t) wijSj(t) = −Si(t + τ ). (121) 1 + exp wijSj(t) − 2 exp 2 n j=1 = n j=1 n j=1 S˜i(t + τ ) = 1 − exp 1 + exp = = = 1 − = 1 − wijSj(t) wijSj(t) n n j=1 j=1 n j=1 n j=1 2 n j=1 1 + exp wijSj(t) − 2 exp n j=1 n j=1 1 + exp wijSj(t) 2 n j=1 n j=1 wijSj(t) wijSj(t) n j=1 − 1 = − 1 1 + exp − wij(−Sj(t)) 1 + exp n j=1 wijSj(t) n j=1 n j=1 Equation (121) implies that there is an one-to-one relation between S˜i(t) and −Si(t). Since S(∞) = S, we have where S˜(∞) = −S. (122) Therefore, starting from −S(t) will converge to −S. Theorem A4. A saturated initial expression state S(0) converges to a saturated equilibrium expression state S if the following condition is satisfied after a finite number k of iteration steps of ( 3 ) starting from S(0) wijSj(t ≥ k) > − ln[ (αi − 1) − wijSj(t ≥ k) < ln[ (αi − 1) − where αi = under the constraint that αi ≥ 2. Proof. First, consider i ∈ J+(S). If S is the fixed point attractor for the initial state S(0), then X = 0 is the fixed point attractor for the initial expression state X(0) = S(0) − S of (112). This is equivalent to finding the condition under which (112) is a contraction mapping. According to the Banach fixed point theorem, the sufficient condition is that the norm of the Jacobian matrix satisfies J ∞ < 1. To prove Theorem A4, we try to seek the largest Ba = {X ∈ n | X < a} where a is the upper bound to have J ∞ < 1, i.e., Ji(t+τ )(X(t)) 1 < 1, for all i ∈ J+(S). Set ⎡ y = e−β exp ⎣ − n j=1 ⎤ n j=1 |wij| 2αiy = (1 + y)2 < 1, i ∈ J+(S), which gives the quadratic equation y2 + 2(1 − αi)y + 1 = (y − y1)(y − y2) > 0, (126) (127) (128) (123) (124) or (125) j=1 n j=1 n j=1 n j=1 n j=1 It is easy to check that y1, y2 are all nonnegative when αi ≥ 2, and y2 < 1. Equation (128) is valid if and only if y is chosen within the two disjoint ranges (−∞, y2) (y1, ∞) if αi > 2, (−∞, 1) (1, ∞) if αi = 2, i.e., either smaller than y2 or larger than y1. This implies that ⎡ e−β exp ⎣ − j=1 ⎤ ⎡ wijXj(t)⎦ = exp ⎣ − wijSj(t)⎦ n j=1 ⎤ > y1 < y2 Therefore, we obtain wijSj(t) > − ln[ (αi − 1) − (134) If S(t = k) satisfies (134), then the trajectory starting from S(t = k) will converge to S. Therefore, S(t ≥ k) will also satisfy (134). Thus, wijSj(t ≥ k) > − ln[ (αi − 1) − (131) (132) (133) (135) (136) Equation (135) proves the condition for i ∈ J+(S) given in Theorem A4. For i ∈ J−(S), (112) becomes Xi(t + τ ) = Si(t + τ ) − (−1) 2 = 1 + eβ exp − wijXj(t) n j=1 , i ∈ J−(S). when S(t) is close to the fixed point attractor S. Thus, we obtain wijSj(t) < − ln[ (αi − 1) + (αi − 1)2 − 1] . (139) (αi − 1) + 1 (αi − 1)2 − 1 (αi − 1) − (αi − 1)2 − 1 [(αi − 1) + (αi − 1)2 − 1] [(αi − 1)− (αi − 1)2 − 1] Note that − ln[ (αi − 1) + (αi − 1)2 − 1] (αi − 1) − (αi − 1)2 − 1 (αi − 1)2 − (αi − 1)2 + 1 and similarly we obtain Set ⎡ y = eβ exp ⎣ − j=1 Then the proof procedure is the same as above except that for i ∈ J−(S) wijSj < 0, i ∈ J−(S) wijSj(t) < ln[ (αi − 1) − (141) wijSj(t ≥ k) < ln[ (αi−1)− (αi − 1)2 − 1] , (142) i.e., the condition for i ∈ J−(S) given in Theorem A4. Theorem A5. A saturated state S satisfying (58) and (59) is a fixed point attractor. Proof. Suppose S is a fixed point of (48). Let B = {Sˆ ∈ n | X = Sˆ −S < } where is chosen sufficiently small such that there is only a single fixed point S within B . X = 0 is a fixed point in B with representation X because S is a fixed point for Sˆ in B . If we can prove that X = 0 is a fixed point attractor in B , then so is S. n j=1 = ln = ln = ln j=1 n n j=1 Subtracting Si on the both sides of (48) yields Xi(t + τ ) = Si(t + τ ) − 1 where β = wijSj − θi ≥ β. The following step of proof is exactly the same as that in Theorem A2, and will not repeat here. Additional file Additional file 1: Supplementary information. Application to the yeast cell-cycle network. Competing interests The authors declare that they have no competing interests. Acknowledgements Support for this work was provided by DOE of USA. (137) (138) (140) = = = = n j=1 j=1 j=1 n j=1 2 j=1 2 2 1. R Thornhill , A Moller , Developmental stability, disease and medicine . Biol. Rev. Camb. Philos. Soc . 72 , 497 - 548 ( 1997 ) 2. H McAdams , A Arkin , It's a noisy business! Genetic regulation at the nanomolar scale . Trends. Genetics . 15 , 65 - 69 ( 1999 ) 3. G von Dassow, E Meir , E Munro , G Odell , The segment polarity network is a robust development module . Nature . 406 , 188 - 192 ( 2000 ) 4. M Morohashi , A Winn , M Borisuk , H Bolouri , J Doyle , H Kitano , Robustness as a measure of plausibility in models of biochemical networks . J. Theor. Biol . 216 , 19 - 30 ( 2002 ) 5. A Eldar , R Dorfman , D Weiss , H Ashe , B Shilo , N Barkai , Robustness sof the BMP morphogen gradient in Drosophila embryonic patterning . Nature . 419 , 304 - 308 ( 2002 ) 6. E Ozbudak , M Thattai , I Kurtser , A Grossman , A van Oudenaarden , Regulation of noise in the expression of a single gene . Nat. Genet . 31 , 69 - 73 ( 2002 ) 7. M Elowitz , A Levine , E Siggia , P Swain P , Stochastic gene expression in a single cell . Science . 297 , 1183 - 1186 ( 2002 ) 8. F Isaacs , J Hasty , C Cantor , J Collins, Prediction of measurement of an autoregulatory genetic module . Proc. Natl. Acad. Sci. USA . 100 , 7714 - 7719 ( 2003 ) 9. W Blake , M Kaern , C Cantor , J Collins, Noise in eukaryotic gene expression . Nature . 422 , 633 - 637 ( 2003 ) 10. N Ingolia , Topology and robustness in the Drosophila segment polarity network . PLoS Biol . 2 ( 6 ), 805 - 815 ( 2004 ) 11. A Wagner , Evolution of gene networks by gene duplications: a mathematical model and its implications on genome organization . Proc Nati. Acad. Sci. USA . 91 , 4387 - 4391 ( 1994 ) 12. A Wagner , Dose evolutionary plastcity evolve? Evolution . 50 ( 3 ), 1008 - 1023 ( 1996 ) 13. S Ciliberti , O Martin , A Wagner , Robustness can evolve gradually in complex regulatory gene networks with varying topology . PLoS Comput. Biol . 3 ( 2 ), 0164 - 0173 ( 2007 ) 14. C Espinosa-Soto , A Wagner , Specialization can drive the evolution of modularity . PLoS Comput. Biol . 6 ( 3 ), 1 - 10 ( 2010 ) 15. Y Kwon , K Cho , Quantitative analysis of robustness and fragility in biological networks based on feedback dynamics . Bioinformatics . 24 ( 7 ), 987 - 994 ( 2008 ) 16. M Siegal , A Bergman , Waddington's canalization revisited: developmental stability and evolution . Proc. Nati. Acad. Sci. USA . 99 , 8 - 10532 ( 2002 ) 17. O Matin , A Wagner , Multifunctionality and robustness trade-offs in model genetic circuits . Biophys. J . 94 , 2927 - 2937 ( 2008 ) 18. S Ciliberti , O Martin , A Wagner , Innovation and robustness in complex regulatory networks . Proc. Nati. Acad. Sci. USA . 104 , 13591 - 13596 ( 2007 ) 19. J Draghi , O Matin , A Wagner , The evolutionary dynamics of evolvability in a gene network model . J. Evov. Biol . 22 , 599 - 611 ( 2009 ) 20. R Azevedo , R Lohaus , S Srinivasan , K Dang , C Burch , Sexual reproduction selects for robustness and negative epistasis in artificial gene networks . Nature . 440 , 87 - 90 ( 2006 ) 21. E Mjoisness , D Sharp , J Reinitz , A connectionist model of development . J. Theor. Biol . 152 , 429 - 453 ( 1991 ) 22. I Mendoza , E Alvarez-Buylla , Dynamics of the genetic regulatory network of Arabidopsis thaliana flower morphogenesis . J. Theor. Biol . 193 , 307 - 319 ( 1998 ) 23. T Feng , B Tirozzi , An analysis on neural dynamics with saturated sigmoidal functions . Comput. Math. Appl . 34 , 71 - 99 ( 1997 ) 24. F Li , T Long , Y Lu , Q Quyang , C Tang , The yeast cell-cycle network is robustly designed . PNAS . 101 ( 14 ), 4781 - 4786 ( 2004 ). Copyright 2004 National Academy of Sciences , USA 25. A Granas , J Dugundji , Fixed Point Theory . (Springer-Verlag, New York, 2003 )