Exploring the Dynamic Organization of Random and Evolved Boolean Networks

: The properties of most systems composed of many interacting elements are neither determined by the topology of the interaction network alone, nor by the dynamical laws in isolation. Rather, they are the outcome of the interplay between topology and dynamics. In this paper, we consider four di ﬀ erent types of systems with critical dynamic regime and with increasingly complex dynamical organization (loosely deﬁned as the emergent property of the interactions between topology and dynamics) and analyze them from a structural and dynamic point of view. A ﬁrst noteworthy result, previously hypothesized but never quantiﬁed so far, is that the topology per se induces a notable increase in dynamic organization. A second observation is that evolution does not change dramatically the size distribution of the present dynamic groups, so it seems that it keeps track of the already present organization induced by the topology. Finally, and similarly to what happens in other applications of evolutionary algorithms, the types of dynamic changes strongly depend upon the used ﬁtness function.


Introduction
Many natural, social or artificial systems can be described as networks of interacting elements, where the dynamical variables are associated to the nodes of the network, and directed links represent the influence of a variable upon another one. While networks can represent different kinds of systems and relationships, let us consider here for the sake of definiteness a dynamical system. A variable x i (t) (i = 1 . . . N), which can take either continuous or discrete values, is associated to every node i of the network at time t, and a deterministic law (e.g., a differential or finite difference equation) determines the time behavior of x i . If this equation contains at least one term which depends upon x j then there is a link from node j to node i.
The study of the properties of networks has become one of the major threads in complexity science. The discovery of the widespread presence of some types of topological relationships, associated to important generic properties (like e.g., the presence of hubs), prompted several studies on topology [1][2][3]. It soon became clear that the properties of these complex systems are neither determined by the topology of the interaction network alone, nor by the dynamical laws in isolation, but they are the outcome of the interplay between topology and dynamics [4][5][6][7][8][9]. This implies that it is difficult to tell in advance what the properties of a complex system are, unless of course it belongs to a class of already well-studied cases.
are markedly different, thus showing that the sharing of a common topology induces significant correlations among the nodes-which cannot of course be observed in purely random avalanches. To the best of our knowledge, this is the first case where this kind of (dynamical) organization in truly random BNs has actually been detected.
Given the importance of critical RBNs, we have also introduced an evolutionary algorithm which ensures that a population of critical Boolean networks remains critical during its evolution. Starting from purely random Boolean networks, the population evolves according to a Genetic Algorithm (GA) in order to maximize a given fitness function; for example, the achievement of attractors with a predefined fraction of active genes. A standard GA easily performs this task by changing the bias of the set of Boolean functions (see Section 2.2 for a definition of the bias) but, in doing so, it leads the networks in the ordered region. Thus, if criticality can provide some other advantages, these are lost in evolution. We have therefore introduced a modified GA which tends to keep the evolved networks critical. It has been presented elsewhere [34], where we have also shown that it is successful in maintaining criticality while achieving its goal, in the case where the goal is a fixed fraction of active nodes (as in the example just discussed), as well as in the case where the target is that of giving rise to avalanches with a predefined fraction of "up" nodes, i.e., nodes whose activation is higher in the perturbed network than in the wild type. The modified GA is briefly recalled in Section 2.3.
Since the RI methodology has been successful in showing some organization of fully random Boolean networks, we have decided to apply it to also study the differences between RBNs and the populations of evolved networks in the two cases above. Note that the evolved networks are no longer fully random, although they are dynamically critical. As will be shown in Section 3.2, the method proves able to identify measurable differences between the RBNs and the evolved networks, thus showing its capacity to react to this form of organization. Similarly to what happens in other applications of evolutionary algorithms to RBNs, the types of changes strongly depend upon the fitness function, as briefly discussed in the final Section 4.
Note that Section 2 summarizes the results of previous works, while Section 3 presents the new results of this paper, which are further discussed in Section 4.

The Relevance Index
The Relevance Index (RI) methodology, which draws its inspiration from methods used in neuroscience [35,36], makes use of indexes based on information theory and is suitable for exploring the organization of complex systems. In particular, we utilize the index called zI (see below), which makes it possible to identify, as components of a system, sets of variables that show a high degree of internal coordination. The RI methodology has been applied with interesting results to several systems: some of them had been artificially designed in order to test the effectiveness of the technique, while others referred to interesting physical, chemical, biological, or socio-economic systems [13][14][15][16].
Let U be the set of discrete variables {x 1 , x 2 . . . x n } describing a system whose status changes along a collections of m observations, and a subset S of U composed of k elements. Its integration I(S), semidefinite positive, is defined as: where H represents the Shannon entropy (respectively, of the involved variables and of the subset as a whole). The integration of a (sub) set of variables could be interpreted as their distance from the independence condition [35]-the greater the distance, the greater the probability that the variables of interest are exchanging information. We can note that the number of variables that can be correlated is arbitrary (obviously, there are limits to this "arbitrariness": for further details, see [35]): in particular, the use of integration allows to (a) overcome the pairwise analysis (a group of four elements is not always the simple merge of six binary interactions) and (b) to highlight some kinds of nonlinear relationships otherwise not visible [11,35]. The integration depends on the size of the group under observation and the size of the alphabet of the elements that compose it (elements that determine the number of degrees of freedom of the system) and on the number m of observations: in order to rank the examined groups, therefore, we need to normalize the integration according to these terms. A useful formula could be found in [12]: where we insert within a z-score form the information that the integration multiplied twice by the number of observations follows a Chi-square distribution with d freedom degree: the L j terms denoting the cardinality of the alphabet of the j-th variable belonging to the subset S to be analyzed. Ideally all possible groups have to be tested and sorted in descending order with respect to the zI index: the independent non overlapping groups that survive a sieve action are the potentially Relevant Sets (RSs in the following) [12]. All the first RSs in the list are then selected to be collectively represented by single variables, which are introduced in place of the original ones to form a new data series: the cycle of analysis, sieving, and compacting repeats several times until the last computed zI index reaches a minimum threshold (the number of all possible subsets of a set grows indeed very rapidly: it is, however, possible to sample the subsets to be tested in various ways (including the use of optimization strategies) [37,38], while the combination of the compaction and iteration procedures helps in allowing the identification of large groups [12]). The z-score typically has the purpose of verifying whether the mean value of a distribution deviates significantly from a certain reference value: in this paper, we will assume as reference thresholds θ 3 = 3.0 and θ 5 = 5.0, which in the case of Gaussian statistics correspond, respectively, to p-values close to 1.3·10 −3 and 2.7·10 −7 .
At the end of the analysis, groups of variables of different sizes (and possibly single not aggregated variables) are obtained. This grouping is based on the observation of the behavior of the (single and collective) variables along the observations: its characteristics, therefore, depend on the system's dynamics.
A RBN, as introduced by Kauffman [23,24,39], is a dynamic system whose components (nodes in an oriented graph, each one representing a gene of a genetic regulatory network) can assume only two levels (0 and 1). If the state of a single node can be represented by x i (t), the state of the whole network can be represented by the vector X(t) = [x 1 (t), x 2 (t), . . . x N (t)]. In "classic" RBNs, each node has the same incoming connectivity, while multiple connections and auto-connections are prohibited. Each node is associated with a Boolean function, which represents the response to the signals (proteins-not represented in the model) coming from the upstream nodes. In this paper, in case of random generation of the Boolean functions, for each entry, we extract a 1 with probability p (a quantity called "bias")-and consequently the frequency of the values equal to 0 is (p − 1).
The topology and Boolean function associated to each node do not change in time, and the network dynamics is discrete and synchronous, so fixed points and cycles are the only possible asymptotic states in finite networks.
The RBN's average degree of connectivity and the type of Boolean functions influence the behavior of the system, which can show ordered or disordered dynamic regimes. In ordered regimes, there are few attractors, whose period scales as a power law with the number of nodes N of the system, while in disordered regimes, the period of the very numerous attractors scales as a more than polynomial function [23,24]. The disordered regime also shows a strong sensitivity to initial conditions, which is not present for the ordered regime.
As mentioned in Section 1, these different behaviors can be related to the value of a single index, often called the Derrida parameter λ: on average, a small perturbation tends to vanish in "ordered" networks (λ < 1), while in chaotic networks (λ > 1) it increases in time, and in critical networks (λ = 1) it tends to maintain its size. The Derrida parameter can be determined by studying the average behavior in time of the distance between two close initial states of a network and taking its limiting value when the initial distance approaches its minimum value (dynamical sensitivity) [20,47,48]. However, direct measurement of the dynamical Derrida parameter in physical biological systems may be prohibitive, due to their transient nature, but it is sometimes possible, under plausible assumptions, to circumvent this problem and to infer the value of λ from other data. This is the case of the effects of gene perturbations, induced by permanently inhibiting the expression of a single gene (knock-out) and by observing the changes of the expression levels of the other genes.
In synthetic networks (by means of simulations), it is possible to evaluate these effects by comparing the behavior of a non-perturbed system ("wild type", briefly WT) with the behavior of a perturbed system, which differs from the previous one in that one of the genes, which are active in the WT system, is silenced (its Boolean function is fixed at False-"knock-out" event). It is possible to observe the two systems over time, and to mark nodes having different values in the different systems as "perturbed" nodes: the number of such nodes (different from the silenced node) provides the size of the "avalanche of changes" (or more briefly "avalanche") to which these nodes belong. In order to avoid the idiosyncrasies of single networks, each perturbation is carried out on a different RBN belonging to the same statistical "ensemble" (same number of nodes N and same bias p). Avalanches are important in this study because they constitute the experimental observations from which some aspects of the dynamic organization of systems subject to perturbation can be deduced [40].

Genetic Algorithm
In this paper, we are interested in highlighting the presence of significant dynamic structures in dynamic systems, and in observing the changes in their characteristics passing from random systems, to systems with organization, up to evolved systems-in our case, evolved by means of Genetic Algorithms (GA in the following) [49].
It should be noted that in this article, we make use of genetic algorithms in a twofold way, since we are interested in their ability (a) to solve problems, and (b) to emulate evolution.
Indeed, evolutionary and in general bio-inspired algorithms-based on the principles of the biological evolution [49,50], or on the self-organized behavior of some groups of animals such as colonies of ants, flocks of birds, and schools of fish [51][52][53][54]-play a fundamental role in complex optimization problems, allowing the development of new and robust techniques. These techniques show characteristics such as adaptation, scalability, speed, autonomy, parallelism, and fault tolerance. As consequence, several algorithms inspired by evolution [55,56] or biological behaviors of groups of animals [52,[57][58][59] have been used to find near-optimal solution to complex optimization problems in several areas. Alongside these by now "classic" themes, in recent years, there has been an interesting increase in the variety of bio-inspired algorithms, inspired by the success of plant lifestyles [60,61], Algorithms 2020, 13, 272 6 of 15 by some flight behaviours [62] or hunting behaviours in mammals [63][64][65][66], or by the presence of hierarchies [67] and opportunistic behaviours [68] in birds. Very often, this type of algorithms turns out to be the best option when an exact optimizer would not be able to give results in admissible times [55]. Bio-inspired meta-heuristic optimization algorithms are turning out to be increasingly used in different applications since they are: (i) based on simple ideas, easily implemented; (ii) typically not strongly dependent on the structure of the optimization problem; (iii) able to find a near-optimal solution; and (iv) easily applicable in different areas [52]. In this context, the GA framework plays a significant role [50].
On the other hand, for our aims, the GA's ability to mimic the natural evolutive processes [49,69] is also very interesting.
Because of the evolution, the final populations are no longer random as in the beginning-although they still preserve some randomness-so in the following, we will indicate them with the term "BN", where the character "R" for "random" has been removed.
The initial populations on which the GA acts are composed by RBNs sharing the same topology, the genotype of each individual being composed of the sequence of the truth tables of the Boolean functions that guide the dynamic response of each individual's node. Due to the arbitrariness of the numbering of nodes, the fixed topology does not limit the set of networks that can be generated. We run a series of experiments with networks whose nodes have two inputs each, the first generation being composed of random systems having the Boolean function bias equal to 0.5-and therefore typically showing a critical dynamic regime (see the previous section).
Given that criticality is supposed to be an important property, it is particularly interesting to consider a modified version of the GA, presented in [34], which preserves the bias of the Boolean functions during inheritance from parents to children. This GA maintains the usual selection procedure (a roulette wheel) and has a slightly modified mutation and crossover. In particular, we use the single-point crossover (where the cut-off point is randomly chosen among those that keep the parental bias as much as possible), while the mutation hits randomly, but in a way to again get closer to the bias of the parents. It has been verified that this modified GA preserves the bias of the initial individuals: additionally, in case of populations composed of dynamically critical individuals, the property of criticality is maintained over the generations [34]. (In this article, we follow the classic criticality measure for RBNs, represented by the Derrida procedure carried out on a high number of random initial conditions [20,47].) By following this approach presented in [12], we aim to obtain two populations of BNs, evolved with two different purposes: in one case, we search for the highest possible number of active nodes on attractors (a non-trivial task, wanting to keep the BNs in critical regime-Fit1 in the following), while in the other, we want avalanches following knock-out events to be composed of nodes whose response to stress is more of an increase than a decrease in activity (Fit2 in the following), a feature present in knock-out events in living systems [70]. (In particular: Fit1 identifies the system attractors, computes the average value of each node among each attractor, and computes the fitness function as the fraction of average values greater than 0.5; Fit2 reaches a random chosen attractor, computes the average value of each node among the attractor, and performs knock-out events on nodes having positive averages, recording the fraction of nodes (downstream from the silenced ones) that present an increase in activity with respect their state in the unperturbed situation: the fraction is the fitness value. See [12] for more details.) As already commented, in each GA run, the population of BN is made up of networks sharing the same topology, the genotype being composed by the sequence of Boolean functions that determine the dynamic response of each node. (Due to the arbitrariness of the position of nodes, the fixed topology does not limit the set of networks that can be generated.) For each fitness, we made 50 GA runs; Table 1 presents the parameter values used during simulations. (In the cases we analyzed, 100 generations were enough to reach a maximum fitness that could not be further improved, while the Algorithms 2020, 13, 272 7 of 15 increase in average fitness also stopped. We made some sample of runs even considerably longer, without observing substantial changes.) Table 1. The parameter values used during simulations. In order to characterize the systems under analysis, when necessary (searching for attractors, determination of the Derrida parameter), we used 10,000 initial conditions. The used BNs are composed by nodes having two inputs each: RBNs with this topology have a critical regime if their bias is equal to 0.5.

System
Parameter Value

Random Avalanches and Avalanches in Random Boolean Networks
A well-known technique in the study of gene regulation is that of inducing a permanent clamping of the activation value of a gene to a fixed value, e.g., by silencing it (knock-out events) [70]. The genes involved in each avalanche can react to the perturbation by increasing or decreasing their activity, situations identified later in the article by the terms UP and DOWN.
We then performed knock-out experiments in 50 RBNs each composed by 50 nodes-following an approach already developed in [12]-blocking to zero, one at a time, each node whose average activity on the attractor was greater than 0.5. In each of the analyzed systems, we obtained from 18 to a maximum of 45 avalanches, with an average on the total of 50 RBNs very close to 30 avalanches; in doing this in each system on average, about 40 nodes were affected by at least one avalanche ( Figure 1, "RBN" system class). (For presentation convenience, in the following, we will use the term "RBN" to indicate in general the framework of the Random Boolean Networks, to indicate a system "class" (the particular collection of RBNs we use in this paper), and to indicate a single exponent belonging to the aforementioned class. In the rare case of ambiguity, we will explicitly manage the different situations.) We then created 50 groups of avalanches, each group being composed by 30 avalanches with a size distribution like the one found in the 50 RBNs; each node belonging to an avalanche can have UP or DOWN value with equal probability (the typical situation for RBNs). In these random avalanche groups, typically all nodes were affected by at least one avalanche (Figure 1, "RND" system class). (Let us discuss first the results regarding these first classes of systems, by leaving to next sessions the introduction and discussion of the two classes of evolved systems (indicated as Fit1 and Fit2 in Figure 1.)) In each RBN (and in each group of the RND class) the nodes constitute the variables to be grouped, depending on their behavior in the m observations, each observation being a complete measure of response to the initial perturbation, that is, an avalanche. The final distribution of the groups identified by the RI analysis based on the index zI (i.e., the final RSs) provides us with clues about the dynamic organization of each system-that is, how the nodes form dynamically integrated groups, influencing in such a way the progress of the perturbations. in doing this in each system on average, about 40 nodes were affected by at least one avalanche (Figure 1, "RBN" system class). (For presentation convenience, in the following, we will use the term "RBN" to indicate in general the framework of the Random Boolean Networks, to indicate a system "class" (the particular collection of RBNs we use in this paper), and to indicate a single exponent belonging to the aforementioned class. In the rare case of ambiguity, we will explicitly manage the different situations.) We then created 50 groups of avalanches, each group being composed by 30 avalanches with a size distribution like the one found in the 50 RBNs; each node belonging to an avalanche can have We can perform the RI analysis by using in succession the two previously introduced threshold levels: in the case of the θ 5 threshold, we are focusing our attention on the most dynamically active elements, while the θ 3 threshold has so far been used as a level below which it is no longer correct to carry out mergers (create larger collective variables) [12]. Figure 2a shows some of the main average characteristics of the obtained RSs' distributions by using the θ 5 threshold. We can note that, despite the fact that the avalanches in RND systems involve a very high number of nodes, the nodes present within a RS are significantly less than the nodes present within a RS in RBN systems. UP or DOWN value with equal probability (the typical situation for RBNs). In these random avalanche groups, typically all nodes were affected by at least one avalanche (Figure 1, "RND" system class). (Let us discuss first the results regarding these first classes of systems, by leaving to next sessions the introduction and discussion of the two classes of evolved systems (indicated as Fit1 and Fit2 in Figure 1.)) In each RBN (and in each group of the RND class) the nodes constitute the variables to be grouped, depending on their behavior in the m observations, each observation being a complete measure of response to the initial perturbation, that is, an avalanche. The final distribution of the groups identified by the RI analysis based on the index zI (i.e., the final RSs) provides us with clues about the dynamic organization of each system-that is, how the nodes form dynamically integrated groups, influencing in such a way the progress of the perturbations.
We can perform the RI analysis by using in succession the two previously introduced threshold levels: in the case of the θ5 threshold, we are focusing our attention on the most dynamically active elements, while the θ3 threshold has so far been used as a level below which it is no longer correct to carry out mergers (create larger collective variables) [12]. Figure 2a shows some of the main average characteristics of the obtained RSs' distributions by using the θ5 threshold. We can note that, despite the fact that the avalanches in RND systems involve a very high number of nodes, the nodes present within a RS are significantly less than the nodes present within a RS in RBN systems. In other words, in RND systems, most of the nodes show absence of correlation with other nodes (a clue indicating absence of dynamic organization). Furthermore, it is possible that some group of nodes is present by chance, but their sizes compared to those of the groups in RBN systems are extremely limited, as evidenced by the comparison between the average size of the maximum and median RS of the two system classes.
The same observations can be made by analyzing the systems with the θ3 threshold (Figure 2b), with the obvious remark that the total number of nodes belonging to at least one RS increases; the average and median RSs' sizes also increase, but the ordering between the classes remains unchanged. The same holds in all our experiments: therefore, for the sake of simplicity, in the following, we will always use the θ5 threshold.
As a final comment of this section, we can note that the propagation of perturbations in random systems is considerably different from having perturbations themselves random: the introduction of a topology-even if it is itself random-introduces an strong element of dynamic organization. In other words, in RND systems, most of the nodes show absence of correlation with other nodes (a clue indicating absence of dynamic organization). Furthermore, it is possible that some group of nodes is present by chance, but their sizes compared to those of the groups in RBN systems are extremely limited, as evidenced by the comparison between the average size of the maximum and median RS of the two system classes.
The same observations can be made by analyzing the systems with the θ 3 threshold (Figure 2b), with the obvious remark that the total number of nodes belonging to at least one RS increases; the average and median RSs' sizes also increase, but the ordering between the classes remains unchanged. The same holds in all our experiments: therefore, for the sake of simplicity, in the following, we will always use the θ 5 threshold.
As a final comment of this section, we can note that the propagation of perturbations in random systems is considerably different from having perturbations themselves random: the introduction of a topology-even if it is itself random-introduces an strong element of dynamic organization.

Static and Dynamic Characteristics in Evolved Systems
The evolved systems present a dynamic organization that is certainly different from that present in random systems. Interestingly, this difference is not always evident using structural analyses alone.
A first difference with respect to random system with structure (the RBN class) can be seen in Figure 1 (Section 3): the increased number of genes active on attractors in systems evolved with Fit1 fitness allows a higher number of knock-out events; nevertheless, the increase in the number of nodes involved in avalanches does not increase in proportion. On the contrary, systems evolved with Fit2 fitness present very few nodes susceptible of knock-out events. These differences in external behavior, however, do not provide clues to identify the internal dynamic organization of these systems.
Remember that all individuals of each GA run share the same initial random topology: thus, the action of the GA involves frequency distribution of the Boolean functions, and/or the adjacency map of all the possible pairs of Boolean functions.
Each of the 16 Boolean functions has 50 chances of appearing in an RBN, and in each group, there are 50 RBNs: indeed, out of the 2500 potentially possible presences, each Boolean function appears less than 250 times. To verify the distance from a uniform distribution situation (156 presences for each Boolean function), it is therefore possible to apply the Chi-square test [71].
The Chi-square test indicates that the distribution of the Boolean functions in the case of Fit1 systems is significantly anomalous (a values of 286 vs. a significance threshold of 25): by analyzing the observed frequency of each Boolean function vs. the expected ones (in case of uniform random distribution), 10 out of 16 functions show a significant deviation at the level of 1% (Figure 3a). The Chi-square test conducted on the Boolean function adjacency matrix (given the observed frequencies of each Boolean function) highlights a further deviation (a value of 513 vs. a significance threshold of 291). The frequencies of the pairs of Boolean functions deviate strongly from the expected ones: that is, the dynamical organization of the system derives from an adjustment of both the frequencies of the Boolean functions and their linking. In Figure 3a, the Boolean functions of Fit1 systems are in order by final bit, so as to highlight the decrease in frequency of the functions having a "0" output in correspondence with an input "1,1", and an increase in the functions having a "1" output in the same case. This observation is compatible with the fact that the asymptotic states of these systems must have a high number of active genes. (This observation is compatible with a remark present in [34], albeit this latter remark provides a lower evidence because of the smaller number of analyzed systems.) The results are markedly different in the case of Fit2 systems. The Chi-square test indicates that the distribution of the Boolean functions is only slightly anomalous (a values of 31 vs. a significance threshold of 25): by analyzing the observed frequency of each Boolean function vs. those predicted in case of random uniform distribution, only 1 out of 16 functions show a significant deviation at the level of 1% (Figure 3b). The Chi-square test conducted on the Boolean function adjacency matrix (given the observed frequencies of each Boolean function) does not evidence deviations (a value of 269 vs. a significance threshold of 291). That is, the static analysis does not evidence almost any significant modification with respect to a random organization.
Interestingly, a more direct observation of the dynamic organization highlights some characteristics of the systems under investigation. Figure 4a has a structure like that of Figure 2a: however, here we have also inserted the values of the Fit1 and Fit2 systems. We can note that the evolved systems maintain some characteristics of the RBN class: the nodes present within a RS are more than the nodes present within a RS in RND systems, as well as the mean and median size of the RSs. This situation confirms and better defines our earlier statement, that even the mere introduction of a topology-even if it is itself random-introduces a strong element of dynamic organization. 291). The frequencies of the pairs of Boolean functions deviate strongly from the expected ones: that is, the dynamical organization of the system derives from an adjustment of both the frequencies of the Boolean functions and their linking. In Figure 3a, the Boolean functions of Fit1 systems are in order by final bit, so as to highlight the decrease in frequency of the functions having a "0" output in correspondence with an input "1,1", and an increase in the functions having a "1" output in the same case. This observation is compatible with the fact that the asymptotic states of these systems must have a high number of active genes. (This observation is compatible with a remark present in [34], albeit this latter remark provides a lower evidence because of the smaller number of analyzed systems.)  The results are markedly different in the case of Fit2 systems. The Chi-square test indicates that the distribution of the Boolean functions is only slightly anomalous (a values of 31 vs. a significance threshold of 25): by analyzing the observed frequency of each Boolean function vs. those predicted in case of random uniform distribution, only 1 out of 16 functions show a significant deviation at the level of 1% (Figure 3b). The Chi-square test conducted on the Boolean function adjacency matrix (given the observed frequencies of each Boolean function) does not evidence deviations (a value of 269 vs. a significance threshold of 291). That is, the static analysis does not evidence almost any significant modification with respect to a random organization.
Interestingly, a more direct observation of the dynamic organization highlights some characteristics of the systems under investigation. Figure 4a has a structure like that of Figure 2a: however, here we have also inserted the values of the Fit1 and Fit2 systems. We can note that the evolved systems maintain some characteristics of the RBN class: the nodes present within a RS are more than the nodes present within a RS in RND systems, as well as the mean and median size of the RSs. This situation confirms and better defines our earlier statement, that even the mere introduction of a topology-even if it is itself randomintroduces a strong element of dynamic organization.
Another important remark is that the action of the GA can significantly change the static structure of the system (Fit1 case), or achieve the dynamic result without producing significant changes at the "architectural" level of the static components of the system (Fit2 case). The used indicators do not show large deviations: thus, it seems that evolution keeps track of the already present organization induced by the topology.
To observe the evolutionary changes, the indicators of Figure 4a are not sufficient, and it is necessary to examine in more detail the final distribution of the sizes of the RSs of the systems Fit1 and Fit2.  Figure 4b shows with the increase of the size of the groups a rapid disappearance of the RSs in the RND systems, and a persistence of the RSs of the other classes of systems. However, the distributions of the evolved systems seem to have different trends from those of the random systems: the size of the RSs of the Fit1 systems shows a small peak, and then decreases faster than the RND systems, while the Fit2 systems show a persistent presence of large groups.
These observations are confirmed by the data present in Figure 5, where the fraction of the groups that exceeds a certain size is shown: the distribution of the sizes of the RSs of the evolved systems deviates from the distribution of the RBN class. An interesting observation is that these deviations can have different directions: towards small groups (systems of the Fit1 class) or larger groups (systems belonging to Fit2 class). Another important remark is that the action of the GA can significantly change the static structure of the system (Fit1 case), or achieve the dynamic result without producing significant changes at the "architectural" level of the static components of the system (Fit2 case). The used indicators do not show large deviations: thus, it seems that evolution keeps track of the already present organization induced by the topology.
To observe the evolutionary changes, the indicators of Figure 4a are not sufficient, and it is necessary to examine in more detail the final distribution of the sizes of the RSs of the systems Fit1 and Fit2. Figure 4b shows with the increase of the size of the groups a rapid disappearance of the RSs in the RND systems, and a persistence of the RSs of the other classes of systems. However, the distributions of the evolved systems seem to have different trends from those of the random systems: the size of the RSs of the Fit1 systems shows a small peak, and then decreases faster than the RND systems, while the Fit2 systems show a persistent presence of large groups. These observations are confirmed by the data present in Figure 5, where the fraction of the groups that exceeds a certain size is shown: the distribution of the sizes of the RSs of the evolved systems deviates from the distribution of the RBN class. An interesting observation is that these deviations can have different directions: towards small groups (systems of the Fit1 class) or larger groups (systems belonging to Fit2 class). The dynamic organization of RBNs involves a number of large RSs that is much higher than that of large RSs randomly present in random avalanches. The action of evolution modified this dynamic organization in different directions for different fitness (amplifying or reducing the dimensions of the RSs).

Conclusions
The properties of many systems composed of many interacting elements are neither determined by the topology of the interaction network alone, nor by the dynamical laws in isolation, but they are the outcome of the interplay between topology and dynamics. (An interesting exception could be that of some systems exhibiting the so-called "swarm intelligence" [54,72], in which the topology of the interactions is variable and depends on the moves of the individual agents-in these cases, the topology can be seen as dependent on both the interaction rules and the history of the system. These systems are, however, outside the present discussion.) In this work we study this relationship by carrying out a systematic comparison on a particularly interesting and understood class of systems, a well-known model of gene regulatory dynamics.
We considered four different groups of structures with increasingly complex dynamical organization (loosely defined as the emergent property of the interactions between topology and dynamics, and identified by us through the distribution of the dynamic groups acting in the system), starting from a completely random situation to a situation with structure, up to organizations with structure and subject to evolution. All the involved systems are in (or very close to) critical dynamic regimes. Whenever possible, we have used analytic methods based on the static structure of the systems under investigation, while it has always been possible to use the dynamic analytic method.
A first significant observation is that the topology per se induces a notable increase in dynamic organization: to the best of our knowledge, this is the first case where this kind of dynamical organization in truly random BNs has actually been detected. It is also interesting that this organization cannot be deduced from observations of static nature, since the RND system has no static structure.
A second noteworthy observation is the fact that the dynamic organization of the classes of systems subjected to evolution does not consist of a large change in the distribution of the RSs: the evolutionary action supplied by the GA keeps track of the already present organization induced by the topology. A suggestive hypothesis on which it is worth working is that this situation may also be valid in other systems, or more in general for the same natural evolution.
Finally, and similarly to what happens in other applications of evolutionary algorithms to RBNs, the types of dynamic changes strongly depend upon the fitness function. Fit1 systems have fewer big RSs than random systems (the RBNs), while Fit2 systems, although involving fewer nodes in the RSs with respect to the other systems, have a greater number of big RSs. It should be noted that the static Figure 5. The fraction of RSs, on the total of all RSs identified in the 50 systems (for each type) under examination, with size equal to at least the value on the X axis, linear (a) and logarithmic (b) scale. The dynamic organization of RBNs involves a number of large RSs that is much higher than that of large RSs randomly present in random avalanches. The action of evolution modified this dynamic organization in different directions for different fitness (amplifying or reducing the dimensions of the RSs).

Conclusions
The properties of many systems composed of many interacting elements are neither determined by the topology of the interaction network alone, nor by the dynamical laws in isolation, but they are the outcome of the interplay between topology and dynamics. (An interesting exception could be that of some systems exhibiting the so-called "swarm intelligence" [54,72], in which the topology of the interactions is variable and depends on the moves of the individual agents-in these cases, the topology can be seen as dependent on both the interaction rules and the history of the system. These systems are, however, outside the present discussion.) In this work we study this relationship by carrying out a systematic comparison on a particularly interesting and understood class of systems, a well-known model of gene regulatory dynamics.
We considered four different groups of structures with increasingly complex dynamical organization (loosely defined as the emergent property of the interactions between topology and dynamics, and identified by us through the distribution of the dynamic groups acting in the system), starting from a completely random situation to a situation with structure, up to organizations with structure and subject to evolution. All the involved systems are in (or very close to) critical dynamic regimes. Whenever possible, we have used analytic methods based on the static structure of the systems under investigation, while it has always been possible to use the dynamic analytic method.
A first significant observation is that the topology per se induces a notable increase in dynamic organization: to the best of our knowledge, this is the first case where this kind of dynamical organization in truly random BNs has actually been detected. It is also interesting that this organization cannot be deduced from observations of static nature, since the RND system has no static structure.
A second noteworthy observation is the fact that the dynamic organization of the classes of systems subjected to evolution does not consist of a large change in the distribution of the RSs: the evolutionary action supplied by the GA keeps track of the already present organization induced by the topology. A suggestive hypothesis on which it is worth working is that this situation may also be valid in other systems, or more in general for the same natural evolution.
Finally, and similarly to what happens in other applications of evolutionary algorithms to RBNs, the types of dynamic changes strongly depend upon the fitness function. Fit1 systems have fewer big RSs than random systems (the RBNs), while Fit2 systems, although involving fewer nodes in the RSs with respect to the other systems, have a greater number of big RSs. It should be noted that the static measurements are able to identify structural changes with respect to random systems (RBN) in Fit1 systems, while they are unable to detect anything significantly different from random systems in the case of Fit2 systems.
Alongside these three noteworthy observations, we could actually note that the fundamental contribution of this work is that of addressing the problem of quantitative identification of the dynamic organization of a system, and that of suggesting for this aim a method (the RI methodology) capable of identifying some interesting observables (the Relevant Sets). The identification of these structures allows a better observation and understanding of the phenomena under investigation, and often as consequence allows asking interesting new questions.
There are many interesting directions for future works: the two most stimulating seem to be the deepening of the possible mechanism of the re-use by evolution of the already existing dynamic organization, and the understanding of the dynamic interaction between the various RSs. Both directions require individual-level analysis.
In the first case, it is possible to keep track of the composition of the RSs in the parent-child sequence that leads from the initial random individual to the final one. In this case, it could be useful to use a learning algorithm without the horizontal exchange of information among individuals belonging to the same generation, as for example the stochastic descent, a strategy already used in RBN case [73], and verifying the distance distribution between these compositions: in the case of low distances, it can be concluded that evolution works based on the already existing dynamic organization, induced by the topology.
The second research line includes the characterization of the information exchange between the different RSs of the evolved individuals-for example by means of measures of information theory-with the consequent possibility of identifying their functional roles.