Robustness Leads Close to the Edge of Chaos in Coupled Map Networks: toward the understanding of biological networks

Dynamics in biological networks are in general robust against several perturbations. We investigate a coupled map network as a model motivated by gene regulatory networks and design systems which are robust against phenotypic perturbations (perturbations in dynamics), as well as systems which are robust against mutation (perturbations in network structure). To achieve such a design, we apply a multicanonical Monte Carlo method. Analysis based on the maximum Lyapunov exponent and parameter sensitivity shows that systems with marginal stability, which are regarded as systems at the edge of chaos, emerge when robustness against network perturbations is required. This emergence of the edge of chaos is a self-organization phenomenon and does not need a fine tuning of parameters.


Introduction
Complex dynamical behaviors on a network can be found in a variety of biological networks, such as gene regulatory networks, neural networks and food-web. Such systems share a common characteristic: observed dynamics are robust against disturbance introduced in its dynamics, as well as against disturbance in its network [1,2,3]. For example, gene expression patterns obtained thorough transcription-translation regulations are kept stable in spite of extrinsic noises (i.e., perturbations in dynamics) and mutations (i.e., perturbations in a network). It seems reasonable to think that robustness against environmental perturbations has been evolutionarily developed for adapting to noisy environments. There are also several advantages in having mutational robustness -it buffers against deleterious mutations. Recently, robustness in biological networks has attracted much attention of many researchers and has been thought to be one of fundamental properties of life [1,2,3,4].
This point of view naturally gives rise to the question: what kind of system emerges when only robustness is required ? The answer to this question will be helpful for understanding the design principle of living systems. In this paper, we investigate a coupled chaotic map network motivated by gene regulatory networks and show that systems at the edge of chaos are selected with only the requirement of robustness against network perturbations.
It has long been hypothesized that living systems favor the edge of chaos, where stability and chaoticity coexist. Originally, Kauffman [5] introduced the Boolean network model (N-K model) as a model of a gene regulatory network, and proposed the hypothesis that living systems prefer the edge of chaos because it allows systems to have complex behaviors [5]. Here we propose an alternative scenario, specifically that the requirement of having robustness against network perturbations drives living systems to the edge of chaos, regardless of whether or not staying at the edge of chaos is beneficial for living systems. In other wards, the edge of chaos can emerge as a byproduct of the robustness.

Model
We propose a coupled map system motivated by gene regulatory networks. Unlike the N-K model, each element in this model has its own dynamics. Assuming that is the gene expression of the i-th gene at time step t, the single gene dynamics are written as x t+1 i = G(x t i ). These dynamics mimic multiple processes in an expression of a single gene. In the presence of N genes, the dynamics of x t i are expressed as where ǫ is a coupling constant. W ij describes the strength of the interaction acting from gene j on gene i; and W ij ∈ [0, 1] satisfies both conditions W ii = 0 and N j W ij = 1 for each i. We call the matrix W , whose ij element is W ij , a network. Here we choose the logistic map g(x, a) = 1 − ax 2 as G(·). We use the model parameters a and ǫ as (a, ǫ) = (1.8, 0.1). This choice indicates that a single disconnected gene exhibits chaotic dynamics. A reason of this choice of G(·) is that a single gene expression is expected to be complex due to multiple processes underlying it We impose an additional constraint on W : the number of input links k to each gene is fixed. We note that in the case of W ij = 1/(N − 1) for all (i, j), the system becomes the globally coupled map (GCM) [6] and it shows highly chaotic behaviors at the parameters (a, ǫ) = (1.8, 0.1). In contrast to the N-K model, variables of this model take continuous values and a linear stability analysis can be applied. In this model, the network W is regarded as a genotype while the attractor of dynamics x t is regarded as a phenotype. Our goal is to design networks under the two different design principles: robustness against phenotypic perturbations (i.e., perturbations in the dynamics of gene x) and robustness against genotypic perturbations (i.e., perturbations on network W ). In both cases, only network W is tuned.

Design of Robust Network Against Perturbations in Dynamics
Let us start with the first design principle, namely robustness against perturbations in dynamics. In other words, we are aiming to design a system with a stable attractor. For this end, we use Lyapunov exponent analysis and a multicanonical Monte Carlo method [7,8].
Once a network W is given, the finite time maximum Lyapunov exponent λ 1 [9] is calculated for the dynamical system in Eq. (1), starting from a given initial state x 0 * i §. We perform the simulation up to T = 1500 and regard the first T ′ = 1000 steps as transient and discard them.
Our aim here is to sample networks with negative λ 1 , which indicates that dynamics of a network is stable. Such networks are expected to be rare for large N, because dynamics tend to be chaotic at the present parameters. We define the probability density of λ 1 as where δ is the Dirac δ-function, and p(W ) is the prior probability density that a network W appears under random sampling. We consider here a network ensemble in which p(W ) is a uniform distribution under the constraints of W ii = 0 and N j W ij = 1, given by where θ(z) is a function that satisfies θ(z) = 1 for z = 0 and θ(z) = 0 for z = 0. If a random sampling method is adopted in order to sample a network with "rare" value of λ 1 with probability D(λ 1 ), 1/D(λ 1 ) samples are required at least. If an annealing § Throughout this study, x 0 * i = sin(i) is used. We confirm that the choice of the initial state does not affect the results. method or an steepest descent method is adopted instead, we would obtain only a network with negative λ 1 but could not estimate D(λ 1 ), which plays a key role in the further analysis. Alternatively, we apply multicanonical Monte Carlo method [7,10,11], which has been used in fields of statistical physics, such as spin glass [12,13] and other studies [14,15]. This method allows us to sample networks with negative λ 1 efficiently and to estimate D(λ 1 ).
Our multicanonical Monte Carlo strategy adopted in this study is to perform random walks in λ 1 space by generating a Markov chain, where each step is biased inversely proportional to the probability D(λ 1 ), and thereby it enables us to obtain a flat histogram in λ 1 space, namely to equally sample λ 1 whose D(λ 1 ) are many orders of magnitude different. In order to generate the Markov chain in multicanonical Monte Carlo, a key quantity is the weight function w(λ 1 ) of λ 1 . If we have w(λ 1 ) that is inversely proportional to D(λ 1 ), networks with various λ 1 value are generated one after another, using the Markov process described in the Appendix (i). As a consequence, a uniform distribution of λ 1 (i.e., a flat histogram of λ 1 ) is obtained. We call these procedures as "random walk in λ 1 space". Details of the algorithm are given in the Appendix (i) and (ii). However, neither w(λ 1 ) nor D(λ 1 ) are known a priori. In this study, the Wang and Landau algorithm [12,13] is used to construct and to tune the weight function w(λ 1 ). Details of the implementation are given in Appendix (ii). Figure 1 shows the calculated densities D(λ 1 ) of λ 1 for the fixed input degree k = 2 − 5. Using density D(λ 1 ) in Fig. 1, the probability that networks with negative λ 1 are observed under random sampling is calculated by P (λ 1 < 0) = 0 λa D(λ 1 )dλ 1 . We estimate P (λ 1 < 0) with k = 2 − 5 and k = N − 1, which are shown in Fig. 2. Each P (λ 1 < 0) shows that a stable attractor becomes increasingly rare as N or k increases, indicating that these systems are in the chaotic phase (we define that a system is in the chaotic phase when only positive values of λ 1 appear as N → ∞). These results are consistent with the behavior of GCM with (a, ǫ) = (1.8, 0.1) [6].

Design of Robust Network Against Perturbations in Network
Using the second design principle, we design systems that are robust against genotypic perturbations (i.e., network perturbations). In other words, we design, using a multicanonical Monte Carlo method, networks W whose trajectory on the attractor hardly changes when a small network perturbation δW is added. We define the parameter sensitivity and use it as a guiding function of the robustness.
Sensitivity analysis using parameter sensitivity has been developed and applied in various fields [16,17]. While most of these studies have dealt with continuous time systems, we define the parameter sensitivity for discrete time systems as follows.
Let us denote the set of elements in W in Eq. (1) by a vector W . When a small network perturbation δW is introduced into W at t = T ′ , the displacement between unperturbed trajectory x t (W ) and perturbed trajectory x t (W + δW ) is approximated by (∂x t /∂W )δW , where ∂x t /∂W is a N×N 2 matrix. We call this matrix the sensitivity matrix ∆ t , and the time evolution of ∆ t is given by where ∂F /∂x is the Jacobian matrix and ∂F /∂W is the parametric Jacobian matrix. It should be noticed that the two trajectories x t (W ) and x t (W + δW ) coincide for t ≤ T ′ , and thus ∆ t = 0 for t ≤ T ′ . The growth rate of the displacement between x t (W + δW ) and x t (W ) with respect to the perturbation vector δW is obtained by ∆ t δŴ , where δŴ = δW /|δW |. The maximum value of |∆ t δŴ | at time step t is given by the maximum singular value σ t 1 of ∆ t matrix. σ t 1 can be obtained by performing the singular value decomposition of ∆ t . σ t 1 diverges for t → ∞ when the maximum Lyapunov exponent of the trajectory is positive. On the other hand, σ t 1 oscillates or converges to a constant value when the maximum Lyapunov exponent is negative. Note that no parameters except for W are perturbed in this study. Once a network W is given, ∆ t and its σ t 1 are estimated for each time step. We define parameter sensitivity γ as the logarithm of an average of σ t 1 along the trajectory: Here, we regard the first T ′ − 1 steps of the trajectory as transient, and discard them. We use T ′ = 1000 and T max = 1500.
Our goal is to sample networks with small γ. However, networks with positive λ 1 tend to have large γ. For this reason, and because almost all networks sampled by random sampling should have positive λ 1 , random sampling is not suitable for the sampling of networks with small γ. Thus, we again apply multicanonical Monte Carlo and perform random walks in λ 1 space with the same w(λ 1 ) estimated above. These random walks facilitate efficient sampling of the networks with small γ, because networks with negative λ 1 are efficiently sampled. We also obtain the two-dimensional density D(λ 1 , γ) of λ 1 and γ as follows: we construct the two-dimensional histogram h(λ 1 , γ) through the random walks, and, after h(λ 1 , γ) is constructed, D(λ 1 , γ) is calculated by Figure 3 shows D(λ 1 , γ) for N = 10 and 20 with k = 2 and 5. These results show that although networks with small γ can take various values for λ 1 , the vast majority of such networks with small γ have negative but near zero λ 1 (see red rectangles in Fig. 3). Figure 5 shows an example of an optimized network with small γ. Although we have examined network topology of optimized networks sampled by multicanonical Monte Carlo, no characteristic difference was found in topology between networks with small γ and those with large γ.
This indicates that, when robustness against network perturbations is optimized (i.e., when γ is minimized), networks with negative but near zero λ 1 will appear with high probability. This appearance of systems with marginal stability can be interpreted as self-organization of the edge of chaos. In this scenario, the system automatically comes close to the edge of chaos without tuning parameters.
In order to confirm that optimization of obustness against network perturbations leads to the emergence of the edge of chaos, we perform simulated annealing. Here, we define a parameter sensitivity Γ without the linear approximation: where δW represents an average over realization of δW . We minimize this Γ using the simulated annealing. The average is taken over 100 samples of δW , and the parameters T ′ = 1000 and T max = 1500 are used. In each step of the simulated annealing, a transition from the current state W to a proposed candidate W ′ is accepted if and only if the ratio exp[−β(Γ(W ′ ) − Γ(W ))] is smaller than a random number uniformly distributed in (0, 1]. Here, temperature 1/β is lowered with the progress of simulation step n. We choose this β as β = 10n/3 (for n < 30000) and β = ∞ (for n ≥ 30000). Note that the function Γ that we aim to minimize fluctuates due to the finite sample size of δW , and thus occasionally an inferior network could be accepted or a suitable network could be rejected, even for β = ∞. In Fig. 4, we plot λ 1 and Γ for networks that are sampled during the simulated annealing. These results indicate that most networks obtained in the last half of the simulations (n ≥ 30000) are in the region −0.05 ≤ λ 1 < 0, and we regard this region as the edge of chaos.

Discussion
In summary, using multicanonical Monte Carlo method, we have observed emergence of the systems at the edge of chaos as a self-organization phenomenon with only the requirement of robustness against network perturbations, which can be interpreted as mutational robustness in the context of the gene regulatory network. We have also performed simulated annealing and confirmed this scenario. We emphasize that no fine tuning of other parameters, such as number of input links k or model parameters (a, ǫ), is needed. The emergence of the edge of chaos with the requirement of mutational robustness is somehow counterintuitive because robustness against network perturbations (mutational robustness) seems to be positively correlated with dynamical stability. The mechanism of the emergence of the edge of chaos, revealed by a multicanonical Monte Carlo method, is as follows: when mutational robustness is required, selected systems need to have λ 1 ≤ 0 because γ for λ 1 > 0 diverges as t → ∞. The density D(λ 1 ) is an increasing function for λ 1 ≤ 0. Therefore, the density of networks becomes largest at λ 1 ∼ 0 under the condition λ 1 ≤ 0 (see Fig. 1 and 3, red rectangles in Fig. 3 indicate the degeneracy of a large numbers of networks with small γ). Due to this degeneracy, systems have a high probability of being at the edge of chaos.
Similar results have been found in recent numerical studies of gene regulatory network [18,19,20], indicating that systems that have the ability to reach a stable fixed point with transient chaotic behavior appear with only the requirement of robustness against genotypic perturbations. These results can be also interpreted as the emergence of the edge of chaos. However, it has not until now been explained why such systems are selected. In this paper, we have proposed a mechanism for the emergence of the edge of chaos, namely that the vast majority of networks that are robust against network perturbations have marginal dynamical stability, and thus networks at the edge of chaos are selected when robustness against network perturbations is required. Note that the converse is not necessarily true: networks with marginal dynamical stability are not always robust against network perturbations. It is reasonable to think that this degeneracy of marginally stable networks appears whenever parameters are set in the chaotic phase, in which chaotic systems are obtained under random construction of systems for N → ∞. Based on the fact that similar results were found in the previous studies [18,19,20,21], most of what we discussed here seems not to depend on the details of the specific model.
Mutational robustness seems to be a natural requirement for living systems. The present study provides the following possible scenario for the emergence of the edge of chaos: the requirement of mutational robustness drives organisms to the edge of chaos, whether or not staying in such a regime is preferable for living systems. In fact, several recent studies have suggested that gene networks of real organisms stay at the edge of chaos [22,23,24,25,26]. We also expect that the scenario is also applicable for the explanation of criticality in neural dynamics [27,28,29], because synaptic connections seem to be designed so that neuron firing patterns do not change radically due to perturbations in the connection. The concept presented here that robustness against network perturbations leads to the edge of chaos may also provide insights into the design of artificial robust networks. λ a < λ 1 < λ b . One step of the Metropolis -Hastings algorithm using w(λ 1 ) is performed as described in the Appendix (i), and then, w(λ 1 ) is modified as w(λ 1 (W new )) → w(λ 1 (W new ))/f , where f is the modification factor. By this modification, w(λ 1 ) of frequently visited bins are made smaller, and thus the appearance of such λ 1 is suppressed. The modification factor is set as f = e at the beginning of the simulation, and the factor is gradually reduced through the simulation to approach to unity in the manner described in [12,13]. We accumulate the number of samples in λ 1 ∼ λ 1 + δλ 1 to make a histogram h(λ 1 ). The Metropolis -Hastings step and the modification of w(λ 1 ) are repeated until the histogram h(λ 1 ) is sufficiently flat in λ a < λ 1 < λ b . This procedure allows us to construct w(λ 1 ) that is inversely proportional to D(λ 1 ). Once such w(λ 1 ) is obtained, w(λ 1 ) is fixed. The Metropolis-Hastings algorithm using this w(λ 1 ) enables a uniform sampling in λ 1 space (i.e., a random walk in λ 1 space) by generating a Markov chain, because each step is biased proportional to D(λ 1 ) −1 and probability of the candidate λ 1 is proportional to D(λ 1 ). Note that in this study we determine λ a and λ b so that P (λ 1 < 0) is precisely calculated; we adopt the criteria that P (λ a )/P (λ 1 = 0) ∼ 10 −2 and P (λ b )/P (λ * ) ∼ 10 −2 , where λ * is the value of λ 1 for the peak of D(λ * ) (see Fig. 1).  Network size dependence of the probability P (λ 1 < 0) of networks with negative λ 1 . The logarithms of P (λ 1 < 0) for k = 2 − 5 decrease linearly or slightly faster than linear as functions of N. The logarithm of P (λ 1 ) for k = N − 1 decreases quadratically.