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 )