PHYSICAL REVIEW E 90, 022814 (2014)

Stability of Boolean networks: The joint effects of topology and update rules Shane Squires,1,2 Andrew Pomerance,3 Michelle Girvan,1,2,4 and Edward Ott1,2,5 1 Department of Physics, University of Maryland, College Park, Maryland, USA Institute for Research in Electronics and Applied Physics, University of Maryland, College Park, Maryland, USA 3 Potomac Research, LLC, Alexandria, Virginia, USA 4 Institute for the Physical Sciences and Technology, University of Maryland, College Park, Maryland, USA 5 Department of Electrical Engineering, University of Maryland, College Park, Maryland, USA (Received 4 April 2014; revised manuscript received 18 July 2014; published 25 August 2014; publisher error corrected 22 October 2014) 2

We study the stability of orbits in large Boolean networks. We treat the case in which the network has a given complex topology, and we do not assume a specific form for the update rules, which may be correlated with local topological properties of the network. While recent past work has addressed the separate effects of complex network topology and certain classes of update rules on stability, only crude results exist about how these effects interact. We present a widely applicable solution to this problem. Numerical simulations confirm our theory and show that local correlations between topology and update rules can have profound effects on the qualitative behavior of these systems. DOI: 10.1103/PhysRevE.90.022814

PACS number(s): 89.75.−k, 05.45.−a, 64.60.aq, 87.16.Yc

I. INTRODUCTION

Systems of interconnected Boolean elements have been successfully used to model the macroscopic behavior of a wide variety of complex systems. Examples include genetic control [1], neural networks [2], ferromagnetism [3], infectious disease spread [4], opinion dynamics [5], and other applications in economics and geoscience [6]. Each of these diverse models share the same basic structure: a set of nodes, each of which has a binary state (0 or 1) at a given integer time t, and a set of update rules that determines the state of each node at time t + 1 given the states of the nodes at time t. The relationships between nodes define a graph, where an edge is drawn from node j to node i if the update rule for node i depends on the state of node j . Depending on the desired application, the model’s graph can be random [1], a lattice [3], fully connected [2], or have other complex topology [7–9]. The states of nodes can be updated deterministically or stochastically, synchronously or asynchronously. Finally, in many cases, update rules are considered to be randomly generated, and, in such cases, they can be drawn from many different ensembles [7,9–12]. One important question about a Boolean network is whether or not it is stable, i.e., whether or not small perturbations of a typical initial state tend to grow or shrink as the system evolves. This question may have important ramifications for systems biology and neuroscience: it has been hypothesized that both gene networks [13] and neuronal networks [14] exist near the critical border separating the stable and unstable regimes. Recently, Pomerance et al. [8] introduced the additional hypothesis that orbital stability of the gene regulatory system may be causally related to cancer. Specifically, they suggested that mutations that promote instability may be a contributing factor for some types of cancers. This hypothesis is consistent with recent results which indicate that a distinguishing feature of cancer cells is extreme variation in gene expression levels [15]. In this paper, we present and numerically verify a general method for studying the stability of large, directed Boolean networks with locally treelike topology. Here, by a locally 1539-3755/2014/90(2)/022814(10)

treelike network we mean that, if two nodes j and i are connected by a short directed path from j to i, it is unlikely that there is another short directed path from j to i [16]. (Such networks are called locally treelike because the set of short paths originating at a given node forms a tree.) This property allows us to make the approximation that the states of two inputs to a node are uncorrelated. Even in cases where the network contains substantial clustering, analyses based on this approximation have been found to yield accurate results [8,17], while making an analytic treatment of the system tractable. Our results offer a means of assessing the stability of a variety of Boolean network systems for which, up to now, no generally effective method has been available. In particular, our approach can be applied to specific complex networks with given topology and essentially arbitrary update rules. Because of this generality, our technique may be useful for modeling systems whose topology and update rules are determined experimentally. For example, one question that may be important in gene regulatory systems is how silencing (or “knocking down”) a particular gene will affect the expression of other genes, as well as the stability of the dynamics of the entire system. Our analysis provides a framework for addressing this kind of question. This is in contrast to previous work, which typically obtained analytical results on the stability of ensembles of networks with given properties, e.g., in-degree distributions or joint in-/out-degree distributions. Our framework can also be used in a more general sense to study how topology, update rules, and correlations between the two affect stability. We demonstrate this with two examples illustrating that the joint effects of network topology and update rules can have profound effects on Boolean network dynamics that cannot be captured by previous theories. The paper is organized as follows. We first introduce the Boolean network model and previous work on stability using the so-called annealed and semiannealed approximations (Sec. II). Next, we show how to apply the semiannealed approximation to a broad class of update rules, and we use this to derive a set of “dynamical biases” that describe the

022814-1

©2014 American Physical Society

SQUIRES, POMERANCE, GIRVAN, AND OTT

PHYSICAL REVIEW E 90, 022814 (2014)

steady-state dynamics of the system (Sec. III). These are then used to analyze the stability of the system (Sec. IV), and we demonstrate agreement of our predictions with numerical simulations for two illustrative examples (Sec. V). The first of these examples is primarily pedagogical, but the second is application oriented, using threshold update rules. II. BACKGROUND A. Boolean networks and stability

Boolean networks are discrete-state dynamical systems in which each of the N nodes of a network has a binary state xi (t) = 0 or 1 at each integer-valued time t, and is updated at the next time t + 1 to a new binary state xi (t + 1) that is determined from the time t states of its network inputs. For now we assume that updates are synchronous, but in Appendix A, we demonstrate that the stability criterion that we obtain is the same whether nodes are updated synchronously or asynchronously. Consider a node i that has Ki network inputs, j1 . . . jKi . The new state of node i is determined by a binary-valued update rule Fi , according to   xi (t + 1) = Fi xj1 (t), . . . ,xjKi (t) . (1) Each Fi may be specified in the form of a “truth table” listing all the 2Ki possible inputs and the corresponding output. Denoting the vector of input states to node i at time t as Xi (t) = (xj1 (t), . . . ,xjKi (t)), we have xi (t + 1) = Fi (Xi (t)).

(2)

The network structure is represented using an adjacency matrix A, where Aij = 1 if there is an edge from j to i (that is, if Fi depends on xj ), and Aij = 0 otherwise. The question we address is whether the dynamics resulting from Eq. (2) are stable with respect to small perturbations. ˜ = We consider two orbits, x(t) = (x1 (t), . . . ,xN (t))T and x(t) (x˜1 (t), . . . ,x˜N (t))T . We define the normalized Hamming distance between these two orbits as the fraction of the nodal states that differ between them, ˜ H (x(t), x(t)) =

N 1  |xi (t) − x˜i (t)|. N i=1

(3)

We will say that node i is “damaged” at time t if xi (t) = x˜i (t), and therefore the normalized Hamming distance may be interpreted as the fraction of nodes that are damaged. To define stability, we assume that N  1 and that the orbits ˜ 0 ))  1. If are close at some initial time t0 , i.e., H (x(t0 ), x(t ˜ H (x(t), x(t)) decreases to zero as the dynamics evolve under Eq. (2), we consider the system to be stable; if it grows to O(1), we consider the system unstable. Our main theoretical result is a criterion for stability that accounts for the joint effects of network topology (i.e., the Aij ) and node dynamics (i.e., the functions Fi ). B. Annealed and semiannealed approximations

The stability of Boolean networks was addressed in the original work of Kauffman [1,13], where Boolean networks were first proposed as a model for genetic dynamics. Since

then, stability has been a principal focus of work on Boolean networks, and the range of analytically tractable cases has expanded dramatically through the introduction of new techniques [8,10,18]. Kauffman assumed a so-called N -K network topology in which Ki was the same at each node, Ki = K, and the K inputs to each node were chosen randomly from amongst the (N − 1) other nodes. Further, for each of the 2K input states Xi , Fi (Xi ) was chosen to be 0 or 1 with probability 1/2. Derrida and Pomeau later generalized this model by introducing a bias parameter 0 < b < 1 such that, for each given input Xi , Fi (Xi ) = 1 with probability b [10]. We refer to such functions as “biased update rules.” Numerically, it was found that such systems undergo a phase transition from a stable phase to an unstable phase as the parameters K and b are varied. Derrida and Pomeau also proposed a method of stability analysis that predicts the critical point for this transition by introducing a related “annealed” system. First, note that the problem that they and Kauffman were interested in was one in which the network (Aij ) and the node dynamics (Fi ) are initially randomly chosen, fixed forever after, and then used to create the dynamics (“quenched” randomness). The annealed problem is different in that the random choices of the network edges and nodal update rules are made at every time step (but note that the same sequence ˜ In contrast with the of update rules is used for both x and x). stability of the quenched system, the stability of the annealed system can be analytically determined; specifically, the system will be stable if 2b(1 − b)K < 1 [10]. Derrida and Pomeau conjectured that, in the large N limit, the stability boundaries for the quenched and annealed situations are approximately the same. Here, we refer to this conjecture as the “annealed approximation.” The annealed approximation has been very well supported by the results of numerical experiments, and it has gradually been extended to systems other than N -K networks with biased update rules. In particular, the annealing approach can be applied to systems with a distribution of in-degrees [19–23], joint in-degree/out-degree distributions [24,25], “canalizing” [26] update rules [7,27–31], and threshold update rules [11,32–35]. While the annealed approximation has provided powerful techniques for studying ensembles of random networks, it cannot be applied to specific networks with complex topological features. Reference [8], on the other hand, developed a “semiannealing” technique that attempts to answer the question of stability based directly on the adjacency matrix of the system under examination. In a semiannealed system, the network topology is fixed and may be complex; only the update rules are selected randomly at each time step. The semiannealed approximation is the hypothesis that the stability of a semiannealed system, which can be computed analytically, is a good approximation for the stability of a typical system with quenched update rules on the same network. Besides allowing for the consideration of complex network topologies, this approach also makes it possible to consider the contributions of individual nodes to the system dynamics. The semiannealed approximation has been used to analyze the stability of complex locally treelike networks with biased update rules, generalized so that the bias bi may vary from node to node [8,36], and has also been applied to canalizing update

022814-2

STABILITY OF BOOLEAN NETWORKS: THE JOINT . . .

PHYSICAL REVIEW E 90, 022814 (2014)

rules [9,36] and multilayer Boolean networks [37]. In each case, the stability of the system is determined by the magnitude of the largest eigenvalue of a modified adjacency matrix, which allows for the consideration of quite general topological features. In fact, Refs. [8,9] considered systems with a variety of complex features, such as degree assortativity, community structure, motifs, and correlation between node degree Ki and node bias bi . Many of these features dramatically affect the stability of the system in ways that cannot be predicted using fully annealed techniques. Numerical results strongly support the validity of the semiannealed analysis [8,9,36,37]. Thus, semiannealed techniques have vastly extended the range of network topologies whose stability can be analyzed. However, they have only been applied to two limited classes of update rules, and a clear next step is the development of a semiannealed analysis that can be applied with fewer constraints on the update rules. This is a difficult problem because, for general update rules, performing a stability analysis requires detailed knowledge of the steady-state behavior of the system. Here, we address this problem by introducing a set of “dynamical biases,” pi , which describes the steady-state behavior of the semiannealed system and can be determined from a set of consistency conditions.

III. GENERALIZED SEMIANNEALING AND DYNAMICAL BIASES

We begin by describing the dynamics of our general semiannealed model. The update rule Fi for each node i will be randomly selected at each time step according to the following procedure. Each node i is assigned an update rule ensemble, Ti , which contains both a set of update rules as well as the probabilities of those rules being selected. Here, we will denote a particular update rule as f and the associated selection probability as Pr[Fi = f ]. For example, consider the case where, for some node i, Ti contains two rules, f1 and f2 , with associated probabilities Pr[Fi = f1 ] = 34 and Pr[Fi = f2 ] = 14 . In this case, the state of i will be updated according to xi (t + 1) = f1 (Xi (t)) with probability 34 , or according to xi (t + 1) = f2 (Xi (t)) with probability 14 . On each time step, a similar random selection is made independently at each network node, using the update rule ensemble for that node. In the problem of stability, the update rule is selected only once for each node at each time ˜ step, and the resulting update rule is used for both x and x. When using the semiannealed approximation to understand the stability of a quenched system, one must choose an appropriate update rule ensemble based on the system under study. In typical applications of Boolean networks, each quenched update rule Fi is initially selected from a random ensemble. In these cases, one should use the same ensemble for Ti , so the choice is straightforward. Otherwise, the selection of the update rule ensemble Ti may depend on the case of interest. We illustrate this with examples in Sec. V. The dynamics of a semiannealed system are somewhat similar to those of probabilistic Boolean networks [38]. They may be described by the probabilities F¯i (Xi ) that the state of node i, given inputs Xi , will be 1 on the next time step. This

is just the average of Fi (Xi ) over the ensemble Ti ,  F¯i (Xi ) = Pr[Fi = f ]f (Xi ),

(4)

f ∈Ti

where we have used the fact that f (Xi ) = 0 or 1. It is important to note that F¯i (Xi ) is solely determined from Ti , i.e., independent of the update rule assignments at other nodes. Thus, computation of F¯i (Xi ) is straightforward. Another important quantity, which is also purely local and therefore can also be calculated directly from Ti , is the probability that a node i will become damaged if its input states on the previous time step are given by Xi and X˜ i for the original and perturbed orbits. We denote this quantity by Di (Xi ,X˜ i ), and it may be calculated using Di (Xi ,X˜ i ) = Pr[xi (t + 1) = x˜i (t + 1)|Xi ,X˜ i ]  = Pr[Fi = f ]|f (Xi ) − f (X˜ i )|,

(5a) (5b)

f ∈Ti

where, once again, we have used the fact that f (Xi ) = 0 or 1. This formalism allows for the treatment of a wide variety of update rules, including but not limited to previously studied cases. As an example, consider the case of biased update rules with node-dependent biases bi (i.e., for any given combined input state, the output is assigned to be 1 with probability bi ). Semiannealing over biased rules can be represented as selecting from the ensemble of all possible update rules with appropriate probabilities. Specifically, the probabilities Pr[Fi = f ] follow the distribution bin (1 − bi )m−n , where m = 2Ki is the number of possible input states, and n is the number of input states Xi that result in f (Xi ) = 1. It is easy to show that this distribution results in F¯i (Xi ) = bi for all Xi , as well as the well-known result that the “sensitivity” of a biased function is Di (Xi ,X˜ i ) = 2bi (1 − bi ) for all X˜ i = Xi [8,10,27]. Alternatively, here (as in most cases), it is not actually necessary to specify the ensemble Ti . In fact, our subsequent analysis only makes use of the functions F¯i and Di , and in many cases, including the example in Sec. V B, it is more convenient to specify them directly than to calculate the probabilities Pr[Fi = f ]. Note that this approach can also be used to describe a quenched system if only one function is in each Ti . However, we will primarily focus on cases in which stochasticity plays a significant role in the dynamics. Typically, this results in the existence of a single ergodic attractor for the system, simplifying the analysis. (A sufficient, but certainly not necessary, condition for this to occur is if 0 < F¯i (Xi ) < 1 for all i and Xi .) The semiannealed approximation asserts that the stability of this attractor in the semiannealed system is the same as that of typical attractors in quenched systems. We will comment on the application of our results to quenched dynamics below. Now we introduce the dynamical biases pi and a set of consistency equations that can be used to determine them. Given an orbit on the attractor of the semiannealed system, we define pi to be the fraction of time that the state xi (t) of node i is 1, or the probability that xi (t) = 1 at a randomly chosen time t. A global average of this quantity has been considered in fully annealed network ensembles (e.g., [29,30,32]).

022814-3

SQUIRES, POMERANCE, GIRVAN, AND OTT

PHYSICAL REVIEW E 90, 022814 (2014)

We first note that pi is determined by the set of probabilities Pr[Xi ] of i receiving each possible input vector Xi , using  Pr[Xi ] Pr[xi = 1|Xi ] pi = Xi

=



Pr[Xi ]F¯i (Xi ),

(6)

Xi

where F¯i is as defined in Eq. (4). Assuming that the network topology is locally treelike, the states of the inputs to node i can be treated as statistically independent [8,17]. Therefore, the probabilities Pr[Xi ] are determined by the biases of i  s inputs. Letting Ji denote the set of indices of nodes that are inputs to i,  [xj pj + (1 − xj )(1 − pj )], (7) Pr[Xi ] = j ∈Ji

where we have used the fact that each xj = 0 or 1. That is, the probability of a state Xi in which an input xj = 0 is weighted by a factor of (1 − pj ), and the probability of a state in which xj = 1 is weighted by a factor of pj . Inserting Eq. (7) into Eq. (6) yields a set of N equations for the N dynamical biases pi . In typical applications, these equations can be solved by choosing a random initial value for each pi and iterating Eqs. (6) and (7) until they converge. Moreover, we find numerically that there is a unique solution that is stable under iteration, and this represents the attractor of the semiannealed system. We will assume this in the analysis below, but here we note that the analysis can be extended in a natural way to cases where this does not occur. For example, some systems (including completely quenched systems) may have more than one attractor, in which case the attractors can be found by iterating multiple initial guesses for pi , and the stability of each attractor can be evaluated individually [39]. IV. STABILITY ANALYSIS

We now consider the stability of the semiannealed system. We define a vector y(t) such that yi (t) is the probability that i is damaged at time t, i.e., yi (t) = Pr[xi (t) = x˜i (t)].

(8)

Marginalizing over Xi and X˜ i in Eq. (8) and inserting Eq. (5a),  Pr[Xi (t),X˜ i (t)]Di (Xi ,X˜ i ). (9) yi (t + 1) = Xi

X˜ i

Because we are considering the question of stability, we ˜ are close to each other in the have assumed that x(t) and x(t) sense of Hamming distance for times t close to t0 . In this case, we can drop terms of O(y 2 ) on the right-hand side of Eq. (9), which corresponds to ignoring the possibility that Xi (t) and X˜ i (t) differ for two or more inputs. Moreover, if X˜ i (t) and Xi (t) are the same, Di = 0 via Eq. (5b), so nothing is contributed to the sum in Eq. (9). Therefore, the only values of X˜ i that contribute substantially to the sum are ones in which X˜ i (t) and j Xi (t) differ for exactly one node j . Hence we define Xi (t) to be a vector that is the same as Xi (t) except that the state of input node j is damaged [x˜j (t) = xj (t)]. Using this notation,

we can rewrite Eq. (9) as     j j Pr Xi (t),Xi (t) Di Xi ,Xi . yi (t + 1) =

(10)

j ∈Ji Xi

Furthermore, because the network is locally treelike, we treat the inputs to node i as uncorrelated, yielding   j Pr Xi (t),Xi (t) = Pr[Xi ] Pr[xj (t) = x˜j (t)] (11) = Pr[Xi ]yj (t). When substituted into Eq. (10), this leads to    j yi (t + 1) = yj (t) Pr[Xi ]Di Xi ,Xi . j ∈Ji

(12)

Xi

Since the second sum is time independent, we can write  Rij yj (t), (13a) yi (t + 1) = j

Rij ≡



 j Pr[Xi ]Di Xi ,Xi ,

(13b)

Xi

where Rij = 0 when there is no edge from j to i. (The secondorder terms in this derivation are discussed further in Appendix B, where we derive an expression for the critical slope of the stability phase transition.) Rij may be interpreted as the probability that damage will spread from node j to node i; in analogy with the terminology of Ref. [27], we call it the “effective activity” of j on i. The average of the normalized Hamming distance over all possible perturbations and realizations of the semiannealed  ˜ dynamics is H (x(t), x(t)) = N1 i yi (t), so the stability of the system is determined by whether or not the elements of y(t) grow with time. This can be understood by writing Eq. (13b) in matrix form, y(t + 1) = R y(t).

(14)

The growth of y is determined by whether or not R has any eigenvalues λ whose magnitude is greater than 1. If so, and the initial perturbation has a nonzero component along the left eigenvector associated with one of these eigenvalues, the expected Hamming distance will initially grow as λt . Therefore, we may use the largest-magnitude eigenvalue of R, which we denote λR , to predict the stability of the system. It can be shown that this eigenvalue is always real and positive. With additional assumptions that are applicable to commonly studied cases, it can also be shown that this eigenvalue has multiplicity 1 and is well separated in magnitude from other eigenvalues, which implies that it can be computed efficiently through power iteration of R on an initial random vector [40]. Therefore, we can classify the stability of the system using λR , according to the rule λR < 1 stable, λR = 1 critical, λR > 1 unstable.

(15)

One major advantage of our approach is that, from a computational perspective, evaluating p, R, and λR is typically much faster than finding the average Hamming distance through simulations. We discuss this further in Appendix C, along with other computational aspects of the above solution.

022814-4

STABILITY OF BOOLEAN NETWORKS: THE JOINT . . .

PHYSICAL REVIEW E 90, 022814 (2014)

Another potential advantage of Eq. (15) is that it may facilitate qualitative understanding of the effects of certain network or update rule features on stability. For example, previous work has shown that introducing assortativity to a network tends to increase the largest eigenvalue of certain types of modified adjacency matrices, and this effect can be estimated analytically [41,42]. This suggests that similar results may be derived for λR in future work. It is also possible that empirical gene regulatory networks will contain important topological or update rule features that may not correspond to commonly studied network metrics. We emphasize that our analysis nonetheless applies to these cases. Moreover, such features may be highlighted by any cases where the direct evaluation of λR gives a substantially different stability prediction than annealed or approximation-based techniques. In these cases, research aimed at understanding the discrepancy may lead to the discovery or study of new network features. Finally, in contrast to fully annealed techniques, in which nodes are distinguished only by their degree signature, semiannealing also allows for the consideration of individual nodes and their contribution to instability in the full system. This enables the treatment of questions such as the optimal method of stabilizing or destabilizing a network through some limited set of possible external interventions. One interesting case is where the intervention strategy is to permanently silence or activate a small set of nodes [i.e., forcing xi (t) = 0 or xi (t) = 1 for some node i]. Another case is when damage can actively be corrected during the dynamics for some nodes. In both cases, the problem of selecting nodes for intervention can be addressed using our framework. This and other related questions are discussed further in Appendix D.

0.4 (a)

largest random smallest

H

0.3 0.2 0.1 0 0 2.5

0.1 (b)

0.3

0.4

0.5

0.3

0.4

0.5

1.5 1 0.5 0

H

0.2

0.1 (c)

0.2

α

largest random smallest

0.1

0 1

V. NUMERICAL RESULTS

We now use the general framework presented above to analyze two cases that illustrate the effects of correlations between local topological features and update rules. In each example, we construct a single network with N = 105 nodes using the configuration model [43]. Because experimental results suggest that gene regulatory networks have short-tailed in-degree distributions and long-tailed out-degree distributions [44,45], we use Poisson-distributed in-degrees and scale-free out-degrees. We chose a mean degree of 4 for the Poisson distribution, which corresponds to an exponent γ ≈ 2.2 for the scale-free distribution. In Figs. 1 and 2, we plot the average Hamming distance H and λR against a tuning parameter for each model. When calculating H , we use a set of quenched update rules, whereas we calculate λR using the semiannealed theory described above. In the figures, we show H for both a single quenched system as well as an average over 50 sets of quenched update rules. For each quenched set of update rules, we calculate H through the following procedure. We generate an initial condition for an orbit x(t), then time-evolve it until a time t0 = 100, where we expect it to have completed any transient behavior. Here, the initial condition is chosen randomly from a uniform distribution [xi (0) = 0 or 1 with equal probability]. ˜ Next, we apply a small perturbation to create the orbit x(t). Specifically, we randomly choose a small fraction ε = 0.01 of

α

largest random smallest

2 λR

0.2

λR

1.5

FIG. 1. (Color online) Normalized average Hamming distance H and λR for a network with XOR, OR, and AND update rules (described in Sec. V A). Filled markers are averaged over 50 quenched realizations of the update rules, while hollow markers show a single quenched realization. Red squares, blue circles, and green triangles represent different correlations between network topology and update rules; see text for details. (a) When the fraction of XOR rules α is used as a tuning parameter, the stability transitions for the three cases are far apart. Locations where λR = 1 are marked by arrows and agree with the observed transitions. (b) λR plotted as a function of α. (c) Plotting H against λR directly, we see that the transition for all three curves occurs at λR = 1.

the components of x(t0 ) and damage their states [i.e., if node i is chosen to be damaged, then x˜i (t0 ) = 1 − xi (t0 )]. Finally, we time-evolve both trajectories in parallel, and we measure ˜ the long-time behavior of H by averaging H (x(t), x(t)) from t = 400 to t = 500. Finally, we repeat this procedure for ten initial conditions x(0), and average H over initial conditions to obtain H . A. Example 1: XOR, OR, and AND update rules

In this example, we illustrate the effect of correlations by assigning a node either a highly “sensitive” update rule (XOR)

022814-5

SQUIRES, POMERANCE, GIRVAN, AND OTT

PHYSICAL REVIEW E 90, 022814 (2014)

update rules, and the remaining nodes are evenly split between OR and AND rules. We consider three cases: (1) XOR is assigned to the αN nodes with the largest in-degree (red or medium gray markers in Fig. 1); (2) XOR is randomly assigned to αN nodes irrespective of their degrees (blue or dark gray markers); or (3) XOR is assigned to the αN nodes with smallest in-degree (green or light gray markers). In all three cases, the remainder are randomly assigned OR or AND. In numerical simulations, all update rule assignments are quenched. In order to find appropriate semiannealing probabilities to use in our theoretical prediction, we note that the initial assignment of XOR is deterministic in cases (1) and (3), but OR and AND are assigned randomly. Therefore, in the theory, we treat the network and the identity of the XOR nodes as fixed, and anneal over the OR and AND nodes, assigning a probability of 12 to choosing either OR or AND on each time step. Other annealing choices are also possible, but we choose this because it most straightforwardly resembles the quenched assignment of update rules above. As can be seen in Fig. 1(a), the values of α at which the three cases become unstable are rather different, thus demonstrating that stability is strongly affected by correlation between the local topological property of nodal in-degree and the sensitive XOR update rule. This agrees with the intuitive expectation that assigning update rules that result in high values of Rij to topologically important nodes will increase instability. In Fig. 1(b), we see that our stability criterion correctly accounts for this behavior. When we replot H against λR , we see that in each case the network becomes unstable at λR ≈ 1, as predicted by the theory. This is also pointedly illustrated by the vertical arrows in Fig. 1(a) marking the values of α at which λR = 1. XOR

H

0.2

(a)

corr. uncorr. anti.

0.1

0 0

0.5 – θ

1

λR

(b) 1 correlated uncorrelated anticorrelated 0.5 0

0.5 – θ

1

0.2

H

(c)

correlated uncorrelated anticorrelated

0.1

0 1

λR

1.2

B. Example 2: Threshold networks

FIG. 2. (Color online) Normalized average Hamming distance H and λR for a threshold network (described in Sec. V B). The notation and results are analogous to Fig. 1, except that threshold rules are used and θ¯ is used as a tuning parameter in (a) and (b). Once again, the theoretical predictions [arrows in panels a and c] agree with the observed transitions.

or a less sensitive update rule (OR or AND), based on the node’s in-degree. That is, the update rule at each node i is randomly drawn from three classes: (a) XOR, whose output is 1 (0) if i has an odd (even) number of inputs that are 1; (b) OR, whose output is 1 if and only if at least one input is 1; or (c) AND, whose output is 1 if and only if all Ki inputs are 1. Following [27], we refer to XOR as highly sensitive because any single damaged input will cause its output to become damaged, so Rij = 1 whenever there is an edge from j to i. On the other hand, if node i has OR or AND as its update rule, damaging node j will damage node i only if every node other than j is 0 or 1, respectively. Thus in cases (b) and (c), Rij depends on the node biases pj , which are obtained by solving Eqs. (6) and (7). For OR, Rij = k (1 − pk ), and for AND, Rij = k pk , where the products are taken over all inputs k that are not equal to j . Figure 1 shows results for a network with these three classes of update rules. We assign a fraction of nodes α to have

We now consider networks with threshold rules, which take the form 

 (16) xi (t + 1) = U wij xj (t) − θi , j

where U denotes the unit step function, θi is a threshold value, and wij is a signed weight whose magnitude reflects the strength of the influence of node j on node i and whose sign indicates whether node j “activates” or “represses” node i (i.e., promotes xi to be 1 or 0). Such threshold rules may be recast as Boolean functions Fi by enumerating all possible Xi and calculating whether the weighted sum of inputs exceeds the threshold θi . Conversely, threshold rules are appropriate for Boolean network applications in which each edge has a fixed “activating” or “repressing” character. The special case of random networks with θi = 0 and wij = ±1 has been considered previously using the annealed approximation [11,32–35,46]. Results for threshold networks are shown in Fig. 2 and are generated as follows. To assign the weight wij for each edge j → i, we first choose a sign, positive or negative, each with probability 12 , to indicate whether the edge will be activating or repressing. Then, the weight is drawn from a normal distribution with mean ±1 (according to the chosen sign) and standard deviation 14 .

022814-6

STABILITY OF BOOLEAN NETWORKS: THE JOINT . . .

PHYSICAL REVIEW E 90, 022814 (2014)

We also consider two additional cases where the weights are either correlated or anticorrelated to a topological property of the network, the product of a node’s in-degree and outdegree, which we call the “degree product.” It has been shown that nodes with high degree product play a crucial role in the stability of Boolean networks [8,25,42]. Here, we generate the correlated and anticorrelated cases by exchanging weights between pairs of edges in the original (“uncorrelated”) assignment of weights. Specifically, we repeat the following procedure. First, we select two random edges j1 → i1 and j2 → i2 in the network. Next, we identify the edge for which i has a higher degree product. Finally, in the correlated (anticorrelated) case, we exchange the values of the two weights if doing so would increase (decrease) the weight of the edge with the higher degree product. We repeat this procedure E/2 times, where E is the number of edges in the network, so that each edge is expected to be considered for one exchange. In this example, we model the case in which the thresholds of different nodes are similar, but not necessarily equal. In the theory, we treat this case by annealing the thresholds θi over a normal distribution with a mean θ¯ and standard deviation 1 σθ = 10 . By Eqs. (16), (4), and (5b), ⎡ ⎛ ⎞⎤  1 F¯i (Xi ) =  ⎣ ⎝ wij xj − θ¯ ⎠⎦ , (17a) σθ j    j  j Di Xi ,Xi = F¯i (Xi ) − F¯i Xi ,

or pi ≈ 1. When the average threshold θ¯ is sufficiently large, most nodes are frozen in the 0 state,  either because they have no input combinations for which j wij xj > θi , or because some inputs are also frozen, so such input combinations do not occur. In this situation, damage is typically only propagated by edges with wij > θi , so most edges have Rij ≈ 0, and the ¯ system  is stable. As θ is decreased, it becomes more common for j wij xj to be close to θi , which causes most nodes to become unfrozen and damage to spread more easily. Finally, at low θ¯ , it is possible for loops or small groups of nodes with mutually activating edges to become frozen in the 1 state, and for such nodes to repress other nodes, which become frozen in the 0 state. Here, the pattern of freezing (for example, whether high-degree or low-degree nodes are more likely to be frozen in the 1 state) depends sensitively on the correlations between degree and edge weight, leading to the differences between the three cases above. The results of Fig. 2 again illustrate the importance of the fact that the stability criterion, Eq. (15), can accommodate correlations between network topology and update rules, in contrast to fully annealed techniques. Note that all three cases in Fig. 2(a) share the same network topology and thresholds, as well as the distribution of weights. Therefore, any theory that does not take into account the joint effects of network topology and update rules would be unable to predict the marked difference in the behavior of these three cases.

(17b)

where (x) is the cumulative normal distribution, x (2π )−1/2 −∞ exp(−t 2 /2)dt. Equation (17) can be used to find pi , Rij , and λR . Here, as in many cases, it is not necessary to enumerate the ensemble of update rules or the associated probabilities, because F¯i and Di can be calculated directly. When finding the average normalized Hamming distance H in the numerical simulations, we choose a quenched value for θi by writing θi = θ¯ + δθi , where δθi is drawn from a normal distribution with mean 0 and standard deviation σθ . In Fig. 2, we show results for both a single quenched set of δθi (hollow markers) as well as an average over 50 quenched sets of δθi (filled markers). In each case, single quenched realizations show similar behavior to the average, in agreement with the semiannealing hypothesis. More striking is the qualitative difference between the anticorrelated case and the correlated and uncorrelated cases. At low thresholds, the anticorrelated network is stable, whereas both the correlated and uncorrelated case are unstable. As the threshold is increased, the anticorrelated network becomes unstable before becoming stable again at large thresholds. This agrees with Fig. 2(b), where we see that in all three cases, λR initially increases with increasing threshold, but the anticorrelated case is the only one where λR < 1 initially. Finally, in Fig. 2(c), we replot the same data for H against λR . We see that in all three cases the stability transition clearly occurs at λR = 1, confirming our analysis. Although the full dynamics of this system can be quite complex, the behavior in Fig. 2 can be understood qualitatively by considering the number of nodes whose states are nearly “frozen” [47,48] (i.e., nodes that rarely change their state after the transient period). In our theory, this corresponds to nodes for which the solution to Eqs. (6), (7), and (17a) has pi ≈ 0

VI. DISCUSSION

We have presented a general framework for addressing the question of orbit stability of large, locally treelike Boolean networks, given arbitrary network topology and update rules. There are four steps in this process: (1) identify the update rule ensemble Ti (or the functions F¯i and Di ) for each node i; (2) calculate the dynamical biases pi of each node i by iterating Eqs. (6) and (7); (3) calculate the effective activities Rij using Eq. (13b); and (4) find the largest eigenvalue of R, λR , which determines the stability of the system through Eq. (15). As illustrated above, the first step requires a judicious selection of which aspects of the update rules should remain quenched, but is typically straightforward thereafter. The second step—and, in general, the use of dynamical biases in the analysis—is essential for all but the simplest classes of update rules, because the probability of damage spreading from one of a node’s inputs depends on the states of the other inputs. For this reason, the analysis presented here is the first that can be applied to specific complex Boolean networks with general update rules, taking into account the steady-state dynamics of the system. As examples of the application of our general stability criterion, we analyzed both a pedagogical case and the case of threshold networks, where all update rules are assumed to be of the form of Eq. (16). Whereas research into the stability of Boolean networks has primarily focused on either topology or update rules alone, our results confirm that correlations between the two can have profound qualitative effects on the dynamical properties of a Boolean network. The dynamical biases pi and the stability indicator λR offer systematic measures for evaluating these effects. As discussed in Sec. IV, our approach also makes it possible to study the effects of interventions that modify either the

022814-7

SQUIRES, POMERANCE, GIRVAN, AND OTT

PHYSICAL REVIEW E 90, 022814 (2014)

topology or update rules of a given Boolean network, because we can treat specific networks with quite general update rules. Further consequences and development of these results remain to be explored. Because the activity matrix R depends on both the dynamical biases p and the damage probability functions Di , introducing any feature to the network or update rules can affect λR in several (sometimes conflicting) ways, as in the case of threshold networks. Further exploration of other features is likely to yield other interesting examples. Another interesting direction of work is to develop techniques similar to those used in Refs. [8,41,42], which would allow pi and λR to be estimated analytically in particular cases without solving Eqs. (6) and (7) numerically. This would allow greater intuitive understanding of the effects of various characteristics of the network or update rules, as well as the interaction between these characteristics. We hope that our results here will draw more attention to these issues, and, more broadly, to the fascinating interplay between topology and update rules in the dynamics of Boolean networks. ACKNOWLEDGMENTS

This work was supported by ARO grant W911NF-12-10101.

where ρ is a diagonal matrix with ρi in each row, I is the identity matrix, and R is the activity matrix. In order to see that Eq. (15) also applies in this case, we note that, at criticality, y(t + 1) = y(t), so that Eq. (A1) reduces to ρ(R − I) y(t) = 0. This has a solution for y = 0 only if λR = 1. Note, however, that in this case, for λR > 1, the growth rate of the Hamming distance will be at a rate of the order of 1/N smaller than the rate of the synchronously updated networks, because N time steps of asynchronous update correspond to one time step of synchronous update. APPENDIX B: CRITICAL SLOPE

The second-order terms in the expansion in Eq. (13a) may be used to derive the critical slope of H near λR = 1. We include a sketch of this derivation because it may be useful for near-critical approximations or for designing networks with extreme behavior near the critical point. We begin by considering input combinations in which two distinct inputs j,k j and k are damaged, which we denote X˜ = Xi in analogy j with our definition of Xi above. Following similar steps as those that led to Eq. (13a) and retaining all terms of O(y 2 ), we obtain   Rij yj + Rij k yj yk , yi = j

APPENDIX A: ASYNCHRONOUS UPDATES

Rij k

Asynchronous updates may arise in discrete state systems for several reasons. Links may have nonuniform delays, δij , that model delays arising from, for example, the chemical kinetics of gene regulation. In this case, the dynamics would be described by a modified version of Eq. (1) in which the state of node i at time t depends on the states of its inputs j at times (t − δij ). Another alternative is a model in which nodes are individually chosen to be updated in a stochastically determined order. Here, we show that the stability condition given in the main text applies not only to the case of synchronous nodal updates, but to asynchronous models as well, including both of these examples. In particular, we consider update times, τ1 < τ2 < · · · < τt < · · · , where the update intervals, (τt+1 − τt ), are arbitrary, incommensurate, and do not influence the analysis. Since the update times are incommensurate, we approximate the deterministic choice of node to update at each time, indexed by integer t, with a stochastic process where node i (and only node i) is chosen to be updated by Eq. (2) with probability ρi . This is also appropriate, of course, for systems that are inherently stochastic. To analyze this case, we must adjust the approximate update equation, Eq. (14). Since node i is chosen independently of the values of the nodes, the joint probability at time step t that node i is chosen for update and that node i differs between the two initial  conditions after the update is given approximately by ρi Rij yj (t). If node i is not chosen for update at this time step, yi does not change. Putting this together, we get, for small t and small initial perturbations, yi (t + 1) ≈ ρi Rij yj (t) + (1 − ρi )yi (t), which we rewrite in matrix form as y(t + 1) = ρ(R − I) y(t) + y(t),

(A1)

j,k

 1 j,k  ≡ Pr[Xi ]Di Xi ,Xi − Rij , 2 X

(B1)

i

where Rij is defined as in the main text, and we have assumed that y has reached its steady state behavior. Note that when j = k, Rij k = 0. Now we may derive the critical slope through a perturbation expansion near the critical point, yi = εHyi1 + ε2 yi2 and λR = 1 + ελ1R , where superscripts for yi and λR refer to the level of the perturbation expansion. From Eq. (14), y1 must be equal to the right Frobenius-Perron eigenvector of R, denoted  v here, which we normalize so that i vi = 1. Inserting the second-order expansion and simplifying, we obtain   Rij yj2 + H 2 Rij k vj vk . (B2) yi2 = H λ1R vi + j

j,k

This expression may be simplified by using the left FrobeniusPerron eigenvector, denoted u. Multiplying through by ui and summing over i, the left-hand side and the second term on the right-hand side of Eq. (B2) cancel to leading order in ε. With the remaining terms, we find that the critical slope mc = H /λ1R is given by  i ui vi . (B3) mc = −  i,j,k Rij k ui vj vk This result may be used numerically to find the critical slope in particular cases. It may also be used to approximate the critical slope analytically, when good approximations for u and v are known, as in Refs. [8,41,42]. APPENDIX C: COMPUTATIONAL CONSIDERATIONS

Because of our use of dynamical biases and the locally treelike approximation, the stability analysis presented above

022814-8

STABILITY OF BOOLEAN NETWORKS: THE JOINT . . .

PHYSICAL REVIEW E 90, 022814 (2014)

is applicable to very large values of N for sparse networks. This offers a tremendous computational improvement over previous theoretical treatments of similar systems; for example, the analysis of probabilistic Boolean networks in Ref. [38] relies upon a state transition matrix of size 2N × 2N , which is intractable for networks with more than a few dozen nodes. We also find that our method is much more efficient than direct computation of the Hamming distance, which requires O(nic t3 E) operations, where nic is the number of initial conditions used, t3 is as defined in Sec. II A, and E is the number of edges. This is because nic and t3 must be relatively large in order to obtain good estimates for H . The most computationally intensive step in our analysis is iterating Eqs. (6) and (7), which requires on the order of 2Ki steps for each i (see below). Calculating Rij for each edge also requires 2Ki steps, but only needs to be performed once. The eigenvalue λR may be found through power iteration, which requires only O(E) operations per iteration (using sparse matrix-vector multiplication) and also converges quickly. Iterating Eqs. (6) and (7) is fast for large-N networks as long as most nodes have Ki < 15, because few iterations are needed for convergence. In cases where the average degree is large, direct application of our method becomes intractable, but in many cases additional simplifying assumptions can be used for such systems. For example, for threshold networks with high in-degree, rather than enumerating all 2Ki possible input states, one could use  the central limit theorem to approximate the distribution of j wij xj . In the case of moderate values of Ki , it is worth noting that, while a na¨ıve approach would use Ki 2Ki steps to calculate the set of Pr[Xi ] for node i, this can be achieved somewhat more efficiently through the following process when Ki > 2. Begin with a list containing the values (1 − pj ) and pj for some input j . Then for each additional input k, copy the list, multiply one copy by (1 − pk ) and the other by pk , then join the two copies. This way, all Pr[Xi ] are generated in approximately 2Ki +1 operations. APPENDIX D: INDIVIDUAL CONTRIBUTIONS TO INSTABILITY

A great deal of information about the effects of network interventions (as discussed in Sec. IV) can be derived from the left and right eigenvectors associated to λR , which may be

[1] [2] [3] [4] [5]

S. A. Kauffman, J. Theor. Biol. 22, 437 (1969). J. J. Hopfield, Proc. Natl. Acad. Sci. USA 79, 2554 (1982). R. J. Glauber, J. Math. Phys. 4, 294 (1963). M. E. J. Newman, Phys. Rev. E 66, 016128 (2002). C. Castellano, S. Fortunato, and V. Loreto, Rev. Mod. Phys. 81, 591 (2009). [6] B. Coluzzi, M. Ghil, S. Hallegatte, and G. Weisbuch, Int. J. Bifurcat. Chaos 21, 3511 (2011). [7] S. Kauffman, C. Peterson, B. Samuelsson, and C. Troein, Proc. Natl. Acad. Sci. USA 100, 14796 (2003). [8] A. Pomerance, E. Ott, M. Girvan, and W. Losert, Proc. Natl. Acad. Sci. USA 106, 8209 (2009).

obtained through power iteration along with λR . We denote these eigenvectors u and v, respectively. For example, to first order, it is only possible to cause Y > 0 by initially damaging at least one node with ui = 0, and only nodes with vi = 0 will have yi (t) > 0 when t is large. These results can also be interpreted using the relationship between the stability of Boolean networks and percolation explored in Ref. [36]: nodes with ui = 0 and vi = 0, respectively, form the giant in-component and out-component of the network of edges for which Rij = 0. It is also possible to use u and v to estimate the contribution of an individual node to damage spreading in the network, using the concept of dynamical importance introduced in Ref. [49]. Let λR denote the change in λR that would result if a node i were unable to spread damage to its outputs. Then

λR /λR ≈ ui vi /uT v. A similar result for edge removal is also available [49]. These results are useful for understanding the flow of damage through a specific network, as well as for stabilization in systems where there is a mechanism for correcting the states of damaged nodes. Another approach to intervening in a network, more applicable to gene regulatory applications such as those discussed in Sec. I, is to silence nodes or edges. In the theory, silencing a node would correspond to forcing pk = 0 for some k, and silencing an edge would correspond to considering only the values of Xi in which xj = 0 in the sums in Eqs. (6), (7), and (13b). The change in p after such an intervention describes its effect on the behavior of the rest of the system (i.e., for gene regulation, the change in the expression levels of other genes), and the resulting change in λR describes the overall effect on the stability of the system. Given the goal of stabilizing an unstable network through a small intervention, one could compare a number of possible interventions and choose the one that results in the smallest value of λR . When doing so, it is best to repeat the analysis of Sec. IV for each modified system, rather than applying a perturbative approximation, because the change in p that results from silencing a node can affect R and λR in complicated ways. Even so, a large number of such interventions can be compared in a computationally efficient manner by storing the values of p and v for the original system. These can then be used as initial conditions for the iterative processes used to find the dynamical biases and λR , decreasing the convergence time for both processes.

[9] A. Pomerance, M. Girvan, and E. Ott, Phys. Rev. E 85, 046106 (2012). [10] B. Derrida and Y. Pomeau, Europhys. Lett. 1, 45 (1986). [11] T. Rohlf and S. Bornholdt, Physica A 310, 245 (2002). [12] F. Li, T. Long, Y. Lu, Q. Ouyang, and C. Tang, Proc. Natl. Acad. Sci. USA 101, 4781 (2004). [13] S. Kauffman, The Origins of Order: Self-Organization and Selection in Evolution (Oxford University Press, Oxford, 1993). [14] W. L. Shew, H. Yang, T. Petermann, R. Roy, and D. Plenz, J. Neurosci. 29, 15595 (2009). [15] H. C. Bravo, V. Pihur, M. McCall, R. A. Irizarry, and J. T. Leek, BMC Bioinf. 13, 272 (2012).

022814-9

SQUIRES, POMERANCE, GIRVAN, AND OTT

PHYSICAL REVIEW E 90, 022814 (2014)

[16] In models with synchronous dynamics, we require an even less restrictive condition, that it is unlikely that there are two short paths of the same length. [17] S. Melnik, A. Hackett, M. A. Porter, P. J. Mucha, and J. P. Gleeson, Phys. Rev. E 83, 036112 (2011). [18] A. Mozeika and D. Saad, Phys. Rev. Lett. 106, 214101 (2011). [19] R. V. Sol´e and B. Luque, Phys. Lett. A 196, 331 (1994). [20] B. Luque and R. V. Sol´e, Phys. Rev. E 55, 257 (1997). [21] J. J. Fox and C. C. Hill, Chaos 11, 809 (2001). [22] M. Aldana and P. Cluzel, Proc. Natl. Acad. Sci. USA 100, 8710 (2003). [23] M. Aldana, Physica D 185, 45 (2003). [24] D.-S. Lee and H. Rieger, J. Theor. Biol. 248, 618 (2007). [25] D.-S. Lee and H. Rieger, J. Phys. A 41, 415001 (2008). [26] In a canalizing rule, there is one input that, when it is in a particular state, determines the output completely, regardless of the states of the other inputs. [27] I. Shmulevich and S. A. Kauffman, Phys. Rev. Lett. 93, 048701 (2004). [28] S. Kauffman, C. Peterson, B. Samuelsson, and C. Troein, Proc. Natl. Acad. Sci. USA 101, 17102 (2004). [29] A. A. Moreira, L. A. Nunes Amaral, Phys. Rev. Lett. 94, 218702 (2005). [30] J. Kesseli, P. R¨am¨o, and O. Yli-Harja, Phys. Rev. E 74, 046104 (2006). [31] T. P. Peixoto, Eur. Phys. J. B 78, 187 (2010). [32] F. Greil and B. Drossel, Eur. Phys. J. B 57, 109 (2007). [33] A. Szejka, T. Mihaljev, and B. Drossel, New J. Phys. 10, 063009 (2008). [34] J. G. T. Za˜nudo, M. Aldana, and G. Mart´ınez-Mekler, in Information Processing and Biological Systems, edited by S. Niiranen and A. Ribeiro (Springer-Verlag, Berlin, 2011), pp. 113–151. [35] R.-S. Wang and R. Albert, Phys. Rev. E 87, 012810 (2013). [36] S. Squires, E. Ott, and M. Girvan, Phys. Rev. Lett. 109, 085701 (2012).

[37] E. Cozzo, A. Arenas, and Y. Moreno, Phys. Rev. E 86, 036115 (2012). [38] I. Shmulevich, E. R. Dougherty, S. Kim, and W. Zhang, Bioinformatics 18, 261 (2002). [39] We do not consider attractors that correspond to periodic or chaotic orbits in the map defined by Eqs. (6) and (7), although they have been shown to exist in unusual cases [30,50], because they are unlikely to correspond to biologically realistic dynamics. [40] These statements follow directly from the Frobenius-Perron theorem when R is a primitive matrix. They are also nearly always true for networks that are above the percolation threshold without strong community structure. In this case, it is safe to assume that the network has only a single macroscopic strongly connected component, whose associated FrobeniusPerron eigenvalue dominates all other eigenvalues. The effects of community structure on the spectra of adjacency matrices are discussed in detail in Ref. [51]. [41] J. G. Restrepo, E. Ott, and B. R. Hunt, Phys. Rev. E 76, 056119 (2007). [42] E. Ott and A. Pomerance, Phys. Rev. E 79, 056111 (2009). [43] M. E. J. Newman, SIAM Rev. 45, 167 (2003). [44] N. Guelzim, S. Bottani, P. Bourgine, and F. K´ep`es, Nat. Genet. 31, 60 (2002). [45] R. Dobrin, Q. K. Beg, A.-L. Barab´asi, and Z. N. Oltvai, BMC Bioinf. 5, 10 (2004). [46] K. E. Kurten, J. Phys. A 21, L615 (1988). [47] H. Flyvbjerg, J. Phys. A 21, L955 (1988). [48] T. Mihaljev and B. Drossel, Phys. Rev. E 74, 046101 (2006). [49] J. G. Restrepo, E. Ott, and B. R. Hunt, Phys. Rev. Lett. 97, 094102 (2006). [50] M. Andrecut and M. K. Ali, Int. J. Mod. Phys. B 15, 17 (2001). [51] S. Chauhan, M. Girvan, and E. Ott, Phys. Rev. E 80, 056114 (2009).

022814-10