Emergent Criticality from Co-evolution in Random Boolean Networks
Min Liu
The co-evolution of network topology and dynamics is studied in an evolutionary Boolean network model that is a simple model of gene regulatory network. We find that a critical state emerges spontaneously resulting from interplay between topology and dynamics during the evolution. The final evolved state is shown to be independent of initial conditions. The network appears to be driven to a random Boolean network with uniform in-degree of two in the large network limit. However, for biologically realized network sizes, significant finite-size effects are observed including a broad in-degree distribution and an average in-degree connection between two and three. These results may be important for explaining properties of gene regulatory networks.
paper.tex
% ****** Start of file apssamp.tex ******
%
%   This file is part of the APS files in the REVTeX 4 distribution.
%   Version 4.0 of REVTeX, August 2001
%
%   Copyright (c) 2001 The American Physical Society.
%
%   See the REVTeX 4 README file for restrictions and more information.
%
% TeX'ing this file requires that you have AMS-LaTeX 2.0 installed
% as well as the rest of the prerequisites for REVTeX 4.0
%
% See the REVTeX 4 README file
% It also requires running BibTeX. The commands are as follows:
%
%  1)  latex apssamp.tex
%  2)  bibtex apssamp
%  3)  latex apssamp.tex
%  4)  latex apssamp.tex
%
\documentclass[letter,twocolumn,showpacs,preprintnumbers,amsmath,amssymb]{revtex4}
%\documentclass[preprint,showpacs,preprintnumbers,amsmath,amssymb]{revtex4}

% Some other (several out of many) possibilities
%\documentclass[preprint,aps]{revtex4}
%\documentclass[preprint,aps,draft]{revtex4}
%\documentclass[prb]{revtex4}% Physical Review B

\usepackage{graphicx}% Include figure files
\usepackage{dcolumn}% Align table columns on decimal point
\usepackage{bm}% bold math

%\nofiles

\begin{document}

%\preprint{APS/123-QED}

\title{Emergent Criticality from Co-evolution in Random Boolean Networks}
% Force line breaks with \\
\author{Min Liu}
\email{mliu3@uh.edu}
 \author{Kevin E. Bassler}%
 \email{bassler@uh.edu}
\affiliation{%
Department of Physics, University of Houston, Houston, TX 77204-5005
}%


\date{\today}% It is always \today, today,
             %  but any date may be explicitly specified

\begin{abstract}
The co-evolution of network topology and dynamics is studied in an
evolutionary Boolean network model that is a
simple model of gene regulatory network. We find that a critical
state emerges spontaneously resulting from interplay between
topology and dynamics during the evolution. The final evolved state
is shown to be independent of initial conditions. The network
appears to be driven to a random Boolean network with uniform
in-degree of two in the large network limit. However, for
biologically realized network sizes, significant finite-size effects
are observed including a broad in-degree distribution and an average
in-degree connection between two and three. These results may be
important for explaining properties of gene regulatory networks.
\end{abstract}

\pacs{87.23.Kg,  89.75.Hc, 05.65.+b, 87.15.Aa}% PACS, the Physics and Astronomy
                             % Classification Scheme.
%\keywords{Suggested keywords}%Use showkeys class option if keyword
                              %display desired
\maketitle

\section{Introduction}
Boolean networks\cite{Kauffman1,Kauffman2,RBN1,RBN2,RBN5,RBN4,RBN3,SOCRBN}
have been extensively studied over the past three decades. They have
applications as models of gene regulatory networks, and also as
models of social and economic systems. 
As ``coarse-grained'' models of genetic networks
they aim to capture much of the observed systematic behavior of the networks while 
simplifying local gene expression to a binary (on/off) state\cite{lessmore}.
Despite their simplicity, recent work has demonstrated that the model
can indeed predict the essential features of the dynamics of a biological
genetic circuit\cite{lessE1,lessE2}. An important feature of
Boolean networks is that they have a continuous
phase transition between so-called ordered and chaotic phases. It has been
argued\cite{edge2,Kauffman3} that gene regulatory networks of living
systems should be at or close to criticality, at the so-called
``edge of chaos" between the two phases, because then they can maintain both the
evolvability and stability. 

Many studies of critical random Boolean networks(RBNs) have considered networks
with homogeneous topology
in which each node has the same number of inputs from other
nodes~\cite{Parisi,drossel2,drossel3,drossel4,RBN6,Socolar,economics,canalization2,canalization}.
Accumulating experimental
evidence\cite{experiment0,experiment1,experiment3}, however, shows
that real genetic networks do not have homogeneous connectivity,
but, instead, are topologically heterogeneous. This diversity of
architecture is, presumably, of great importance for the stability
of living cells. Studies of RBNs with heterogeneous topology have
analytically determined the location of the ordered to chaotic phase
transition in the large network limit and also
demonstrated how different kinds of topology affects the
stability of the dynamics\cite{analytic,Fox,Aldana1,Aldana2}. 
These studies emphasized
the importance of criticality and the influence of the network's
topology on its dynamics, but they did not attempt to explain
how critical networks with heterogeneous topology come to exist.

Recent analysis of real gene regulatory
networks\cite{experiment2} has uncovered the possibility that
the interactions between genes can change
in response to diverse stimuli.
The resulting changes in the network topology
can be far greater than what is expected simply from random mutation.
However, few general
principles are known about the evolution of network
topology. In an effort to determine what some of those principles may be
Bornholdt and Rohlf\cite{Bornholdt}
studied a simple model of neural networks known as
random threshold networks (RTNs). 
In their study, the topology of the RTN evolved according to a rule that
depends on the local dynamics of the network.
According to their rule, active nodes, whose binary state changes in time, tend to lose
links, while inactive node, whose binary state is fixed, tend to gain new links.
They observed that with this rule for changing network topology,
in the limit of large networks, 
the RTNs evolved to a critical network with average connectivity
$\overline{K}=2$. The study, therefore, discovered an
interdependence between
a network's dynamics and its topology.  

Motivated by 
these recent findings concerning gene regulatory networks and the behavior of
RTNs, here we
investigate 
the effect of a similar evolutionary rule on a simple
model of a genetic regulatory system.
In particular, 
we study an
evolutionary RBN model.
However, there is a fundamental difference between RTNs and RBNs.
In RTNs the dynamics of each node is controlled by the same threshold function, while in
RBNs the dynamics of each node is controlled by a randomly chosen Boolean function.
Therefore, in addition to evolving the topology of the network,
we also allow the Boolean functions used by the nodes to change. 
In the context of a model of gene regulation the rule we use to evolve the network topology 
assumes that
there is some selection pressure on an individual gene due to its activity. The
evolution causes genes that are in a frozen state of regulation, and, thus, are almost
nonfunctional, to gain functionality, while it reduces the functionality of genes that
are actively regulated. 
Thus this study investigates the co-evolution of network structure and network
dynamics.
We show that, independent of the initial topology of the
network, the RBN evolves to a critical network with a finite number of
nodes and which has a heterogeneous topology. Two different variants
of our model are considered, and our principal conclusions are the
same for both of them. Perhaps our most important result concerns
the finite-size effects of the model. As the size of the network
increases, the distribution of in-degree connectivity becomes
increasingly narrow and sharply peaked at a value of $K=2$. However,
for biologically realized network sizes, we show that the final
evolved critical state has a broad distribution of in-degree
connectivity with an average value between 2 and 3. Both of these
features also occur in real networks, suggesting that many of the
topological features of real networks may be due to their finite
size.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}[t]
\includegraphics[scale=0.52]{RBN.eps}% Here is how to import EPS art
\caption{\label{rbn} A directed graph $G(7,3)$ represents a
homogeneous RBN with 7 nodes and in-degree connectivity of 3. Black
and white represent binary states `1' and `0' respectively. The
state vector of network is $\mathbf{\Sigma}(t) = (0,1,0,0,1,1,0)$.
The arrow on each link indicates the direction of information flow
in the network. }
\end{figure}
\section{Definition of Model}

\subsection{Dynamics of Random Boolean Networks }

A generalized RBN consists of $N$ randomly interconnected nodes,
$i=1, \cdots, N$, each of which has $K_i$ in-degree connections from
nodes that regulate its behavior. The simplest Boolean network model
is a homogeneous RBN in which each node has a same number of input
nodes. In this case, the connections between nodes is described by a
random directed graph $G(N,K)$ consisting of $N$ nodes with uniform
in-degree $K$. Figure~\ref{rbn} illustrates an example of $G(7,3)$.
Each node has a Boolean dynamical state at time $t$,  $\sigma_i(t) =
0 $ or $1$. The state of each node at time $t+1$ is a function of
all states of its $K_i$ regulatory nodes at time $t$. Hence, the
discrete dynamics of the network is given by
\begin{equation}\label{dynamics}
\sigma_i(t+1) = f_i\left(\sigma_{i_1}(t),\sigma_{i_2}(t),\cdots
,\sigma_{i_{K_i}}(t)\right)
\end{equation}
where $i_1$,$i_2$,$\cdots$,$i_{K_i}$ are those input nodes
regulating node $i$. The function $f_i$ is a Boolean function of
$K_i$ variables that determines the output of node $i$ for all of
the $2^{K_i}$ possible sets of input. 
Note that random Boolean functions can be generated with an ``interaction bias'' $p$
by setting the output value of the function for each set of input to be $1$ with probability
$p$.  
The bias $p$ can be interpreted as a biochemical reaction parameter.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Talk about the attractor and state space next.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Given the Boolean state of each node $i$ at time $t$, $\sigma_i(t)$,
the state vector of network is defined as 
$\mathbf{\Sigma}(t)=(\sigma_1(t),\cdots,\sigma_N(t))$. 
The path
that $\mathbf{\Sigma}(t)$ takes over time $t$ is a dynamical
trajectory in the phase space of system.  
Because the
dynamics defined in Eq.~\ref{dynamics} is deterministic and the
phase space is finite, all dynamical trajectories eventually
become periodic. That is, after some possible transient behavior,
each trajectory will repeat itself forming a cycle given by,
\begin{equation}\label{attractor}
\mathbf{\Sigma}(t) = \mathbf{\Sigma}(t + \Gamma).
\end{equation}
The periodic part of the trajectory is the attractor of the
dynamics, and the minimum $\Gamma > 0$ that satisfies
Eq.~\ref{attractor} is the period of the attractor.


%In general, multiple attractors of different periods coexist in the
%phase space. The initial state $\mathbf{\Sigma}(0)$ determines which
%attractor the network finally reaches. The set of all initial states
%reaching the same attractor forms a basin in the phase space. 
Two phases exist in RBNs, chaotic and ordered, characterized by their
dynamical behavior\cite{RBN1,RBN2}. One important way of
distinguishing the two phases is to measure the
distribution of its attractor periods beginning with random initial
states. For RBNs in the chaotic phase the distribution of attractor
periods is sharply peaked near an average value that grows
exponentially with system size $N$, and for RBNs in the ordered
phase the distribution of attractor periods is sharply peaked near
an average value that is nearly independent of $N$. Critical RBNs,
however, have a broad power law distribution of attractor
periods\cite{Parisi}.


\subsection{Co-evolution in Random Boolean Networks}
In our model, the evolutionary changes of the topology of the
network are driven by the dynamics of the network, and the functions
that control the dynamics of network simultaneously evolve due to
the changes in the network topology. Thus, this study investigates
the co-evolution of network topology and dynamics. Similar to the
one used in Ref.~\cite{Bornholdt}, the topology-evolving rule is
simply that a frozen gene grows a link while an active gene loses a
link. The dynamical functions can be changed in either an annealed
or a quenched way. The detailed algorithm is defined as follows:
\smallskip
\begin{enumerate}
    \item Start with a homogeneous RBN, $G(N,K_0)$ with uniform in-degree
    connectivity $K_i=K_0$ for all $N$, and generate a random Boolean function $f_i$ for each node $i$.
    \item Choose a random initial system state $\mathbf{\Sigma}(0)$.
    Update the state using Eq.~\ref{dynamics}
    and find the dynamical attractor. See the appendix for a description
of the algorithm used to find the attractor.
    \item Choose a node $i$ at random and determine its average
    activity $\overline{O}(i)$ over the attractor.
    \begin{equation}\label{criterion}
    \overline{O}(i) = \frac{1}{\Gamma} \sum_{t=T}^{T+\Gamma-1} \sigma_i(t)
\end{equation}
where $T$ is a time large enough so that the periodic attractor has
been reached, and $\Gamma$ is the period of the attractor. If
$\overline{O}(i)=1$ or $0$, then its state does not change over the
duration of the attractor; it is frozen. Alternatively, if
$0<\overline{O}(i)<1$, then node $i$ is active during the attractor.
\item Change the network topology by rewiring the connections to the
node chosen in the previous step. If it is frozen, then a new incoming
link from a randomly selected node $j$ is added to it. If it is
active, then one of its existing links is randomly selected and removed.
Note that this rewiring changes $K_i$.
\item The Boolean functions of network are regenerated.
Two different methods have been used:
\begin{itemize}
    \item Annealed model: A new Boolean function is generated for every
    node of the network.
    \item Quenched model: A new Boolean function is generated only for the chosen
   node $i$, while the others remain what they were previously.
\end{itemize}

\item Return to step 2.
\end{enumerate}

The time scale for an evolutionary change of the networks, steps
$2-6$ above, is called an epoch. As we will see, using this rewiring
rule the network topology evolves from a homogeneous one to a
heterogeneous one.
For simplicity, all random Boolean functions are generated with $p=1/2$, and 
therefore all Boolean functions with the same in-degree are equally likely to be generated. 

\begin{figure}[t]
\includegraphics[scale=0.35]{Caverage_K.eps}% Here is how to import EPS art
\caption{\label{AverageK}(color online) (a). Evolution of the ensemble averaged
in-degree connectivity in the annealed model for networks of size
$N=30$. The networks in each ensemble
initially start from different uniform connectivity,
$K_0 = 2, 3, 4$, and 5, but reach a same
statistical steady state $\langle \overline{K}\rangle=3.06$. Each
ensemble contained 15,000 realizations of the
network. (b). Evolution of ensemble averaged in-degree connectivity
for networks of three different size $N=30,50$, and 100 in the annealed
model.}
\end{figure}

\section{Simulation and Results }
We have simulated both the annealed and the quenched variants of the
model. Both variants give very similar results, and our principal
findings are the same for both variants. Therefore,  we will present
here mainly results from the annealed variant. Graph (a) of
Fig.~\ref{AverageK} shows the evolution of the average in-degree
connectivity $\overline{K}=\frac{1}{N}\sum_{i=1}^N K_i$ for networks
of size $N=30$ in the annealed variant of the model. Four curves are
shown. They show the results obtained by beginning with networks
with different uniform connectivity $K_0=2$, 3, 4, and 5. Each curve
is the average of 15,000 independent realizations of the network
evolution. This ensemble average is indicated by the angular
brackets. Each different realization in an ensemble begins with a
different random network and with a different random initial state
vector. Remarkably, despite the difference in initial conditions,
all four curves collapse after about 200 epochs, and they all
approach the same final statistical steady state that has an
average in-degree connectivity $\langle \overline{K}\rangle =
3.06$. This occurs without tuning and suggests that the final
evolved topology of the network is independent of the initial
topology of the network. The steady state value of $\langle \overline{K}
\rangle$ depends on the size of the system as shown in graph (b) of
Fig.~\ref{AverageK}. Starting with networks that all have the same
initial uniform connectivity $K_0=4$, but which have different size
$N=30,50$, and 100, we find that larger networks evolve to steady states
with smaller values of $\langle \overline{K} \rangle$. Very
similar results are obtained for the quenched version of the model.
For example, for networks with $N=30$ the steady state value of the
average connectivity is $\langle \overline{K}\rangle=3.08$.

\begin{figure}[t]
\includegraphics[scale=0.32]{distn200.eps}% Here is how to import EPS art
\caption{\label{degree200}(color online) Distribution of in-degree (square) and
out-degree (circle) connectivity in the annealed model. The size of
the networks is $N=200$.} \vspace{0.4in}
\includegraphics[scale=0.32]{Indist.eps}% Here is how to import EPS art
\caption{\label{degreesize}(color online) Distribution of in-degree connectivity
in the annealed model. The network sizes are $N=200$, $400$, and
$1000$.}
\end{figure}
We have also calculated the in-degree and out-degree connectivity
distributions, $P(K_{in})$ and $P(K_{out})$, of the evolved RBNs in
the steady state. Initially, all nodes of the network have a uniform
in-degree $K_0$, meaning that the in-degree distribution is a
discrete delta function $P(K_{in}) = \delta_{K_{in},K_0}$, and the
out-degree distribution is a binomial distribution. However, through
the evolutionary rewiring of the network both the in-degree and
out-degree distributions change. The in-degree and out-degree
distributions in the steady state of the annealed version, are shown
in Fig.~\ref{degree200} for $N=200$. They are both right skewed
bell-shaped distributions peaking at $K=2$. The out-degree
distribution remains a binomial distribution but the average
connectivity changes. The in-degree distribution, although it has
the same average connectivity as the out-degree distribution, is
more sharply peaked. As the size of the network grows, the in-degree
distribution becomes increasingly narrow and peaked at the value
$K_{in}=2$, as shown in Fig.~\ref{degreesize}. Based on this
observation, we conjecture that the distribution tends to converge
into a discrete delta function $\delta_{K_{in}, 2}$ in the large
network limit $N\rightarrow \infty$, indicating that the network
becomes a homogeneous RBN in that limit.

In order to probe the dynamical nature of evolved steady states we
computed the distribution $P(\Gamma)$ of steady state attractor
period $\Gamma$ in the ensemble of RBNs simulated. The distribution
has a broad, power-law behavior for both the annealed and quenched
variants of the model. Figure~\ref{powerlaw} shows the results for
networks with $N=200$. As long as $N$ is about 30 or larger, results
for other size networks are similar. As discussed above, this
power-law distribution indicates that the networks have critical
dynamics. Also in the figure, the straight line has a slope of 1.0.
Thus, the critical exponent describing the power-law is
approximately 1.0. This value of the exponent is obtained for all
system sizes studied in both the quenched and annealed versions of
the model. In short, we find that a robust criticality emerges in
the evolutionary RBNs.
\begin{figure}[t]
\includegraphics[scale=0.32]{Attractor_QA.eps}% Here is how to import EPS art
\caption{\label{powerlaw}(color online) Power law distribution of steady state
attractor period $\Gamma$ in both annealed (circle) and quenched
(square) models for $N=200$ system. The dashed straight line has a
slope of $1.0$. }
\end{figure}

Given the steady state value $\langle\overline{K} \rangle=2$ in
the large network limit $N\rightarrow\infty$, we studied the
finite-size effects in the model. As shown in Fig.~\ref{finitesize},
the values of $\langle \overline{K}(N)\rangle$ for finite $N$ obey
the scaling function
\begin{equation}\label{}
   \langle\overline{K}(N)\rangle - 2 = A N^{-\beta} .
\end{equation}
Fitting the data to this function, we find that the coefficient is
$A=2.50 \pm 0.06$ and the exponent is $\beta=0.264\pm 0.005$. Thus
the value of
$\langle\overline{K}(N)\rangle$ is always larger than 2 for finite
$N$. Note that steady state values of the average connectivity in
random threshold networks have a similar scaling form
\cite{Bornholdt}, but, in that case $A=12.4\pm 0.5$ and $\beta =
0.47 \pm 0.01$.

\section{Discussion and Conclusions}
The mechanism we find here that leads to the emergent critical state
has some similarity to self-organized criticality
(SOC)\cite{SOC1,SOC2}, but is different. SOC is the tendency of
driven dissipative dynamical systems to organize themselves into a
critical state far from equilibrium through avalanches of activity
of all sizes. In our particular model, the evolutionary Boolean
network is a dissipative dynamical system because multiple different
states may map into the same attractor so that information is lost.
Similar to SOC systems, our model is driven subject to two competing
rules, and the network organizes itself into a steady state that
results from a dynamical balance of the competition between those
rules. Moreover, the critical state is robust irrespective of
initial connectivity in both the quenched or annealed versions of
the model. The emergent critical state acts like a global attractor
in the evolution process. However, unlike the mechanism of
traditional SOC, but similar to the mechanisms that have been shown to
lead to criticality in random threshold networks \cite{Bornholdt} and
in homogeneous Boolean networks \cite{economics}, the
self-organizing mechanism here is based on a topological phase
transition in dynamical networks.
\begin{figure}[t]
\includegraphics[scale=0.32]{finitesize.eps}% Here is how to import EPS art
\caption{\label{finitesize}(color online) Finite-size effects in the annealed
model. The data shown are for systems of different sizes $N=30,
50, 80, 100, 150, 200$, and 400. The dashed straight line has a slop of
$0.26$.}
\end{figure}

Our results indicate that, with the rewiring rules we use, networks
in the limit $N\rightarrow \infty$ will evolve to have a homogeneous
in-degree connectivity of 2. However, we find that finite size
networks evolve to have a broadly distributed heterogeneous
in-degree connectivity. The average in-degree connectivity of the
evolved networks is between 2 and 3 for biologically realized
network sizes. This result may be important for explaining the
observed structure of real gene regulatory networks. Real genetic
networks have a number of nodes $N$ ranging from near 100 to
thousands, and typically have a heterogeneous connectivity with an
average in-degree slightly larger than 2. For example, a recent
experiment studying the gene regulatory network of \emph{S.
cerevisiae} \cite{experiment2} found that the network contains 3420
genes and has an average in-degree connectivity $\overline{K}_{in} =
2.1$.


Finally, we note that real genetic networks exhibit an approximately
scale-free out-degree distribution while the in-degree distribution
is exponentially decaying\cite{experiment2,Albert}. In our numerical
study, we similarly obtain an in-degree distribution that decays
faster than the out-degree distribution, but our model does not
produce a scale-free like out-degree distribution. Therefore, in
order to be more realistic the model needs to be extended by adding
other factors that will capture this feature. We do not believe
though that such extensions will alter the principal conclusions of
this paper concerning the importance of finite size effects in the
evolution of network topology in real genetic networks.

\begin{acknowledgments}
We thank Maximino Aldana, Julio Monte, Toshimori Kitami  and Adam
Sheya for stimulating discussions. This work was supported
by the NSF through grant \#DMR-0427538, and by SI International
through the AFRL under contract \#FA8756-04-C-0258.
\end{acknowledgments}

\appendix*
\section{}
In the work presented here and in our previous studies\cite{economics,canalization,canalization2}
the following algorithm was used to determine the dynamical attractor.

Begin by using the initial state
$\mathbf{\Sigma}(0)$
as the ``checkpoint'' state.
Then for each time $T$, $0<T \leq T_1$, the state is updated using Eq.~1.
If 
$\mathbf{\Sigma}(T) = \mathbf{\Sigma}(0)$, 
then the attractor is found, it has period $T$,
and the search ends.
Note that the new state is compared only to the
checkpoint state and no other previous states.

If after $T_1$ updates no attractor is found then the checkpoint state
is changed to
$\mathbf{\Sigma}(T_1)$.
For each time $T$, $T_1<T\leq T_2$, the state is again updated using Eq.~1.
If 
$\mathbf{\Sigma}(T) = \mathbf{\Sigma}(T_1)$, 
then the attractor is found, it has period $T-T_1$,
and the search ends.

If after $T_2$ updates the attractor has still not been found then the checkpoint state
is changed to
$\mathbf{\Sigma}(T_2)$.
For each time $T$, $T_2<T\leq T_2 + T_{max}$, the state is again updated using Eq.~1.
If 
$\mathbf{\Sigma}(T) = \mathbf{\Sigma}(T_2)$, 
then the attractor is found, it has period $T-T_2$,
and the search ends.

Finally, if after $T_2+T_{max}$ updates the attractor has still not been found then the state
$\mathbf{\Sigma}(T_2+T_{max})$
is used as the new checkpoint state.
If 
$\mathbf{\Sigma}(T) = \mathbf{\Sigma}(T_2+T_{max})$, 
then the attractor is found, it has period $T-T_2-T_{max}$,
and the search ends.

If no attractor is found after this procedure, then we stop. In this case, the average output
state in Eq.~3 is calculated assuming that the attractor period is $\Gamma = T_{max}$ and that
the final checkpoint state
$\mathbf{\Sigma}(T_2+T_{max})$
is on the attractor.

This algorithm finds all attractors that have period less than or
equal to $T_{max}$ and that have a transient time to reach the
attractor from the initial state less than or equal to
$T_2+T_{max}$. For the results presented in this paper we used
$T_1=100$, $T_2=1000$, and $T_{max}=100 000$.

\bibliography{citation}% Produces the bibliography via BibTeX.

\end{document}

paper.bbl
\begin{thebibliography}{34}
\expandafter\ifx\csname natexlab\endcsname\relax\def\natexlab#1{#1}\fi
\expandafter\ifx\csname bibnamefont\endcsname\relax
  \def\bibnamefont#1{#1}\fi
\expandafter\ifx\csname bibfnamefont\endcsname\relax
  \def\bibfnamefont#1{#1}\fi
\expandafter\ifx\csname citenamefont\endcsname\relax
  \def\citenamefont#1{#1}\fi
\expandafter\ifx\csname url\endcsname\relax
  \def\url#1{\texttt{#1}}\fi
\expandafter\ifx\csname urlprefix\endcsname\relax\def\urlprefix{URL }\fi
\providecommand{\bibinfo}[2]{#2}
\providecommand{\eprint}[2][]{\url{#2}}

\bibitem[{\citenamefont{Kauffman}(1969)}]{Kauffman1}
\bibinfo{author}{\bibfnamefont{S.~A.} \bibnamefont{Kauffman}},
  \bibinfo{journal}{J. Theor. Biol.} \textbf{\bibinfo{volume}{22}},
  \bibinfo{pages}{437} (\bibinfo{year}{1969}).

\bibitem[{\citenamefont{Kauffman}(1990)}]{Kauffman2}
\bibinfo{author}{\bibfnamefont{S.~A.} \bibnamefont{Kauffman}},
  \bibinfo{journal}{Physica D} \textbf{\bibinfo{volume}{42}},
  \bibinfo{pages}{135} (\bibinfo{year}{1990}).

\bibitem[{\citenamefont{Derrida and Pomeau}(1986)}]{RBN1}
\bibinfo{author}{\bibfnamefont{B.}~\bibnamefont{Derrida}} \bibnamefont{and}
  \bibinfo{author}{\bibfnamefont{Y.}~\bibnamefont{Pomeau}},
  \bibinfo{journal}{EuroPhys. Lett.} \textbf{\bibinfo{volume}{1}},
  \bibinfo{pages}{45} (\bibinfo{year}{1986}).

\bibitem[{\citenamefont{Derrida and Weisbuch}(1986)}]{RBN2}
\bibinfo{author}{\bibfnamefont{B.}~\bibnamefont{Derrida}} \bibnamefont{and}
  \bibinfo{author}{\bibfnamefont{G.}~\bibnamefont{Weisbuch}},
  \bibinfo{journal}{J. Phys.} \textbf{\bibinfo{volume}{47}},
  \bibinfo{pages}{1297} (\bibinfo{year}{1986}).

\bibitem[{\citenamefont{Bastolla and Parisi}(1998)}]{RBN5}
\bibinfo{author}{\bibfnamefont{U.}~\bibnamefont{Bastolla}} \bibnamefont{and}
  \bibinfo{author}{\bibfnamefont{G.}~\bibnamefont{Parisi}},
  \bibinfo{journal}{Physica D} \textbf{\bibinfo{volume}{115}},
  \bibinfo{pages}{203} (\bibinfo{year}{1998}).

\bibitem[{\citenamefont{Albert and Barab\'{a}si}(2000)}]{RBN4}
\bibinfo{author}{\bibfnamefont{R.}~\bibnamefont{Albert}} \bibnamefont{and}
  \bibinfo{author}{\bibfnamefont{A.~L.} \bibnamefont{Barab\'{a}si}},
  \bibinfo{journal}{Phys. Rev. Lett.} \textbf{\bibinfo{volume}{84}},
  \bibinfo{pages}{5660} (\bibinfo{year}{2000}).

\bibitem[{\citenamefont{Aldana et~al.}()\citenamefont{Aldana, Coppersmith, and
  Kadanoff}}]{RBN3}
\bibinfo{author}{\bibfnamefont{M.}~\bibnamefont{Aldana}},
  \bibinfo{author}{\bibfnamefont{S.}~\bibnamefont{Coppersmith}},
  \bibnamefont{and} \bibinfo{author}{\bibfnamefont{L.~P.}
  \bibnamefont{Kadanoff}}, \eprint{in \emph{Perspectives and Problems in
  Nonlinear Science}, edited by E. Kaplan, J.E. Marsden, and K.R. Sreenivasan
  (Springer-Verleg, Berlin, 2003)}.

\bibitem[{\citenamefont{Luque et~al.}(2001)\citenamefont{Luque, Ballesteros,
  and Muro}}]{SOCRBN}
\bibinfo{author}{\bibfnamefont{B.}~\bibnamefont{Luque}},
  \bibinfo{author}{\bibfnamefont{F.~J.} \bibnamefont{Ballesteros}},
  \bibnamefont{and} \bibinfo{author}{\bibfnamefont{E.~M.} \bibnamefont{Muro}},
  \bibinfo{journal}{Phys. Rev. E} \textbf{\bibinfo{volume}{63}},
  \bibinfo{pages}{051913} (\bibinfo{year}{2001}).

\bibitem[{\citenamefont{Bornholdt}(2005)}]{lessmore}
\bibinfo{author}{\bibfnamefont{S.}~\bibnamefont{Bornholdt}},
  \bibinfo{journal}{Science} \textbf{\bibinfo{volume}{310}},
  \bibinfo{pages}{449} (\bibinfo{year}{2005}).

\bibitem[{\citenamefont{Albert and Othmer}(2003)}]{lessE1}
\bibinfo{author}{\bibfnamefont{R.}~\bibnamefont{Albert}} \bibnamefont{and}
  \bibinfo{author}{\bibfnamefont{H.~G.} \bibnamefont{Othmer}},
  \bibinfo{journal}{J. Theor. Biol.} \textbf{\bibinfo{volume}{223}},
  \bibinfo{pages}{1} (\bibinfo{year}{2003}).

\bibitem[{\citenamefont{Li et~al.}(2004)\citenamefont{Li, Long, Lu, Quyang, and
  Tang}}]{lessE2}
\bibinfo{author}{\bibfnamefont{F.}~\bibnamefont{Li}},
  \bibinfo{author}{\bibfnamefont{T.}~\bibnamefont{Long}},
  \bibinfo{author}{\bibfnamefont{Y.}~\bibnamefont{Lu}},
  \bibinfo{author}{\bibfnamefont{Q.}~\bibnamefont{Quyang}}, \bibnamefont{and}
  \bibinfo{author}{\bibfnamefont{C.}~\bibnamefont{Tang}},
  \bibinfo{journal}{Proc. Natl. Acad. Sci. U.S.A.}
  \textbf{\bibinfo{volume}{101}}, \bibinfo{pages}{4781} (\bibinfo{year}{2004}).

\bibitem[{\citenamefont{Packard}()}]{edge2}
\bibinfo{author}{\bibfnamefont{N.~H.} \bibnamefont{Packard}}, \eprint{in
  \emph{Dynamics Patterns in Complex Systems}, edited by J. A. S. Kelso, A. J.
  Mandell, and M. F. Shlesinger(World Scientific, Singapore, 1988)}.

\bibitem[{\citenamefont{Kauffman}()}]{Kauffman3}
\bibinfo{author}{\bibfnamefont{S.~A.} \bibnamefont{Kauffman}},
  \eprint{\emph{The Origins of Order} (Oxford University Press, New York,
  1993)}.

\bibitem[{\citenamefont{Bastolla and Parisi}(1997)}]{Parisi}
\bibinfo{author}{\bibfnamefont{U.}~\bibnamefont{Bastolla}} \bibnamefont{and}
  \bibinfo{author}{\bibfnamefont{G.}~\bibnamefont{Parisi}},
  \bibinfo{journal}{J. Theor. Biol} \textbf{\bibinfo{volume}{187}},
  \bibinfo{pages}{117} (\bibinfo{year}{1997}).

\bibitem[{\citenamefont{Greil and Drossel}(2005)}]{drossel2}
\bibinfo{author}{\bibfnamefont{F.}~\bibnamefont{Greil}} \bibnamefont{and}
  \bibinfo{author}{\bibfnamefont{B.}~\bibnamefont{Drossel}},
  \bibinfo{journal}{Phys. Rev. Lett.} \textbf{\bibinfo{volume}{95}},
  \bibinfo{pages}{048701} (\bibinfo{year}{2005}).

\bibitem[{\citenamefont{Drossel}(2005)}]{drossel3}
\bibinfo{author}{\bibfnamefont{B.}~\bibnamefont{Drossel}},
  \bibinfo{journal}{Phys. Rev. E} \textbf{\bibinfo{volume}{72}},
  \bibinfo{pages}{016110} (\bibinfo{year}{2005}).

\bibitem[{\citenamefont{Drossel et~al.}(2005)\citenamefont{Drossel, Mihaljev,
  and Greil}}]{drossel4}
\bibinfo{author}{\bibfnamefont{B.}~\bibnamefont{Drossel}},
  \bibinfo{author}{\bibfnamefont{T.}~\bibnamefont{Mihaljev}}, \bibnamefont{and}
  \bibinfo{author}{\bibfnamefont{F.}~\bibnamefont{Greil}},
  \bibinfo{journal}{Phys. Rev. Lett.} \textbf{\bibinfo{volume}{94}},
  \bibinfo{pages}{088701} (\bibinfo{year}{2005}).

\bibitem[{\citenamefont{Samuelsson and Troein}(2003)}]{RBN6}
\bibinfo{author}{\bibfnamefont{B.}~\bibnamefont{Samuelsson}} \bibnamefont{and}
  \bibinfo{author}{\bibfnamefont{C.}~\bibnamefont{Troein}},
  \bibinfo{journal}{Phys. Rev. Lett.} \textbf{\bibinfo{volume}{90}},
  \bibinfo{pages}{098701} (\bibinfo{year}{2003}).

\bibitem[{\citenamefont{Socolar and Kauffman}(2003)}]{Socolar}
\bibinfo{author}{\bibfnamefont{J.~E.~S.} \bibnamefont{Socolar}}
  \bibnamefont{and} \bibinfo{author}{\bibfnamefont{S.~A.}
  \bibnamefont{Kauffman}}, \bibinfo{journal}{Phys. Rev. Lett.}
  \textbf{\bibinfo{volume}{90}}, \bibinfo{pages}{068702}
  (\bibinfo{year}{2003}).

\bibitem[{\citenamefont{Paczuski et~al.}(2000)\citenamefont{Paczuski, Bassler,
  and Corral}}]{economics}
\bibinfo{author}{\bibfnamefont{M.}~\bibnamefont{Paczuski}},
  \bibinfo{author}{\bibfnamefont{K.~E.} \bibnamefont{Bassler}},
  \bibnamefont{and} \bibinfo{author}{\bibfnamefont{A.}~\bibnamefont{Corral}},
  \bibinfo{journal}{Phys. Rev. Lett.} \textbf{\bibinfo{volume}{84}},
  \bibinfo{pages}{3185} (\bibinfo{year}{2000}).

\bibitem[{\citenamefont{Bassler et~al.}(2004)\citenamefont{Bassler, Lee, and
  Lee}}]{canalization2}
\bibinfo{author}{\bibfnamefont{K.~E.} \bibnamefont{Bassler}},
  \bibinfo{author}{\bibfnamefont{C.}~\bibnamefont{Lee}}, \bibnamefont{and}
  \bibinfo{author}{\bibfnamefont{Y.}~\bibnamefont{Lee}},
  \bibinfo{journal}{Phys. Rev. Lett.} \textbf{\bibinfo{volume}{93}},
  \bibinfo{pages}{038101} (\bibinfo{year}{2004}).

\bibitem[{\citenamefont{Bassler and Liu}()}]{canalization}
\bibinfo{author}{\bibfnamefont{K.~E.} \bibnamefont{Bassler}} \bibnamefont{and}
  \bibinfo{author}{\bibfnamefont{M.}~\bibnamefont{Liu}}, \eprint{in \emph{Noise
  in Complex Systems and Stochastic Dynamics III, Proc. of SPIE Vol. 5845},
  edited by L. B. Kish, K. Lindenberg, and Z. Gingl (SPIE, Bellingham, WA,
  2005), P. 104}.

\bibitem[{\citenamefont{Thieffry et~al.}(1998)\citenamefont{Thieffry, Huerta,
  Perez-Rueda, and Collado-Vides}}]{experiment0}
\bibinfo{author}{\bibfnamefont{D.}~\bibnamefont{Thieffry}},
  \bibinfo{author}{\bibfnamefont{A.~M.} \bibnamefont{Huerta}},
  \bibinfo{author}{\bibfnamefont{E.}~\bibnamefont{Perez-Rueda}},
  \bibnamefont{and}
  \bibinfo{author}{\bibfnamefont{J.}~\bibnamefont{Collado-Vides}},
  \bibinfo{journal}{Bioessays} \textbf{\bibinfo{volume}{20}},
  \bibinfo{pages}{433} (\bibinfo{year}{1998}).

\bibitem[{\citenamefont{Lee et~al.}(2002)\citenamefont{Lee, Rinaldi, Robert,
  Odom, and \emph{et al.}}}]{experiment1}
\bibinfo{author}{\bibfnamefont{T.~I.} \bibnamefont{Lee}},
  \bibinfo{author}{\bibfnamefont{N.~J.} \bibnamefont{Rinaldi}},
  \bibinfo{author}{\bibfnamefont{F.}~\bibnamefont{Robert}},
  \bibinfo{author}{\bibfnamefont{D.~T.} \bibnamefont{Odom}}, \bibnamefont{and}
  \bibinfo{author}{\bibfnamefont{Z.~B.~J.} \bibnamefont{\emph{et al.}}},
  \bibinfo{journal}{Science} \textbf{\bibinfo{volume}{298}},
  \bibinfo{pages}{799} (\bibinfo{year}{2002}).

\bibitem[{\citenamefont{Tong et~al.}(2004)\citenamefont{Tong, Lesage, Bader,
  Ding, and \emph{et al.}}}]{experiment3}
\bibinfo{author}{\bibfnamefont{A.~H.~Y.} \bibnamefont{Tong}},
  \bibinfo{author}{\bibfnamefont{G.}~\bibnamefont{Lesage}},
  \bibinfo{author}{\bibfnamefont{G.~D.} \bibnamefont{Bader}},
  \bibinfo{author}{\bibfnamefont{H.}~\bibnamefont{Ding}}, \bibnamefont{and}
  \bibinfo{author}{\bibfnamefont{H.~X.} \bibnamefont{\emph{et al.}}},
  \bibinfo{journal}{Science} \textbf{\bibinfo{volume}{303}},
  \bibinfo{pages}{808} (\bibinfo{year}{2004}).

\bibitem[{\citenamefont{Luque and Sol\'e}(1997)}]{analytic}
\bibinfo{author}{\bibfnamefont{B.}~\bibnamefont{Luque}} \bibnamefont{and}
  \bibinfo{author}{\bibfnamefont{R.~V.} \bibnamefont{Sol\'e}},
  \bibinfo{journal}{Phys. Rev. E} \textbf{\bibinfo{volume}{55}},
  \bibinfo{pages}{257} (\bibinfo{year}{1997}).

\bibitem[{\citenamefont{Fox and Hill}(2001)}]{Fox}
\bibinfo{author}{\bibfnamefont{J.~J.} \bibnamefont{Fox}} \bibnamefont{and}
  \bibinfo{author}{\bibfnamefont{C.~C.} \bibnamefont{Hill}},
  \bibinfo{journal}{Chaos} \textbf{\bibinfo{volume}{11}}, \bibinfo{pages}{809}
  (\bibinfo{year}{2001}).

\bibitem[{\citenamefont{Aldana}(2003)}]{Aldana1}
\bibinfo{author}{\bibfnamefont{M.}~\bibnamefont{Aldana}},
  \bibinfo{journal}{Physica D} \textbf{\bibinfo{volume}{185}},
  \bibinfo{pages}{45} (\bibinfo{year}{2003}).

\bibitem[{\citenamefont{Aldana and Cluzel}(2003)}]{Aldana2}
\bibinfo{author}{\bibfnamefont{M.}~\bibnamefont{Aldana}} \bibnamefont{and}
  \bibinfo{author}{\bibfnamefont{P.}~\bibnamefont{Cluzel}},
  \bibinfo{journal}{Proc. Natl. Acad. Sci. U.S.A}
  \textbf{\bibinfo{volume}{100}}, \bibinfo{pages}{8710} (\bibinfo{year}{2003}).

\bibitem[{\citenamefont{Luscombe et~al.}(2004)\citenamefont{Luscombe, Babu, Yu,
  Snyder, and Gerstein}}]{experiment2}
\bibinfo{author}{\bibfnamefont{N.~M.} \bibnamefont{Luscombe}},
  \bibinfo{author}{\bibfnamefont{M.~M.} \bibnamefont{Babu}},
  \bibinfo{author}{\bibfnamefont{H.}~\bibnamefont{Yu}},
  \bibinfo{author}{\bibfnamefont{M.}~\bibnamefont{Snyder}}, \bibnamefont{and}
  \bibinfo{author}{\bibfnamefont{S.~A. T.~M.} \bibnamefont{Gerstein}},
  \bibinfo{journal}{Nature} \textbf{\bibinfo{volume}{431}},
  \bibinfo{pages}{308} (\bibinfo{year}{2004}).

\bibitem[{\citenamefont{Bornholdt and Rohlf}(2000)}]{Bornholdt}
\bibinfo{author}{\bibfnamefont{S.}~\bibnamefont{Bornholdt}} \bibnamefont{and}
  \bibinfo{author}{\bibfnamefont{T.}~\bibnamefont{Rohlf}},
  \bibinfo{journal}{Phys. Rev. Lett.} \textbf{\bibinfo{volume}{84}},
  \bibinfo{pages}{6114} (\bibinfo{year}{2000}).

\bibitem[{\citenamefont{Bak et~al.}(1987)\citenamefont{Bak, Tang, and
  Wiesenfeld}}]{SOC1}
\bibinfo{author}{\bibfnamefont{P.}~\bibnamefont{Bak}},
  \bibinfo{author}{\bibfnamefont{C.}~\bibnamefont{Tang}}, \bibnamefont{and}
  \bibinfo{author}{\bibfnamefont{K.}~\bibnamefont{Wiesenfeld}},
  \bibinfo{journal}{Phys. Rev. Lett.} \textbf{\bibinfo{volume}{59}},
  \bibinfo{pages}{381} (\bibinfo{year}{1987}).

\bibitem[{\citenamefont{Bak et~al.}(1988)\citenamefont{Bak, Tang, and
  Wiesenfeld}}]{SOC2}
\bibinfo{author}{\bibfnamefont{P.}~\bibnamefont{Bak}},
  \bibinfo{author}{\bibfnamefont{C.}~\bibnamefont{Tang}}, \bibnamefont{and}
  \bibinfo{author}{\bibfnamefont{K.}~\bibnamefont{Wiesenfeld}},
  \bibinfo{journal}{Phys. Rev. A} \textbf{\bibinfo{volume}{38}},
  \bibinfo{pages}{364} (\bibinfo{year}{1988}).

\bibitem[{\citenamefont{Albert}(2005)}]{Albert}
\bibinfo{author}{\bibfnamefont{R.}~\bibnamefont{Albert}}, \bibinfo{journal}{J.
  Cell Sci.} \textbf{\bibinfo{volume}{118}}, \bibinfo{pages}{4947}
  (\bibinfo{year}{2005}).

\end{thebibliography}
