Finite size effects and symmetry breaking in the evolution of networks of competing Boolean nodes
Abstract
The effects of the finite size of the network on the evolutionary dynamics of a Boolean network are analyzed. In the model considered, Boolean networks evolve via a competition between nodes that punishes those in the majority. Previous studies of the model have found that large networks evolve to a statistical steady state that is both critical and highly canalized, and that the evolution of canalization, which is a form of robustness found in genetic regulatory networks, is associated with a particular symmetry of the evolutionary dynamics. Here it is found that finite size networks evolve in a fundamentally different way than infinitely large networks do. The symmetry of the evolutionary dynamics of infinitely large networks that selects for canalizing Boolean functions is broken in the evolutionary dynamics of finite size networks. In finite size networks there is an additional selection for input inverting Boolean functions that output a value opposite to the majority of input values. These results are revealed through an empirical study of the model that calculates the frequency of occurrence of the different possible Boolean functions. Classes of functions are found to occur with the same frequency. Those classes depend on the symmetry of the evolutionary dynamics and correspond to orbits of the relevant symmetry group. The empirical results match analytic results, determined by utilizing Pólya’s theorem, for the number of orbits expected in both finite size and infinitely large networks. The reason for the symmetry breaking in the evolutionary dynamics is found to be due to the need for nodes in finite size networks to behave differently in order to cooperate so that the system collectively performs as well as possible. The results suggest that both finite size effects and symmetry are important for understanding the evolution of realworld complex networks, including genetic regulatory networks.
pacs:
87.23.Kg, 05.65.+b, 87.14.Gg, 89.75.HcI Introduction
Boolean networks Kauffman (1969, 1993); Aldana et al. (2003); Drossel (2008) have been studied extensively over the past three decades. They consist of a directed graph in which the nodes have binary output states that are determined by Boolean functions of the states of the nodes connected to them with directed inlinks. They have applications as models of gene regulatory networks as well as models of physical, social, and economic systems. As “coarsegrained” models of genetic networks they aim to capture the essential features of the dynamical behavior of the real networks while simplifying local gene expression to a binary (on/off) state Bornholdt (2005). Recent work has demonstrated that, despite their simplicity, Boolean networks can indeed describe many of the important features of the dynamics of biological genetic circuits Albert and Othmer (2003); Li et al. (2004); Shmulevich et al. (2005); Lau et al. (2007). For example, it has been shown that a Boolean network model of the segment polarity gene regulatory networks that control embryonic segmentation in Drosophila melanogaster can reproduce the wildtype gene expression pattern and ectopic patterns due to various mutants Albert and Othmer (2003).
Motivated by the fact that evolution plays a crucial role in forming the regulatory relations among genes, a number of models that evolve the structure and dynamics of Boolean networks have been studied Kauffman and Smith (1986); Bornholdt and Sneppen (1998); Paczuski et al. (2000); Bornholdt and Rohlf (2000); Bornholdt and Sneppen (2000); Luque et al. (2001); Bassler et al. (2004); Bassler and Liu (2005); Liu and Bassler (2006); Szejka and Drossel (2007); Rohlf (2007); Braunewell and Bornholdt (2007). These evolutionary Boolean network (EBN) models generally seek to determine the properties of the networks that result from the evolutionary mechanism being considered. For example, many of the studies have focused on the topology of links of the networks that result from evolutionary mechanisms that rewire the links. Other studies have focused on finding evolutionary mechanisms that result in networks that have dynamics that are robust against various types of perturbations, or that result in networks that are in a “critical” state poised between ordered and “chaotic” dynamical behavior.
An example of an EBN is the model of competing Boolean nodes first introduced in Ref. Paczuski et al. (2000), and studied subsequently in Refs. Bassler et al. (2004); Bassler and Liu (2005). In this model the nodes of a Boolean network compete with each other in a variant of the Minority game Challet and Zhang (1997). The network evolves by changing the Boolean function of the node that loses the game to a new randomly chosen Boolean function. In the principal variant of the model that has been studied, only the Boolean functions of the nodes evolve, not the links. Although, it should be noted that changing the Boolean functions used by the nodes effectively results in a rewiring of the links Reichhardt and Bassler (2007). In the original paper on the model, it was shown that the network selforganizes to a statistically steady, nontrivial critical state with this evolutionary mechanism. Later it was discovered that the critical state the network evolves to is highly canalized Bassler et al. (2004). Canalization Waddington (1942) is a type of network robustness known to exist in genetic regulatory networks Rutherford and Lindquist (1998); Quietsch et al. (2002); Bergman and Siegal (2003). It exists when certain expression states of a subset of the genes that regulate the expression of a gene control the expression of the gene regardless of whether or not the other genes that otherwise affect its expression are being expressed. In this case, the states of the subset of regulatory genes that control the expression of a gene are “canalizing” inputs to the gene. Canalization is thought to be an important property of developmental biological systems because it buffers their evolution, allowing greater underlying variation of the genome and its regulatory interactions before some deleterious variation can be expressed phenotypically Wagner (2005).
The canalized nature of the evolved steady state was demonstrated by preforming an ensemble of similar simulations Bassler et al. (2004). Each simulation in the ensemble involved a different random network and began with a different initial condition. After they were run long enough to allow the networks to evolve to a steady state, the average frequency that the various possible Boolean functions occurred was measured and averaged over the ensemble of simulations. For Boolean networks of size nodes in which each node has inlinks, it was found that the possible Boolean functions of three variables organize into 14 different classes in which all of the functions in each class occur with approximately the same frequency. It was then found that the various classes of functions mostly could be distinguished by the fraction of canalizing inputs that their functions have. Moreover, the classes whose functions have larger fractions of canalizing inputs, that is, the ones that are more canalizing, occurred with larger frequency.
The reason that there are 14 different classes of Boolean functions is due to the symmetry properties of the evolutionary dynamics Reichhardt and Bassler (2007). The set of Boolean functions of inputs maps onetoone onto the set of configurations of the dimensional Ising hypercube in which each of the vertices of the hypercube have a binary state. The different classes correspond to the group orbits of the “Zyklenzeiger” group Harrison (1963) which is the hypercubic, or hyperoctahedral symmetry group combined with parity. In this case, parity means simultaneously inverting the binary states of each of the vertices. A group orbit is the set of configurations that map into each other through applications of a group’s symmetry operations. The number of orbits a group has can be calculated using Pólya’s theorem, and using Pólya’s theorem it has been shown that there are 14 different group orbits of the 3 dimensional Zyklenzeiger group Reichhardt and Bassler (2007). Thus, the function classes correspond to orbits of the symmetry group of the evolutionary dynamics, and for large networks that group is the 3d Zyklenzeiger group.
In this paper, the effects of the finite size of the network on the selforganized evolution of networks of competing Boolean nodes are studied. Other EBN models Bornholdt and Sneppen (2000); Liu and Bassler (2006) have important finite size effects, indicating that small networks may behave in a fundamentally different way than large networks. Most previously reported results on this model have been for relatively large networks with nearly 1000 nodes. Here, we find that small networks evolve very differently than large networks. In particular, it will be shown that the competition of nodes in very small networks does not result in the evolution of a canalized network. Instead, the evolutionary mechanism preferentially selects for “input inverting” Boolean functions that have output states opposite to the majority of the inputs they receive. We will also see empirically that for finite size networks the Boolean functions have 46 different classes, instead of 14. Finally, it will be shown that this change from 14 to 46 classes is due to a type of symmetry breaking in the evolutionary dynamics of finite size networks. Pólya’s theorem will be used to show that, indeed, 46 classes are expected from the reduced symmetry of the finite size evolutionary dynamics.
The remainder of this paper is organized as follows. In Sec. II, the definition of random Boolean networks (RBNs) and of canalizing Boolean functions is given, a few important properties of RBNs are noted, and the algorithm used for the evolution of Boolean networks based on a competition between nodes is stated. Section III presents empirical results obtained from simulations of finite sized networks. In Sec. IV, Pólya’s theorem is used to calculate the number of function classes in finite size networks based on a group of symmetry operations that is reduced from that which is relevant to the evolution of infinitely large networks. Section V explains that the symmetry breaking occurs because the behavior of nodes in finite size networks must differ from those in infinitely large networks in order for them to cooperate so that they collectively perform as efficiently as possible. Finally, a summary of our results and a discussion of what they may imply about the importance of finite size and symmetry in the evolutionary dynamics of realworld complex networks is discussed in Sec. VI.
Ii The model
ii.1 Random Boolean Networks
In general, a RBN consists of nodes, , each of which has a dynamical Boolean state, or , that is a function of the Boolean states of other random nodes that regulate its behavior. The regulatory interactions are therefore described by a directed graph. As simple models of genetic networks, nodes can be interpreted as genes and their Boolean states indicate whether they are “not being expressed” or are “being expressed”. Different interpretations apply in other cases. For example, as agentbased models of financial markets, nodes can be interpreted as traders and their Boolean states indicate whether they are “buying” or “selling”. For synchronously updating RBNs, which are the only ones considered here, the state of each node at time is a function of all state of its regulatory nodes at time . Therefore, the discrete dynamics of the network is governed by
(1) 
where ,,, are those nodes that regulate node . The function is a Boolean function of inputs that determines the output of node i for all possible sets of input values. In this study, for simplicity, we consider only homogeneous RBNs in which each node has three input nodes, that is, for all . No selflinks, or multiple inlinks from the same node are allowed. Figure 1 illustrates an example of the type of homogeneous RBNs we consider.
Note that random Boolean functions can be generated with an “interaction bias” by setting the output value for each set of input to be ‘1’ with probability , and ‘0’ with probability . In modeling gene regulatory networks, the bias can be interpreted as a biochemical reaction parameter that defines the probability that a gene becomes “active” due to a particular set of regulatory inputs.
Given the Boolean state of each node at time , , the system state of network is defined as . The path that takes over time is a dynamical trajectory in the phase space of system. Note that there are different system states that are possible. For synchronous updating of the nodes, because the dynamics defined in Eq. 1 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
(2) 
for some time . The periodic part of the trajectory is the attractor of the dynamics, and the minimum that satisfies this equation is the period of the attractor.
The study of the dynamics of RBNs has a long history Derrida and Pomeau (1986); Derrida and Weisbuch (1986); Flyvbjerg (1988); Bastolla and Parisi (1997, 1998); Fox and Hill (2001); Aldana (2003); Aldana and Cluzel (2003); Socolar and Kauffman (2003); Samuelsson and Troein (2003); Kaufman et al. (2005); Greil and Drossel (2005); Drossel (2005); Drossel et al. (2005); Moreira and Amaral (2005). It is known that two distinct phases of dynamical behavior, “chaotic” and “ordered”, exist for RBNs that have Boolean functions that are randomly chosen depending on value of the interaction bias parameter Derrida and Pomeau (1986); Derrida and Weisbuch (1986); Flyvbjerg (1988). For such homogeneous RBNs with , the “chaotic” phase occurs for values of near 0.5, and the “ordered” phase occurs for values of near either 0 or 1. One way of distinguishing the two phases is to measure the distribution of its attractor periods beginning with random initial states Bastolla and Parisi (1997). For RBNs in the chaotic phase the distribution of attractor periods is sharply peaked near an average value that grows exponentially with system size , and for RBNs in the ordered phase the distribution of attractor periods is sharply peaked near an average value that is nearly independent of . There is a continuous transition between the two phases, the socalled “edge of chaos.” At the phase transition, RBNs are “critical” and have a broad, powerlaw distribution of attractor periods.
ii.2 Canalization and Symmetry
Canalization occurs in Boolean networks when the Boolean functions that define the dynamics of the nodes are canalyzing. A Boolean function is canalyzing if its output is fully determined by a specific value of one, or more, of its inputs, regardless of the value of the other inputs. Such an input value, or set of input values, that controls the output of the Boolean function is a canalizing input Kauffman (1969, 1993). How canalyzing a Boolean function of inputs is can be quantified by a set of numbers , which are defined the fraction of the different possible sets of input values that are canalyzing Kauffman (1993). Note that for Boolean functions that have inputs in total, there are different possible sets of input values. For example, in a Boolean function with inputs, there are 6 different possible single () input values, that is, two possible values for each of the three inputs.
An alternative way to understand canalization can be seen by mapping the Boolean functions onto configurations, or “colorings”, of the Ising hypercube Reichhardt and Bassler (2007). Each Boolean function of inputs maps onetoone onto a configuration of the dimensional Ising hypercube. The Ising hypercube is a hypercube which has each vertex labeled, or “colored”, either ‘0’ or ‘1’. In this representation, each of the possible sets of input values corresponds to a different vertex of the hypercube whose location in dimensional space is given by the Boolean vector constructed from its input values, and the coloring of the vertex corresponds to the output value of the function for that set of input values. An example of an Ising hypercube representation for a Boolean function is shown in Fig. 2, where the two possible colors of the vertices, white and black, correspond to the ‘0’ and ‘1’ output values respectively. Note that the vertices can be labeled either with the binary number constructed from the sequence of binary inputs, as in Fig. 2(a), or equivalently with the corresponding decimal number, as in Fig. 2(b). Note that there are 256 different Boolean functions, each corresponding to a different coloring of the Ising cube.
Note also that each of the different Boolean functions of inputs, or their equivalent colorings of dimensional Ising hypercubes, can be uniquely distinguished by an integer defined as
(3) 
where if vertex is white and if it is black. For instance, the Boolean function in Fig. 2 has .
The representation of Boolean functions as colorings of Ising hypercubes can greatly facilitate the recognition of canalizing functions Reichhardt and Bassler (2007). For a dimensional Ising hypercube, the fractions of canalizing inputs of a Boolean function are the fraction of its dimensional hypersurfaces that are homogeneously colored, that is, all vertices on the hypersurface have the same color. In general, there are such hypersurfaces possible. For example in Fig. 2, the whole cube is not homogeneously colored and so , it has one homogeneously colored face and so , and it has homogeneously colored edges and so .
With the Ising hypercube representation of Boolean functions it is easy to see that certain functions are related to others by symmetry. For example, if you permute the order of the labels of the inputs to a function, the resulting function and corresponding hypercube coloring will be related to the original by symmetry. Furthermore, since the ordering of the input labels is arbitrary the two symmetric functions should be indistinguishable in the evolutionary dynamics of the network. This observation suggests that the symmetry properties of Boolean functions is important in the dynamics EBNs. In particular, any symmetry in the dynamics for evolving the Boolean functions should be reflected in the symmetry of the evolved steady state, for instance in the symmetry in the frequency at which the different functions occur.
ii.3 Evolutionary Game
We consider the following variant of the Minority game for evolving the Boolean functions of a RBN Paczuski et al. (2000).

Begin with an unbiased RBN, that is, one whose functions are randomly chosen such that all of the different possible Boolean functions are equally likely to be chosen by each node, and with a random initial state.

For each update on the attractor, determine whether ‘0’ or ‘1’ is the output state of the majority of nodes, and give a point to each node that is part of the majority.

Determine which node is in the majority most often over the length of the attractor. That is, determine the node with the largest number of points. That node loses the game. If two, or more, nodes are tied with the largest score, pick one of them randomly to be the loser.

Replace the function of the losing node with a new randomly chosen unbiased Boolean function.

Return to step 2.
In practice, if the dynamical attractor is longer than some limiting time, , then score is kept only over that limited time. One iteration of the game is called an “epoch”. Exactly one evolutionary change is made to the network each epoch. Note that only the Boolean functions of the nodes evolve, the topology of the networks links does not change. The essential features of the game are (1) frustration Arthur (1994); Challet and Zhang (1997), since most nodes lose each time step, (2) negative reinforcement, since losing behavior is punished, and (3) extremal dynamics Bak and Sneppen (1993), since only the worst performing node’s Boolean function is changed. As mentioned above, previous studies involving large networks with nodes have found that this game causes the networks to evolve to a critical steady state that is highly canalized.
Iii Simulations and Results
Here we report results of simulations of ensembles of different size networks playing this game of competing Boolean nodes. Again, all of the networks simulated were homogeneous random networks with inputs per node. The sizes of the networks studied range from to . Note that is the smallest network with an odd number of nodes that can be formed without self or multiple connections. For each we found that the game caused the networks to evolve to a statistical steady state, and we measured the frequency that each of the 256 different functions occurred in the evolved steady state. This was accomplished by simulating, for each , an ensemble of independent realizations. Each realization started with an independent random network with random links and different unbiased random functions, and with an independent random initial state. The simulation of each realization was run for epochs to allow the network to evolve to the steady state. At the end of each simulation the functions used by each node were recorded, and then used to calculate the average frequency of each function for the ensemble of realizations.
We used the algorithm described in Ref. Liu and Bassler (2006) to find the dynamical attractor each epoch. In searching for the attractor, the maximum attractor period allowed, , varied with . For it was updates, and this limit was reached in about half the evolutionary game’s epochs. For all other , . As decreased, the period limit was reached in a decreasing fraction of the epochs, until only a small fraction, about 10%, reached the limit for . For the maximum possible length of an attractor is , since is the number of different system states. Therefore, for networks of this size, the limit was never reached. If is varied from these values, the quantitative results reported below do change, but the main qualitative findings concerning the symmetry breaking in the dynamics of finite size networks do not.
Figure 3 compares the results for large networks with nodes with those for small networks with nodes. Note that all simulations start with an unbiased RBN. Thus at the beginning of the simulations all Boolean functions are equally likely to be chosen by the nodes and therefore occur with a relative frequency of 1.0. However, in the steady state, clearly, some of the 256 possible functions occur more often than others. Furthermore, there are sets, or classes, of functions that occur with approximately the same relative frequency.
For the networks with nodes, as observed previously Bassler et al. (2004), there appear to be 14 different classes of functions. In order to make this clear, in the figure the functions numbers have been reordered so that all functions in the same class are grouped together. Furthermore, all functions in a class have the same canalization properties and the classes are sorted in descending order of how canalizing the functions in the class are. Note that ordered this way, the frequency of the functions is monotonically decreasing. Therefore, the more canalizing a function is the more likely it is to occur. The evolutionary process selects for canalization.
The results are quite different for networks with nodes. In this case, functions within the 14 classes found for large networks can now occur with very different frequencies. However, again there are classes of functions that occur with the same frequency. To see this note that the functions within each of the 14 classes are ordered according to the value of their frequency with the functions occurring with higher frequency first. Many of the 14 large network function classes split into multiple classes, and a total of 46 different function classes are observed for small networks. The complete reordering of the Boolean functions from their decimal numbers to the numbers used in Figs. 3 and 4 is given in Appendix A.
Note also for , that, in contrast to , there is not a monotonic decrease in frequency of the functions with increasing , and therefore the evolutionary process does not select for canalization. Instead, the evolutionary dynamics selects for “input inverting” functions that output the value opposite to the majority of their input values. How input inverting a Boolean function is can be quantified by its value of
(4) 
which is the weighted sum of the number of inverting outputs minus the weighted sum of the noninverting outputs, where the weight is the strength of the majority in the inputs. The larger is the more input inverting the function is. Appendix B lists the canalization and inversion properties of each class of functions. Every function in each of the 46 classes has the same value of .
The behavior of finite sized networks can be understood in more detail from the results shown in Fig. 4 which shows the relative frequency of functions in the evolved steady state for networks of size , 19, and 39. As network size decreases there is a crossover in the results of the evolutionary process from selecting for canalizing functions to selecting for input inverting functions. For all three sizes of networks, there are 46 classes of functions, and in each there is combination of preference for canalization and input inversion. For networks with nodes the most canalizing functions, and 1, and the least canalizing functions, and 255, occur most often and least often, respectively. However, for networks with nodes the most input inverting function, , and the least input inverting function, , occur most often and least often, respectively.
The splitting of the large network classes can perhaps best be seen by focusing on the behavior of the functions in the large network class , which are functions through 135 in the figure. The Ising cube representation of these 8 Boolean functions are shown in Fig. 5. Because for function , for functions through 131, for functions through 134, and for function , the 8 functions split into 4 different classes depending on the input inversion properties of the functions.
Although the preference for input inverting functions explains much of the frequency distribution in the evolved steady state, it should be noted that it does not explain all features of the behavior of small networks. For example, functions in class , and 11, occur slightly more often than those in class , through 17, despite the fact that those in are more input inverting. Also functions in the large network class , through 127, split into two classes, and , even though all functions in are explicitly not input inverting because their outputs are the same if their input values are simultaneously inverted. The Ising cube representation of these 8 Boolean functions are shown in Fig. 6.
The way that the splitting of large network classes quantitatively depends on network size is shown in the Fig. 7, where the difference in the relative frequency of functions and 135 is plotted as a function of . The inset shows the same data plotted on a loglog scale. The frequency difference appears to vanish either as a power law or slightly faster as increases. Analogous behavior occurs for the finite size splitting of the other function classes. Thus, for any finite value of , even for , there is some splitting and there are actually 46 function classes. Only in the limit of infinite does the splitting vanish and the number of function classes reduces to 14.
Iv Symmetry breaking
Many of the empirical simulation results presented in the previous section are caused by a breaking of symmetry in the evolutionary dynamics of finite size networks. As discussed earlier, each of the 256 Boolean functions of inputs maps onetoone onto a coloring of the 3d Ising hypercube, or, equivalently, the Ising cube. The symmetry of the evolutionary dynamics causes classes of functions to evolve “symmetrically” with the same frequency. Thus, the function classes correspond to orbits of the symmetry group of the evolutionary dynamics Reichhardt and Bassler (2007). In this case, a group orbit is the set of colorings that map into each other through applications of the group’s symmetry operations. The number of orbits a group has can be calculated using Pólya’s theorem Pólya (1937).
Pólya’s theorem states that the number of orbits a permutation group acting on the dimensional Ising hypercube has is
where is the number of operators , is the set of colorings, and is the set of colorings that are left invariant by . To evaluate , the size of , first express each operator in terms of its cycle structure, which is given by its cycle index , where . The cycle structure of describes how it permutes the vertices of the hypercube, and its cycle index indicates it has cycles of length 1, cycles of length 2, …, and cycles of length .
For example, consider a 2dimensional square. The four vertices can be labeled 0, 1, 2, and 3 moving clockwise. One permutation is induced by rotating the square clockwise through 180 degree. The cycle structure of this operator is , which means replace 0 by 2 and 2 by 0, and 1 by 3 and 3 by 1. The corresponding cycle index is .
For each operator without parity, the number of colorings left invariant, , is equal to , where . Operators with parity must be treated slightly differently. Parity means simultaneously reversing the binary state, or color, of each vertex. No colorings are left invariant by an operator with parity that has one or more cycles of length 1, and thus the cycle index of any such operator is defined to be zero. For each operator with parity, the number of colorings left invariant, , is equal to , where and is the Heaviside step function. Note that this means that some operations that include parity may leave no colorings invariant.
Alternatively, the number of orbits, , can be calculated by first finding the “cycle polynomial” of by adding the cycle indices for each together and dividing by . Then, because there are two possible colors of each vertex, replace each in the cycle polynomial with 2 and evaluate the expression. The result will be .
For infinitely large networks, the group describing the evolutionary dynamics is the 3d “Zyklenzeiger” group, which is the cubic group combined with parity Reichhardt and Bassler (2007). The canalization properties of any Boolean function are preserved by the operations of this group. This group has 96 operations. The operations include the 48 cubic symmetry operations plus each of those operations combined with parity. The cycle polynomial for this group is , and the number of group orbits is
Note that this is precisely how many function classes are found empirically in the large network limit.
However, for finite sized networks 46 classes are found empirically. What symmetry group describes the evolutionary dynamics of finite sized networks? To answer this question recall that, as shown in Fig. 4 and discussed in the previous section, the 46 classes appear from splitting many of the 14 classes that occur in the large network limit. This indicates that the symmetry group describing the dynamics of finite size networks is a subgroup of the one for infinitely large networks. Thus, the symmetry of the dynamics in the large network limit is reduced, or broken, in the dynamics of finite size networks.
Also recall from the empirical results that when one of the 14 large network classes splits, it splits into classes that often can be distinguished by the input inversion properties of the functions. This was illustrated in Fig. 5 which shows how the 8 functions in the large network class form four finite size network classes , , , and depending on their input inversion properties . While not all of the classes that split from the same large network class can be distinguished by the value of that their functions have, all the functions in each of the 46 finite size network classes have the same value of . Thus, the evolutionary dynamics of finite size networks is consistent with the symmetry of input inversion.
However, not all of the symmetry operations in the 3d Zyklenzeiger group that describes the evolutionary dynamics of infinitely large networks are consistent with input inversion. Only those that correspond to changes in the Ising cube that are symmetric about the (000)(111) diagonal axis of the cube are consistent with input inversion. Thus, of the 96 operations that describe the symmetry of the evolutionary dynamics of infinitely large networks, only 12 describe the symmetry of the evolutionary dynamics of finite size networks. Figure 8 illustrates 6 of these 12 operations by showing the vertex configuration that results from the operation. (a) shows the configuration that results from the identity operation. (b) and (c) show the configurations that result from rotating the identity cube 120 and 240 degree counterclockwise about the axis , respectively. (d) shows the configuration that results from reflecting the identity configuration about the surface containing vertices 0, 1, 6, and 7. Similarly, (e) and (f) show the configurations that result from reflecting the identity configuration about the surfaces containing vertices 0, 2, 5, and 7, and vertices 0, 3, 4, and 7, respectively. The configurations resulting from the other 6 operations, (g) through (l), can be obtained by inverting the 6 configurations shown in Fig. 8 about the center of the cube and applying the parity operation. For example, in (a), inversion means that exchanging vertex 0 with 7, vertex 3 with 4, vertex 1 with 6, and vertex 2 with 5.
Operation  Cycle Structure  Cycle Index 

(a)  
(b)  
(c)  
(d)  
(e)  
(f)  
(g)  
(h)  
(i)  
(j)  
(k)  
(l) 
These 12 operations form a group. This group is input inversion symmetric. The number of orbits of this group can also be calculated using Pólya’s theorem. Table I lists the cycle structure and cycle index for each of the operations in this group. Adding the cycle indices together, we find that the cycle polynomial for this group is , and that the number of group orbits is
This is exactly what is found empirically. Thus, the input inversion symmetry group describes the evolutionary dynamics of finite size networks. The precise form of the evolutionary dynamics is complicated
It is important to note that, although, as discussed earlier, not all of the classes that split from the same large network class can be distinguished by the value of that their functions have and that class splits even though all of the functions in the class are explicitly not inverting, the selection for input inverting functions by the evolutionary dynamics is still responsible for the splitting in these cases. Because the dynamics of finite size networks selects for input inverting functions, and because not all of the operations of the symmetry group for infinitely large networks are input inversion symmetric, the symmetry of the evolutionary dynamics of infinitely large networks is broken. Once this happens the number of group orbits, or, equivalently, function classes, increases. Since functions in different classes are not related by the symmetry of the evolutionary dynamics there is no reason that they should occur with the same frequency. It is complicated to calculate the quantitative values of the frequencies of the functions in the evolved steady state. They depend on more than just their fractions of canalizing inputs or their input inversion value . Their frequencies also depend on dynamical interactions of functions on different nodes. Therefore, in general, it should be expected that unless functions are in the same orbit of the symmetry group of the evolutionary dynamics that they not occur with the same frequency. However, if functions are in the same orbit, then they are constrained by the symmetry of the evolutionary dynamics to occur with the same frequency.
V Why does the symmetry break?
The remaining question to answer is why the evolutionary dynamics of finite size networks break the symmetry of the evolutionary dynamics of infinitely large networks. To answer this question note that the evolutionary game the nodes are playing is a variant of the minority game Challet and Zhang (1997); Paczuski et al. (2000) that effectively punishes nodes when their output is in the majority and rewards nodes whose output is in the minority. It has been shown that this effect of the evolutionary dynamics can result in the emergence of cooperation and organization of nodes throughout the system in an attempt to collectively perform as efficiently as possible Arthur (1994); Challet and Zhang (1997); Paczuski et al. (2000); Bassler et al. (2004); Lo et al. (2000); Hart et al. (2000); Anghel et al. (2004). The result being that at any time the number of nodes whose output is ‘0’ is almost equal to, or balanced by, the number of nodes whose output is ‘1’. This is true for any size network. However, the way that this balance is achieved in finite size networks is different from the way it is achieved in infinitely large networks.
The balance is easier to achieve in large networks than it is in small networks because the system has greater flexibility to balance the two states in order to achieve the global cooperation. In large networks every one of the different Boolean functions is likely to be chosen by some node in every epoch in each realization. Moreover, any local structure, or motif Milo et al. (2002); ShenOrr et al. (2002), of the network that occurs, including the Boolean functions of the associated nodes, is likely to be simultaneously occurring somewhere else in the network with opposite parity. Thus, the balance between nodes with output ‘0’ and nodes with output ‘1’ can be achieved simply by this “stochastic” balance that occurs naturally in large networks. For infinitely large networks this stochastic balance is exact.
For finite size networks, though, the balance can not be exactly achieved in this way. There will always be statistical fluctuations that prevent it from happening, and, for small networks, since each node chooses only one function in any epoch, there are not enough functions being used for the balance to be achieved stochastically. The lack of diversity in functions forces the evolutionary dynamics to place some extra selection pressure on these small number of chosen functions in order for the system to balance the output of the nodes. Because of this, it selects for input inverting functions.
To see why input inverting functions help balance the output of the nodes, consider the extreme case of . In many epochs the attractor that is found has period . That is, it is a fixed point, and all nodes always output the same Boolean value. Clearly, in this case, the best way for a function to work to achieve a collective balance in the network’s output is for it to output the value opposite to the majority of its input values, that is, for it to be input inverting. If a node’s function is not input inverting, then its output will always be in the majority and the node may lose the game and have its function changed. Input inverting functions are also favored in epochs with longer period attractors, but not as strongly. This mechanism for selecting input inverting functions also works in larger networks, but becomes decreasingly important as grows and the nodes become able to balance the collective behavior stochastically. This is why, as we found empirically, the symmetry breaking decreases as increases, and vanishes in the limit of infinitely large networks when exact collective balance can be achieved stochastically.
It is therefore in epochs where the attractor period is short that the selection for input inverting functions occurs. This contrasts with how canalization is selected for. It has previously been demonstrated that canalization is selected for by large networks in epochs in which attractors with long periods are found Bassler et al. (2004). Thus, in the evolved steady state of finite size networks, since a broad distribution of attractor periods are found, selection for both canalizing functions and input inverting functions occurs. However, as we also found empirically, as increases and the average length of an attractor period increases, the selection pressure for input inverting functions decreases while the selection pressure for canalizing functions increases.
Vi Conclusions
In this study of Boolean networks that evolve because of a competition between nodes for limited resources, we have found that the finite size of a network significantly affects the way that the network evolves and what the emergent properties of the system are. In small networks, the evolutionary dynamics selects for input inverting functions, but as the size of the network increases the evolutionary dynamics shifts from selecting for input inverting functions to selecting for canalizing functions until for infinitely large networks only canalization is selected for. This result can be understood in terms of the symmetry of the evolutionary dynamics. For networks in which each node receives 3 inputs, if the network is infinitely large, then the evolutionary dynamics has the symmetry of the 3d Zyklenzeiger group, which is the cubic group combined with parity, but, if the network is of finite size, then the evolutionary dynamics has the reduced symmetry of the subgroup consistent with input inversion. Thus, the symmetry of the evolutionary dynamics of infinite size networks is broken in the evolutionary dynamics of finite size networks. The symmetry of the evolutionary dynamics was revealed by the distribution of the frequency of occurrence of the various Boolean functions, which shows that there are classes of functions that all occur with the same frequency, and thus evolve symmetrically. Fourteen function classes are found for infinitely large networks, while 46 are found for finite size networks. The function classes correspond to orbits of the symmetry group of the evolutionary dynamics, and the number of orbits for both of the symmetry groups were calculated using Pólya’s theorem and found to match the above numbers. The reason for the symmetry breaking is that the way the nodes must behave to collectively perform as well as possible requires that the evolutionary dynamics be different for infinitely large and finite size networks, and this difference requires that they select for different properties in the Boolean functions.
Our previous work on the coevolution of topology and dynamics of Boolean networks similarly found that finite size can affect the evolution of a network. In that study, which considered a different type of evolutionary dynamics, we found that finite size networks evolve to have a broadly distributed indegree connectivity, whereas infinitely large networks evolve to have a homogeneous indegree connectivity Liu and Bassler (2006). Bornholdt and Rohlf also found important finite size effects in another evolutionary Boolean network model Bornholdt and Rohlf (2000).
As mentioned in the introduction, Boolean networks have been used to model a wide range of biological, physical, social, and economic systems, including both gene regulatory networks and financial markets. Many of these systems evolve through some evolutionary dynamics, and can be modeled by an EBN. Although we have focused on just one particular EBN model in this study, our results show that symmetry and finite size can have important effects on the evolutionary process. The properties of the evolved state can indicate much about the evolutionary process that leds to it. Real networks, even if they can’t be modeled by Boolean networks, are, of course, always of finite size, and symmetries exist in many realworld evolutionary processes. Therefore, our results suggest that finite size effects and symmetry may be important in the evolution of the complex networks describing many realworld systems.
Acknowledgements.
We thank Bogdan Danila for his careful reading of and comments on drafts of this paper. This work was supported by the NSF through grant No. DMR0427538.Appendix A The 46 Function Classes of Boolean functions
The classification of the 256 Boolean functions with inputs is shown in the following table. The name of each class consists of two parts. The first part is letters, such as A and Ca, that indicate the 14 classes found in the evolution of infinitely large networks. The classification scheme was also used in Ref. Bassler et al. (2004). The second part is numbers that indicate the multiple classes that one of the 14 classes split into during the evolution of finite size networks. The total number of finite size network classes is 46. For instance, class for infinitely large networks splits into 6 classes ( through ). The functions are distinguished by both their decimal numbers that uniquely determine the corresponding coloring of the Ising cube, and by their reordered numbers that sequentially places functions that are in the same class.
Class  Function  Function 
A  0, 255  
B1  1, 127  
B2  223, 4, 16, 247, 2, 191  
B3  254,128  
B4  239, 251, 32, 8, 253, 64  
Ca1  5, 119, 95,3, 63,17  
Ca2  48, 80, 10, 245, 68, 34, 221, 207, 187, 12, 175, 243  
Ca3  252, 238, 250, 192, 160, 136  
Cb1  15, 85, 51  
Cb2  204, 240, 170  
D1  87, 55, 21, 19, 7, 31  
D2  117, 69, 47, 81, 35, 13, 79, 59, 93, 11, 49, 115  
D3  143, 14, 179, 50, 213, 84  
D4  241, 205, 171, 112, 76, 42  
D5  138, 140, 186, 174, 220, 206, 208, 244, 162, 176, 196, 242  
D6  248, 168, 224, 234, 236, 200  
E1  20, 215, 18, 183, 159, 6  
E2  65, 111, 9, 123, 33, 125  
E3  130, 144, 246, 132, 190, 222  
E4  96, 72, 235, 40, 237, 249  
Fa1  129, 126  
Fa2  219, 36, 66, 189, 231, 24  
Fb1  23  
Fb2  43, 77, 113  
Fb3  142, 178, 212  
Fb4  232  
Fc1  27, 83, 39, 29, 53, 71  
Fc2  58, 116, 197, 177, 46, 92, 209, 114, 141, 139, 163, 78  
Fc3  184, 172, 202, 216, 228, 226  
G1  133, 62, 131, 118, 94, 145  
G2  103, 61, 37, 67, 91, 25  
G3  181, 211, 70, 199, 26, 52, 157, 38, 28, 155, 167, 82  
G4  137, 110, 124, 122, 161, 193  
G5  100, 56, 74, 227, 98, 173, 185, 229, 217, 203, 44, 88  
G6  218, 230, 194, 188, 152, 164  
Ha  102, 60, 195, 90, 153, 165  
Hb1  30, 86, 54, 147, 135, 149  
Hb2  45, 99, 101, 89, 57, 75  
Hb3  210, 166, 154, 156, 180, 198  
Hb4  169, 108, 106, 120, 201, 225  
I1  22, 151  
I2  146, 134, 182, 158, 214, 148  
I3  121, 109, 41, 107, 73, 97  
I4  233, 104  
J1  150  
J2  105 
Appendix B Canalization and Input Inversion Properties of the Boolean Functions
The following table shows the fractions of canalizing inputs and the input inversion value for each of the 256 Boolean functions with inputs. They are grouped according to their 46 classes that occur in the evolution of finite size networks.
Class  Function  
A  1  1  1  0  
B1  0  3/4  6  
B2  0  1/2  3/4  2  
B3  0  1/2  3/4  
B4  0  1/2  3/4  
Ca1  0  1/3  2/3  8  
Ca2  0  1/3  2/3  0  
Ca3  0  1/3  2/3  
Cb1  0  1/3  2/3  8  
Cb2  0  1/3  2/3  
D1  0  1/6  7/12  10  
D2  0  1/6  7/12  6  
D3  0  1/6  7/12  2  
D4  0  1/6  7/12  
D5  0  1/6  7/12  
D6  0  1/6  7/12  
E1  0  1/6  1/2  4  
E2  0  1/6  1/2  4  
E3  0  1/6  1/2  
E4  0  1/6  1/2  
Fa1  0  0  1/2  0  
Fa2  0  0  1/2  0  
Fb1  0  0  1/2  12  
Fb2  0  0  1/2  4  
Fb3  0  0  1/2  
Fb4  0  0  1/2  
Fc1  0  0  1/2  8  
Fc2  0  0  1/2  0  
Fc3  0  0  1/2  
G1  0  0  5/12  2  
G2  0  0  5/12  2  
G3  0  0  5/12  2  
G4  0  0  5/12  
G5  0  0  5/12  
G6  0  0  5/12  
Ha  0  0  1/3  0  
Hb1  0  0  1/3  4  
Hb2  0  0  1/3  4  
Hb3  0  0  1/3  
Hb4  0  0  1/3  
I1  0  0  1/4  6  
I2  0  0  1/4  
I3  0  0  1/4  2  
I4  0  0  1/4  
J1  0  0  0  0  
J2  0  0  0  0 
References
 Kauffman (1969) S. A. Kauffman, J. Theor. Biol. 22, 437 (1969).
 Kauffman (1993) S. A. Kauffman, The origins of order (Oxford University Press, New York, 1993).
 Aldana et al. (2003) M. Aldana, S. Coppersmith, and L. P. Kadanoff, in Perspectives and Problems in Nonlinear Science, edited by E. Kaplan, J. Marsden, and K. Sreenivasan (SpringerVerleg, Berlin, 2003).
 Drossel (2008) B. Drossel, in Reviews of Nonlinear Dynamics and Complexity (to be published), edited by H. G. Schuster (WileyVCH, Berlin, 2008), eprint arXiv:0706.3351.
 Bornholdt (2005) S. Bornholdt, Science 310, 449 (2005).
 Albert and Othmer (2003) R. Albert and H. G. Othmer, J. Theor. Biol. 223, 1 (2003).
 Li et al. (2004) F. Li, T. Long, Y. Lu, Q. Quyang, and C. Tang, Proc. Natl. Acad. Sci. U.S.A. 101, 4781 (2004).
 Shmulevich et al. (2005) I. Shmulevich, S. A. Kauffman, and M. Aldana, Proc. Natl. Acad. Sci. U.S.A. 102, 13439 (2005).
 Lau et al. (2007) K. Y. Lau, S. Ganguli, and C. Tang, Phys. Rev. E 75, 051907 (2007).
 Kauffman and Smith (1986) S. A. Kauffman and R. G. Smith, Physica D 22, 68 (1986).
 Bornholdt and Sneppen (1998) S. Bornholdt and K. Sneppen, Phys. Rev. Lett. 81, 236 (1998).
 Paczuski et al. (2000) M. Paczuski, K. E. Bassler, and A. Corral, Phys. Rev. Lett. 84, 3185 (2000).
 Bornholdt and Rohlf (2000) S. Bornholdt and T. Rohlf, Phys. Rev. Lett. 84, 6114 (2000).
 Bornholdt and Sneppen (2000) S. Bornholdt and K. Sneppen, Proc. R. Soc. Lond. B 267, 2281 (2000).
 Luque et al. (2001) B. Luque, F. J. Ballesteros, and E. M. Muro, Phys. Rev. E 63, 051913 (2001).
 Bassler et al. (2004) K. E. Bassler, C. Lee, and Y. Lee, Phys. Rev. Lett. 93, 038101 (2004).
 Bassler and Liu (2005) K. E. Bassler and M. Liu, in 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).
 Liu and Bassler (2006) M. Liu and K. E. Bassler, Phys. Rev. E 74, 041910 (2006).
 Szejka and Drossel (2007) A. Szejka and B. Drossel, Eur. Phys. J. B 56, 373 (2007).
 Rohlf (2007) T. Rohlf (2007), eprint arXiv:0708.1637.
 Braunewell and Bornholdt (2007) S. Braunewell and S. Bornholdt (2007), eprint arXiv:0707.1407.
 Challet and Zhang (1997) D. Challet and Y. C. Zhang, Physica A 246, 407 (1997).
 Reichhardt and Bassler (2007) C. J. O. Reichhardt and K. E. Bassler, J. Phys. A: Math. Theor. 40, 4339 (2007).
 Waddington (1942) C. Waddington, Nature 150, 563 (1942).
 Rutherford and Lindquist (1998) R. Rutherford and S. Lindquist, Nature 396, 336 (1998).
 Quietsch et al. (2002) C. Quietsch, T. A. Sangster, and S. Lindquist, Nature 417, 618 (2002).
 Bergman and Siegal (2003) A. Bergman and M. Siegal, Nature 424, 549 (2003).
 Wagner (2005) A. Wagner, Robustness and Evolvability in living systems (Princeton University Press, New Jersey, 2005).
 Harrison (1963) M. A. Harrison, SIAM 11, 806 (1963).
 Derrida and Pomeau (1986) B. Derrida and Y. Pomeau, EuroPhys. Lett. 1, 45 (1986).
 Derrida and Weisbuch (1986) B. Derrida and G. Weisbuch, J. Phys. 47, 1297 (1986).
 Flyvbjerg (1988) H. Flyvbjerg, J. Phys. A 21, L955 (1988).
 Bastolla and Parisi (1997) U. Bastolla and G. Parisi, J. Theor. Biol 187, 117 (1997).
 Bastolla and Parisi (1998) U. Bastolla and G. Parisi, Physica D 115, 203 (1998).
 Fox and Hill (2001) J. J. Fox and C. C. Hill, Chaos 11, 809 (2001).
 Aldana (2003) M. Aldana, Physica D 185, 45 (2003).
 Aldana and Cluzel (2003) M. Aldana and P. Cluzel, Proc. Natl. Acad. Sci. U.S.A 100, 8710 (2003).
 Socolar and Kauffman (2003) J. E. S. Socolar and S. A. Kauffman, Phys. Rev. Lett. 90, 068702 (2003).
 Samuelsson and Troein (2003) B. Samuelsson and C. Troein, Phys. Rev. Lett. 90, 098701 (2003).
 Kaufman et al. (2005) V. Kaufman, T. Mihaljev, and B. Drossel, Phys. Rev. E 72, 046124 (2005).
 Greil and Drossel (2005) F. Greil and B. Drossel, Phys. Rev. Lett. 95, 048701 (2005).
 Drossel (2005) B. Drossel, Phys. Rev. E 72, 016110 (2005).
 Drossel et al. (2005) B. Drossel, T. Mihaljev, and F. Greil, Phys. Rev. Lett. 94, 088701 (2005).
 Moreira and Amaral (2005) A. A. Moreira and L. A. N. Amaral, Phys. Rev. Lett. 94, 218702 (2005).
 Arthur (1994) W. B. Arthur, Am. Econ. Rev. 84, 406 (1994).
 Bak and Sneppen (1993) P. Bak and K. Sneppen, Phys. Rev. Lett. 71, 4083 (1993).
 Pólya (1937) G. Pólya, Acta. Math. 38, 145 (1937).
 Lo et al. (2000) T. S. Lo, P. M. Hui, and N. F. Johnson, Phys. Rev. E 62, 4393 (2000).
 Hart et al. (2000) M. Hart, P. Jefferies, P. M. Hui, and N. F. Johnson, Eur. Phys. J. B 20, 547 (2000).
 Anghel et al. (2004) M. Anghel, Z. Toroczkai, K. E. Bassler, and G. Korniss, Phys. Rev. Lett. 92, 058701 (2004).
 Milo et al. (2002) R. Milo, S. S. Orr, S. Itzkovitz, N. Kashtan, D. Chklovskii, and U. Alon, Science 298, 824 (2002).
 ShenOrr et al. (2002) S. ShenOrr, R. Milo, S. Mangan, and U. Alon, Nature Genetics 31, 64 (2002).