Eur. Phys. J. B 77, 233–237 (2010)
DOI: 10.1140/epjb/e2010-00176-4
Regular Article
THE EUROPEAN
PHYSICAL JOURNAL B
Regulatory networks and connected components of the neutral
space
A look at functional islands
G. Boldhausa and K. Klemmb
Bioinformatics Group, Dept. of Computer Science, University of Leipzig, Germany
Received 14 December 2009 / Received in final form 1st April 2010
Published online 17 June 2010 – c© EDP Sciences, Societa` Italiana di Fisica, Springer-Verlag 2010
Abstract. The functioning of a living cell is largely determined by the structure of its regulatory network,
comprising non-linear interactions between regulatory genes. An important factor for the stability and
evolvability of such regulatory systems is neutrality – typically a large number of alternative network
structures give rise to the necessary dynamics. Here we study the discretized regulatory dynamics of the
yeast cell cycle [Li et al., PNAS, 2004] and the set of networks capable of reproducing it, which we call
functional. Among these, the empirical yeast wildtype network is close to optimal with respect to sparse
wiring. Under point mutations, which establish or delete single interactions, the neutral space of functional
networks is fragmented into ≈4.7×108 components. One of the smaller ones contains the wildtype network.
On average, functional networks reachable from the wildtype by mutations are sparser, have higher noise
resilience and fewer fixed point attractors as compared with networks outside of this wildtype component.
1 Introduction
Neutrality [1] is crucial for robustness and evolvability [2]
of biological systems. It describes the fact that the map-
ping from genotypes to phenotypes is not invertible. A
given phenotype can be encoded by more than one geno-
type. As Wagner [2] writes, “most problems the living have
solved have an astronomical number of equivalent solu-
tions, which can be thought of as existing in a vast neutral
space”.
Computational studies of biopolymers revealed the ex-
istence of neutrality in the relation between sequence and
spatial structure. RNA molecules and proteins are gen-
erated as a chain (sequence) of nucleic bases and amino
acids respectively. The number of sequences folding into
one and the same functionally relevant spatial structure is
found to be large. It is growing exponentially with the size
of the molecule [3,4]. Together with an adjacency given by
single mutations, the phenotypically equivalent genotypes
form the neutral network (or neutral graph). The proper-
ties of this graph, in particular its connectivity, determine
the robustness of the given genotype under mutations and
its evolvability towards new phenotypes.
Going from single molecules to the level of the whole
organism, the phenotype is not given by the set of its
 This work has been originally presented at the European
Conference on Complex Systems, held in Warwick from 21 to
25 September 2009.
a e-mail: boldhaus@bioinf.uni-leipzig.de
b e-mail: klemm@bioinf.uni-leipzig.de
molecule structures alone: the dynamics that arises as the
result of activating and suppressing interactions between
molecules is crucial. This set of interactions is captured as
a regulatory network [5] and gives rise to a temporal se-
quence of chemical concentration vectors that are respon-
sible for the division of a single cell or the development of
an embryo. Again, the mapping from genotypes (interac-
tion networks) to phenotypes (temporal sequences) is not
injective, i.e. several network structures are able to pro-
duce the regulatory dynamics of a given phenotype [6–9].
Here we apply the neutral graph concept to a dynamical
model [10] of cell cycle regulation in the organism of the
yeast species Saccharomyces cerevisiae (budding yeast).
In Section 2 we introduce the model dynamics and the
wiring of the wildtype network. The ensemble of functional
networks that yield dynamics equivalent to the wildtype
is analyzed in Section 3, finding the neutral graph to be
disconnected. In Section 4 we focus on statistical proper-
ties of the subset of networks that are reachable from the
wildtype. After a remark (Sect. 5) on the computation of
network statistics, Section 6 offers a discussion and open
questions.
2 Cell cycle network and Boolean dynamics
During the process of cell division, a eukaryotic cell grows
and divides into two daughter cells. A cell cycle consists
of four distinct and separate phases named G
1
, S, G
2
and
M . In the G
1
(“growth”) phase, the cell commits itself
for cell division under certain conditions. In particular,
a necessary cell size must have been reached. A copy of

234 The European Physical Journal B
Fig. 1. (Color online) (a) The Cell Cycle Network of the yeast wildtype has 11 nodes connected with activating (green) and
inhibiting (red) interactions. Self-suppression is indicated by yellow loops. (b) A sequence of 13 states defines a cell cycle, as
produced by the network in (a). (c) A different network (mutant) performs the same sequence of states. As the wildtype, this
mutant has 34 interactions. However, 19 entries in the interaction matrix differ from the wildtype.
the genetic information is produced in the S (“synthesis”)
phase. The G
2
(“gap”) phase precedes the actual cell di-
vision in the M phase (“mitosis”).
Here we are interested in the network of molecules
(cyclins, inhibitors and degraders of cyclins and transcrip-
tion factors) regulating this process. We consider the reg-
ulatory network of the mono cellular eukaryotic organ-
ism Saccharomyces cerevisiae (budding yeast). Its genome
comprises 13 million base pairs and 6275 genes, of which
approximately 800 are involved in the cell cycle dynam-
ics [11]. The dynamics is controlled by a core of 11 key
regulators with 34 directed interactions [10], shown in
Figure 1a, which we denote as the wildtype network. In-
teractions are captured by a matrix A. If node j has an
activating effect on node i, the corresponding matrix ele-
ment is aij = +1, while inhibition is coded as aij = −1. In
case of no direct influence from j to i, we have aij = 0. Li
et al. [10] model the regulatory dynamics with a Boolean
approach [12,13] where each node i takes state values
Si(t) ∈ {0, 1} when being inactive/active at time t. In the
time-discrete dynamics, nodes are updated synchronously,
based on their weighted input sum hi(t) =
∑
j aijSj(t).
The state at the next time step is obtained by applying
the threshold update rule
Si(t + 1) =
⎧
⎨
⎩
1, hi(t) > 0
0, hi(t) < 0
Si(t), hi(t) = 0.
(1)
From an initial condition S(1), representing the real start-
ing state of the cell cycle, the dynamics produces the
sequence of state vectors S(1), S(2), . . . , S(13), shown in
Figure 1b. The state S(13) = G
1
is a fixed point of the dy-
namics. The system remains in this state until node Cln3
is externally activated. In the real system the external ac-
tivation indicates that the cell size is sufficient for another
division.

G. Boldhaus and K. Klemm: Regulatory networks and connected components of the neutral space... 235
0 20 40 60 80 100 120
number of connections
10
0
10
5
10
10
10
15
10
20
10
25
10
30
10
35
h
is
to
g
ra
m
Fig. 2. (Color online) Histograms of the number of interactions
over functional networks. Positive (green dashed curve), neg-
ative (red dot-dashed), and total connections (black solid) of
almost all functional networks exceed the corresponding counts
in the wildtype network (arrows).
3 Functional networks and the neutral graph
Broadening our treatment of regulatory networks, we con-
sider the set of all networks with interaction matrices
over 11 nodes with entries aij ∈ {−1, 0, +1}. We call a
network functional if it produces the state transitions of
the cell cycle sequence in Figure 1b. Thus, the wildtype
network is functional. However, there are further func-
tional networks. Out of the set of all 32
11
≈ 5.4 × 1057
networks, approximately 5.11 × 1034 are functional [14].
Figure 1c shows an example of a functional network dif-
ferent from the wildtype. Figure 2 shows the statistics for
the number of interactions (arcs) present in functional net-
works. The wildtype network is sparse in comparison with
the average functional network. However, there are func-
tional networks that are even sparser than the wildtype.
These findings analogously hold when activating and in-
hibiting interactions are counted separately. Interestingly,
functional networks have generally more suppressing than
activating interactions, as is the case for the wildtype.
A structure to reflect mutations on the set of functional
networks is the neutral graph. Its nodes are the functional
networks. Functional networks A and B are adjacent (con-
nected by an edge) in the neutral graph if A is turned into
B by a single mutation. According to our definition a mu-
tation is a replacement of one entry in the interaction ma-
trix. The Hamming distance between two networks is the
number of entries in which their interaction matrices dif-
fer. In order to avoid confusion with the networks of inter-
action we employ the term neutral graph as a synonym for
the more commonly used neutral network. An important
property of a neutral graph is its connectedness. A muta-
tional walk from network A to network B is a sequence of
single point mutations that turns A into B without pass-
ing through non-functional networks. The neutral graph
is connected if such a mutational walk exists for each pair
of functional networks. We find that the neutral graph
0 1×10
26
2×10
26
3×10
26
4×10
26
5×10
26
component size
10
6
10
7
10
8
cu
m
u
la
ti
v
e
h
is
to
g
ra
m
Fig. 3. (Color online) Cumulative size distribution of con-
nected components of the neutral graph (falling curve). The
component containing the wildtype has size 5.66× 1025 (verti-
cal line).
considered here is disconnected. One cannot pass from all
functional networks to all others by sequences of muta-
tions that preserve functionality. In fact, mutual reach-
ability between functional networks is rare. The neutral
graph falls into ≈4.7 × 108 connected components with
sizes distributed between ≈6.1× 1024 and ≈4.4× 1026, as
shown in Figure 2. The component of the wildtype com-
prises around 5.66 × 1025 functional networks.
4 The wildtype component
In this section we extend the analysis of the neutral graph.
We focus on a comparison between functional networks in
the wildtype component and all functional networks. Fig-
ures 4a–4c shows how the number of (a) negative, (b) pos-
itive and (c) all interactions is distributed. All three plots
reveal a significant statistical difference between networks
in the wildtype component and the set of all functional
networks. Networks in the wildtype component are sparse
compared with the average functional network.
Geometric information of the neutral graph is provided
in Figure 4d in terms of the Hamming distance of func-
tional networks from the wildtype. Functional networks in
the wildtype component are closer to the wildtype than
the average functional network is. Still the most remote
networks in the wildtype component are found at dis-
tance 77 from the wildtype. Despite its moderate size, the
wildtype component pervades a large part of the network
space.
Shifting attention from the structural to the dynam-
ical properties of the functional networks, let us analyze
the resilience of the dynamics against perturbations. As
a measure of resilience we use the G
1
basin size [10], i.e.
the number of states from which the dynamics eventually
reaches the fixed point G
1
. Clearly, the basin contains at
least the 13 states in the cell cycle sequence. As shown
by the distributions in Figure 4e, actual G
1
basin sizes in
functional networks contain many more states. Compared
with all functional networks, basin sizes of networks in

236 The European Physical Journal B
0 20 40 60 80 100 120
number of negative arcs
0
0.05
0.1
0.15
fr
ac
ti
o
n
o
f
n
et
w
o
rk
s
0 20 40 60 80 100 120
number of positive arcs
0 20 40 60 80 100 120
total number of arcs
0
0.05
0.1
fr
ac
ti
o
n
o
f
n
et
w
o
rk
s
0 20 40 60 80 100 120
distance from wildtype
0 1000 2000
basin size of G
1
fixed point
0.000
0.001
0.002
2040 2048
0.000
0.001
0.002
(a) (b)
(c) (d)
(e)
Fig. 4. (Color online) Comparison of statistics between all
functional networks (solid green curves) and functional net-
works in the wildtype component (dashed blue curves) of the
neutral graph. (a) Histogram of negative interactions; (b) his-
togram of positive interactions and (c) histogram of total num-
ber of interactions in functional networks; (d) histogram of
Hamming distances (minimal number of mutations) from the
wildtype; (e) distribution of basin sizes of the G
1
fixed point.
The inset shows a zoom into the histogram for very large basin
sizes. For comparison, a vertical dotted line in each panel gives
the value of the wildtype network itself. Histograms in pan-
els (a–d) are exact. Histograms in (e) were obtained by uni-
form sampling of 107 functional networks each from the whole
neutral graph and from its wildtype component, respectively.
the wildtype component concentrate at higher values. The
most frequently observed basin size is 2047 for networks in
the wildtype component, cf. the inset of Figure 4e. How-
ever, we have not found a functional network where the
G1 basin contained all 2048 states.
Moreover, the distributions in the number of fixed
points of functional networks show a striking difference
between the wildtype component and the whole neutral
graph. Figure 5a displays geometric distributions in both
cases. However, networks in the wildtype component have
a significantly narrower distribution of fixed points. In-
terestingly, dynamic attractors (limit cycles) with more
than one state show practically the same statistics in the
wildtype component as in the whole neutral graph, cf.
Figure 5b.
0 50 100
number of fixed points
10
-6
10
-5
10
-4
10
-3
10
-2
10
-1
10
0
cu
m
u
la
ti
v
e
h
is
to
g
ra
m
0 5 10
number of limit cycles
(a) (b)
Fig. 5. (Color online) The number of attractors of all func-
tional networks (green curves) and the functional networks in
the neutral graph component containing the wildtype (blue
curves). (a) Distribution of the number of fixed points. (b) Dis-
tribution of the number of limit cycles (attractors of length at
least 2). The wildtype itself has 7 fixed points (vertical dashed
line) and no limit cycles.
5 Computational aspects
Computation of histograms over functional networks is fa-
cilitated by the fact that rows of functional network matri-
ces combine independently. The set of all functional matri-
ces M is a product over sets Mi of functional row vectors
for each node i. In fact, each row of the matrix has its
own neutral graph. The Cartesian product [15] of these is
the neutral graph of the whole system. As noted by Lau
et al. [14], the set of network matrices performing a given
state sequence has a simple combinatorial structure. One
can check independently for each node i if it takes the
required state at each time step t. The states taken by
node i only depend on the i-th row and not on the whole
matrix. Thus, a functional network can be constructed by
independently combined functional row vectors into a ma-
trix. The set Mi of functional row vectors is obtained by
testing each of the 311 ≈ 2 × 105 possible vectors over
{−1, 0, 1} for each node i.
For calculating the histogram of the number of in-
teractions over all functional networks (Fig. 4c), we first
calculate this histogram gi individually for each node i.
The number of functional vectors in row i with exactly x
interactions is
gi(x) = |{r ∈ Mi|k(r) = x}| (2)
where k(r) is the number of non-zero entries in a row
vector r ∈ Mi.
These N histograms are now iteratively combined into
the histogram hi over the first i rows according to
hi(z) =
z
∑
x=0
gi(x)hi−1(z − x) (3)
with initialization h
1
= g
1
. Then hN is the histogram over
the complete matrices having N rows. Histograms are ob-
tained analogously for other observables like the number
of positive and negative interactions and the Hamming
distance. In fact, the iterative combination of histograms

G. Boldhaus and K. Klemm: Regulatory networks and connected components of the neutral space... 237
is applicable for all matrix observables that are decompos-
able into observables of single rows.
The number of attractors and the size of basins (Fig. 5)
do not fall into this class of row-decomposable observables.
Thus sampling is used to obtain statistics of these observ-
ables. For drawing a matrix A from the uniform distri-
bution on M , we draw the i-th row vector of A from the
uniform distribution on Mi for all i ∈ {1, . . . , N}. Before
starting the sampling, the whole sets Mi need to be com-
puted and stored once. In the present case, these sets are
sufficiently small with less than 2× 104 elements for all i.
For each functional matrix A in the sample, basin sizes
and attractor numbers are obtained by complete enumer-
ation of the Boolean state space having 2N = 2048 ele-
ments.
6 Discussion and outlook
We have analyzed the neutral graph (also called neutral
network) of discrete regulatory networks reproducing the
cell cycle sequence of budding yeast [10]. The neutral
graph falls into many connected components. Networks in
different components of the neutral graph are not acces-
sible to each other through a sequence of mutations that
retains cell cycle functionality. Our finding contrasts with
the connected neutral graphs in the work by Ciliberti et
al. in a similar type of discrete regulatory networks [7,8].
There, function is defined as the eventual arrival at a pre-
defined fixed point from a given initial condition. In the
present study, the exact sequence of states leading to the
fixed point is part of the required phenotype. We hypothe-
size that the fragmentation of the neutral graph is caused
by increasing functional constraints.
Further analysis has revealed that functional networks
accessible from the empirical wildtype are structurally and
dynamically distinct from other functional networks. Net-
works in the wildtype component are more sparsely wired
and their dynamics is more resilient to perturbations, as
compared to the average of all functional networks.
Thus, networks in the wildtype component have prop-
erties similar to the wildtype itself. This is remarkable
since most networks in the wildtype component are dis-
tant from the wildtype, having only a few interactions in
common.
Future investigations could establish conditions for
the connectedness of the neutral graph. To what extent is
the fragmentation of the neutral graph caused by the
strong discretization of interaction strengths? Allowing
finer adaptations would lead to less fragmented neutral
graphs. In the extreme (though chemically unrealistic)
limit of continuously evolving interaction strengths, the
set of all functional network matrices is convex and thus
connected.
We thank Anke Busch, Nadine Menzel, and Markus Riester
for valuable comments on the manuscript. This work has been
funded by the VolkswagenStiftung.
References
1. M. Kimura, The Neutral Theory of Molecular Evolution
(Cambridge University Press, Cambridge, UK, 1983)
2. A. Wagner, Robustness and Evolvability in Living Systems
(Princeton University Press, 2005)
3. P. Schuster, W. Fontana, P.F. Stadler, I.L. Hofacker, Proc.
Roy. Soc. Lond. B 255, 279 (1994)
4. A. Babajide, I.L. Hofacker, M.J. Sippl, P.F. Stadler, Fold
Des. 2, 261 (1997)
5. E. Davidson, M. Levin, Proceedings of the National
Academy of Sciences of the United States of America 102,
4935 (2005)
6. S. Bornholdt, K. Sneppen, Phys. Rev. Lett. 81, 236 (1998)
7. S. Ciliberti, O.C. Martin, A. Wagner, Proceedings of the
National Academy of Sciences of the United States of
America 104, 13591 (2007)
8. S. Ciliberti, O.C. Martin, A. Wagner, PLoS Computational
Biology 3, e15 (2007)
9. F. Stauffer, J. Berg, EPL 88, 48004 (2009)
10. F. Li, T. Long, Y. Lu, Q. Ouyang, C. Tang, Proceedings
of the National Academy of Sciences of the United States
of America 101, 4781 (2004)
11. P.T. Spellman, G. Sherlock, M.Q. Zhang, V.R. Iyer,
K. Anders, M.B. Eisen, P.O. Brown, D. Botstein,
B. Futcher, Mol. Biol. Cell. 9, 3273 (1998)
12. S. Kauffman, J. Theor. Biol. 22, 437 (1969)
13. B. Drossel, Reviews of Nonlinear Dynamics and
Complexity (Wiley-VCH, 2008), Chap. Random Boolean
Networks, Vol. 1, pp. 69–99
14. K.Y. Lau, S. Ganguli, C. Tang, Phys. Rev. E 75, 051907
(2007)
15. W. Imrich, Product Graphs: Structure And Recognition
(Wiley Interscience Series in Discrete Mathematics, 2000)