Temporal, Structural, and Functional Heterogeneities Extend Criticality and Antifragility in Random Boolean Networks

Most models of complex systems have been homogeneous, i.e., all elements have the same properties (spatial, temporal, structural, functional). However, most natural systems are heterogeneous: few elements are more relevant, larger, stronger, or faster than others. In homogeneous systems, criticality—a balance between change and stability, order and chaos—is usually found for a very narrow region in the parameter space, close to a phase transition. Using random Boolean networks—a general model of discrete dynamical systems—we show that heterogeneity—in time, structure, and function—can broaden additively the parameter region where criticality is found. Moreover, parameter regions where antifragility is found are also increased with heterogeneity. However, maximum antifragility is found for particular parameters in homogeneous networks. Our work suggests that the “optimal” balance between homogeneity and heterogeneity is non-trivial, context-dependent, and in some cases, dynamic.


Introduction
Phase transitions are found in a wide range of phenomena [1], from brain dynamics [2] to jazz ensembles [3]. Usually, critical dynamics are found near phase transitions [4][5][6]. Criticality exhibits an equilibrium between order and chaos, a balance connecting robustness, information storage, variability, computation, evolvability, and adaptability. Criticality is found in a variety of complex systems and models. Antifragility is a more recent concept [7] and can be described as a property of systems that benefit from noise, perturbations, or disorder. Antifragility has been applied in risk analysis [8], molecular biology [9], and computer science [10], among others [11].
Random Boolean networks (RBNs) were originally proposed as models of genetic regulatory networks [12][13][14]. They are useful when specific topologies are unknown. Moreover, since RBNs are general models, these archetypes offer us the possibility to explore the space of feasible living and computational systems. In classical RBNs (CRBNs), connections and functions are chosen randomly, but all nodes usually have the same number of connections and the same bias in their functions. CRBNs can exhibit three dynamical regimes: ordered, critical, and chaotic. The criticality and antifragility of RBNs depend on many different factors [15]. These can be exploited to guide the self-organization and evolvability of RBNs towards the critical regime [16].
Cellular automata are particular cases of Boolean networks where the state of a variable is determined by its spatial neighbors [17][18][19]. When modeling ferromagnets with this kind of model, it is possible to find an indication that heterogeneity modifies the dynamics.

Random Boolean Networks
Although Boolean networks were first proposed in 1969 by Stuart A. Kauffman [12], they became popular models only decades later [46]. Boolean networks are a generalization of cellular automata. Although one might think that they are a crude simplification of genetic regulatory networks, it turns out that there are several cases where this general model does, in fact, correctly describe expressed and suppressed genetic patterns.
A directed graph D consists of a set V of elements a, b, c, . . . called the nodes of D and a set A of ordered pairs of nodes (a, b), (b, c), . . . called the arcs of D. We use the symbol ab to represent the arc (a, b). If ab is in the arc set A of D, then we say that a is an in-neighbor of b, and also that b is an out-neighbor of a. We say that D is k-in regular (k ≥ 1) if every node has exactly k in-neighbors. The out-degree of a node a is the number of out-neighbors of a. Let D = (V, A) be a k-in directed graph; a family ( f a ) a∈V of functions f a : {0, 1} k −→ {0, 1} is called a Boolean network on D. Figure 1 is an example of a Boolean network on a 2 − in directed graph with 9 nodes. Illustration of a random Boolean network with N = 9 nodes and K = 2 inputs per node (self-connections are allowed). The node rules are commonly represented by lookup tables, which associate a 1-bit output (the node's future state) to each 2 K possible K-bit input configuration. The out-column is commonly called the "rule" of the node.
A Boolean network is called random if the assignment a → f a is made at random by sampling independently and uniformly from the set of all the 2 2 K Boolean functions with K inputs. A function φ : V −→ {0, 1}, a → φ a , is called a state of the random Boolean network on D. The value φ a is called the state of a. The updating function F(φ) of a state φ is the function F : V −→ {0, 1}, a → φ a , defined as Given the complexity exhibited by real genetic networks, Kauffman assumed the following four conditions in his model:

1.
Each gene is regulated by exactly K other genes.

2.
The K genes that regulate each node are chosen randomly using a uniform probability.

4.
Each gene is expressed (i.e., its Boolean value is 1) with probability p and is unexpressed with probability 1 − p.
The first two imply structural homogeneity, the third temporal homogeneity, and the fourth functional homogeneity. Although these conditions might seem to limit the model, a number of interesting dynamic behaviors have been found. In particular, Derrida and Pomeau showed analytically that there exists a dynamical phase transition controlled by the parameters K and p [48], assuming the thermodynamic limit, i.e., N → ∞. In fact, for each value of p, there exists a critical connectivity K c = [2p(p − 1)] −1 such that all perturbations in the initial state vanish if K < [2p(p − 1)] −1 and a small perturbation in the initial state propagates throughout the system for K > [2p(p − 1)] −1 . The first case is referred to as the ordered phase and the second case as the chaotic phase. Figure 2 shows such a phase transition for the classical version of random Boolean networks.  [48]. The solid curve represents the critical connectivity K c = [2p(1 − p)] −1 . When K < K c , the network exhibits ordered dynamics (small shaded area below K c ), while for K > K c , the network exhibits chaotic dynamics (large shaded area above K c ). Therefore, depending on different parameters, the dynamics of RBNs can be classified as ordered, critical (near the phase transition), and chaotic. Figure 3 shows examples of these dynamics for different K values. It is important to emphasize that the model described so far is homogeneous. As can be seen in Figure 2, if criticality is only found near phase transitions, then most of the parameter space will have "undesirable" dynamics. To extend criticality (and also antifragility) in RBNs, we will add different types of heterogeneity to the system, which are described below.
Thus, a heterogeneous random Boolean network is a generalized Boolean network whose connectivities are chosen using a particular probability distribution. In order to obtain an RBN with a heterogeneous structure, the methodology proposed in [34] will be followed and the Poisson and exponential distributions will be used. Strictly speaking, both distributions are heterogeneous. However, the exponential is "more heterogeneous" than the Poisson. Therefore, we consider the former as the heterogeneous structural case and the latter as the homogeneous structural case. In all simulations using structural heterogeneity, the connectivity value (K) was used as the parameter for both the exponential distribution and the Poisson distribution. These two distributions were used to obtain the out-degree for each node.
According to the classification proposed by Gershenson [49], classical RBNs have a synchronous (homogeneous) temporality, while deterministic generalized asynchronous RBNs have an asynchronous (and therefore, heterogeneous) nature. The heterogeneous updating function of a state φ of a random heterogeneous Boolean network on D is the functionF(φ) : where t is called the discrete time step and K + a is the out-degree of a. Thus, each node is updated every number of time steps equal to its out-degree, so the more nodes one node affects, the slower it will be updated. Until now, structural and temporal heterogeneity have been introduced; the impact of these two on the dynamics of RBNs has already been studied in other works. To conclude this subsection, we present a new type of heterogeneity, which we call functional heterogeneity.
As mentioned above, one of the conditions that Kauffman considered when setting up his model was that each gene is expressed with probability p (and therefore, is unexpressed with probability 1 − p). This value of p is the same for all nodes and is fixed in the sense that it does not change throughout the dynamics; in fact, it can be shown that, if p = 0.5, then K c = 2. It is at this value that critical dynamics are observed (near the transition between chaos and order). Functional heterogeneity consists of assigning a probability p a to each node in the network. For this, we take the values from some probability distribution whose mean is a fixed value, in this work 0.5. Thus, each gene is expressed with a different likelihood, and as we will see later, this enriches the dynamics.

Criticality
Complexity, like life or consciousness, has no fully accepted definition [50]. Probably, if one asks 100 complexity scientists for the definition of complexity, they will receive 150 different answers. There is a broad variety of notions and measures of complexity proposed under different contexts [51]. Regardless of the area of study, a partial definition of complexity is expected to be able to capture the transfer of information between components and relate it directly to other abstractions such as emergence or self-organization [52,53].
Inspired by [54], in this paper, we make use of a complexity measure based on Shannon's entropy [55]. Let b be the length of the alphabet (for all the cases considered in this paper, b = 2) and K = 1/ log 2 b a normalization factor; the Shannon entropy is given by the following expression: where p i is the probability of symbol i appearing in a string. Actually, chaotic dynamics are characterized by a high I, while ordered (static) dynamics are characterized by a low I. Following [56], since the critical dynamics is characterized by having a balance between order and chaos, the complexity will be calculated through the product between I and its complement: where the constant 4 is added to normalize the measure to [0, 1] [57]. More precisely, the complexity of each node is calculated from its temporal series, and then, the complexities of all nodes are simply averaged to calculate the complexity of the network. This measure is maximized at phase transitions, and its value decreases gradually, so it is easier to calculate compared to other complexity measures, such as Fisher's information. To qualitatively estimate the amount of criticality in a random Boolean network, we adopt the methodology presented in [34] and construct the complexity curves with respect to average connectivity K. Thus, if the maximum expected complexity is shifted, then the criticality is extended in the sense that the area under the curve is broadened. In other words, critical-like dynamics are found for a broader range of K values. As shown in [34], when structural or temporal heterogeneity are added to a network, then criticality increases.

Antifragility
According to Taleb, antifragility is the ability to improve the capacity of a system in the face of external disturbances [7]. Antifragility should not be confused with resilience or robustness, as these two allude to the system maintaining its functionality or properties after perturbations, whereas antifragility necessarily implies a benefit of the system in the face of disturbances. In Greek mythology, we find antifragility in the Lernaean Hydra, whose virtue of regenerating two heads for every one that was amputated reflects the idea of antifragility. In the natural world, a canonical example is found in the immune system, whose ability to strengthen itself when exposed to pathogens can be seen as antifragile.
Although the concept of antifragility has been applied in many areas of knowledge, few antifragility measures have been developed to date. Following [58], we define the antifragility of an RBN as the product between network "complexity gain" ‫ס∆‬ and network perturbation ∆x. Suppose you have a random Boolean network with N nodes. We express the network perturbations due to external stressors as the change of node states in the RBN. From the N nodes, we flip the states of X randomly chosen nodes (from 1 to 0 or from 0 to 1, without distinction), where the perturbations are added with a frequency O during the simulation runtime T. In other words, perturbations are added whenever the time step t is divisible by O (t mod O = 0). To exemplify let X = 3, O = 5, and T = 100. This implies that the states of three randomly chosen nodes in each configuration are inverted every five time steps until the simulation runtime is 100. By comparing the state transitions of the original network and its perturbed network, we can observe how the perturbations propagate in time [58]. Thus, if T is the total simulation time, the degree of perturbations is defined as follows: where 0 ≤ ∆x ≤ 1.
To define the level of complexity gain in an RBN, we can interpret each node in the network as an agent. Since the higher the complexity, the better the balance between adaptability and robustness, we say that the network has gained complexity when it reaches a high level of complexity. Using the complexity measure described above (Equation (3)), the degree of complexity gain in the network is quantified as: where C 0 is the complexity of the network before adding disturbances and C is the complexity of the network after adding disturbances. Since the value of the complexity ranges between 0 and 1, it follows that −1 ≤ ‫ס∆‬ ≤ 1. Using (4) and (5), we define the fragility of the network as A = −∆‫∆·ס‬x.

Results and Discussion
To have simple notation, we use "Ho" for homogeneous, "He" for heterogeneous, "S" for structural, "T" for temporal, and "F" for functional. First, we compare the complexity curves of the eight cases that we described above, and then, we study the variations of that outcome using networks with different numbers of nodes. Then, we visualize the transitions under certain conditions for each case. We also vary some functional parameters and one temporal parameter to study which combination of them extends criticality even further. These include the standard deviation, the mean, and the domain of the functional distribution, and also, we vary the temporal update. Finally, we explore what is the effect of different heterogeneities on antifragility. Details on how the code works can be found in Appendix A of this article. Figure 4 shows the complexity vs. connectivity curves for the eight cases, which are the combinations of structural, temporal, and functional homogeneity/heterogeneity. For the temporal heterogeneity, we considered an out-degree updating. The distribution used for the functional heterogeneity is triangular with a mean of 0.5, and its domain is the whole interval [0, 1]. The curves that reach lower values of complexity represent the totally homogeneous case (HoSHoTHoF), while the curves that reach higher values of complexity are the threefold heterogeneity case (HeSHeTHeF).

Homogeneity vs. Heterogeneity
In CRBNs (3Ho), the size N of the network does not seem to affect the results that much, but this is not the case, as heterogeneity is added. The K needed to ensure that the triple heterogeneity (3He) maximizes criticality depends on the number of nodes. For N = 100 (Figure 4a), the curves look smooth, and the minimum K at which 3He is greater than all cases is approximately K = 6. In the case N = 75 (Figure 4b), the order of the curves is preserved with respect to the previous case, but the curves look flattened. Moreover, the K required for 3He to be maximal is close to K = 7. Although, for N = 50 (Figure 4c), the order of the curves is also preserved, the K required for 3He to be the upper curve is now a value close to K = 8. Finally, we observed that, at N = 25 (Figure 4d), the curves are less smooth. Furthermore, by the trend of the HeTHeSHeF (3He) curve, we can state that the hierarchy is still fulfilled, but now, a larger K would be needed to better visualize it.
If we focus on a single type of heterogeneity, structural (HoTHeSHoF) leads to the highest increase of complexity C, followed by temporal (HeTHoSHoF), while functional (HoTHoSHeF) has a lower maximum C than 3Ho, but tends to have a higher C for larger K.
As the number of nodes (N) increases, the minimum connectivity (K) for the triple heterogeneity to be beneficial decreases. Thus, there is a relationship between the number of nodes in the network and the critical connectivity K c for heterogeneous RBNs. Furthermore, heterogeneity seems to be additive (the more types of heterogeneity, the more criticality). Although, for very large K, we could claim that the 3He maximally extends criticality, we ran into different limitations. The first is algorithmic, since obtaining the complexity vs. connectivity curves for larger K is computationally expensive. This implies we are limited to lower ranges of K. Since the 3He curve rises for higher K, then the area under the C curve (a measure of criticality that could be used) for other cases is larger than for 3He. Even if we try to truncate these complexity curves, it may be that 3He is not the highest C for all values of N. Therefore, it is necessary to find another way to prove that the triple heterogeneity extends criticality the most. Moreover, the fact that we need a very large K to obtain maximum criticality given a specific configuration implies that, if a complex system is under certain conditions (e.g., suppose resources are scarce or limited), then a certain amount of homogeneity will be better in order to achieve higher criticality. In the view of the above and to give a qualitative proof with the advantages of the triple heterogeneity case, for each combination, we present in Figures 5 and 6 the transitions of the states of the nodes in the network for a fixed K.  Figure 4. Average complexity C of RBNs as the average connectivity K is increased for the eight cases. The labels describe the combinations used; He stands for heterogeneity, Ho for homogeneity, T for temporal, S for structural, and F for functional. The dissimilarities of these curves for different numbers of nodes in the network are shown. With N = 100 nodes (a), it is observed that, approximately from K = 6, the triple heterogeneity (3He) dominates the other cases. For N = 75 (b), the order of the curves is preserved, but the value of K needed to demarcate the dominance of 3He is higher, reaching a value close to K = 7. With N = 50 nodes (c), the curves look much flatter; despite the increase of noise in the curves, the order of the curves is still preserved, but now, an approximate value of K = 8 is needed to visualize that 3He outperforms the other cases. At N = 25 (d), too much noise accumulates in the curves; however, by following the trend of the 3He curve, it can be stated that the order will still be preserved, but now, a much larger K is required to display it. For all curves, a step of ∆K = 0.2 was used, giving a total of 45 K-connectivities. For each of these connectivity values, 1000 networks were generated, and their complexities were averaged. To obtain these results, 2000 time steps were used.
The transitions for the functionally homogeneous cases (N = 100, K = 10) are shown in Figure 5. In HoSHoTHoF (Figure 5a), we observe a clear similarity with white noise. In fact, the complexity for this case is lower than for the other cases. From HoSHeTHoF (Figure 5b) and HeSHoTHoF (Figure 5c), we can state that, in this scenario, structural heterogeneity (HeS) contributes a greater extent of criticality than temporal heterogeneity (HeT). This is despite the fact that HeT seems to give much more structure to the HoSHoTHoF white noise. Finally, by combining the effect of HeS and HeT, we obtain HeSHeTHoF (Figure 5d), and this results in transitions with more structure. Note how the complexity of this last case (C = 0.8417) is larger than the sum of the two previous cases (C = 0.2732 + 0.3603 = 0.6335). The above suggests that heterogeneity is only qualitatively (not numerically) additive; the greater the heterogeneity, the greater the complexity is (and hence, the more likely a higher criticality threshold). In the upper part of the subfigures, we indicate the distributions used for each case. "Poi" stands for Poisson, "Exp" for exponential, "Out" for out-degree, and "Non" for none. At the bottom of each subfigure, the average complexity (per iteration) for each case is shown. The iterations for HoSHoTHoF (a) resemble white noise; this coincides with the fact that the complexity obtained is very low. Adding only temporal heterogeneity (b) gives a little more structure, but the complexity is still very low. Adding only structural heterogeneity (c) does not seem to gain as much structure in the iterations, but still, the complexity is higher than in HoSHeTHoF. Adding both temporal and structural heterogeneity (d) gives more structure to the dynamics, and the complexity increases more. To generate these images, we considered a network with 100 nodes, 100 steps each (time flows downwards), and a connectivity of K = 10.
The dynamic transitions for the functionally heterogeneous cases are shown in Figure 6 (also N = 100, K = 10). The notation is essentially the same as in the previous case, with the difference that we now add "Tri" to refer to the use of the triangular distribution in the functional case. A higher complexity value is observed in HoSHoTHeF (Figure 6a) than in HoSHeTHoF (Figure 5b) and HeSHoTHoF ( Figure 5c); this implies that functional heterogeneity (HeF) outperforms HeS and HeT in extending criticality separately (remember that we focus only on K = 10). However, HeF alone does not outperform the sum of HeS and HeT. On the other hand, both HeSHoTHeF ( Figure 6b) and HoSHeTHeF (Figure 6c) are consistent with previous observations, in which HeS provides a greater extension of criticality with respect to HeT; to verify the above, it suffices to compare the complexities of these two cases. Finally, we observed that 3He (Figure 6d) has higher complexity with respect to all the previous cases, confirming once again that the more heterogeneities, the higher the complexity is. It is important to mention that these calculations are for a relatively large K. However, for smaller K (and different values of N), the above is not always true. In addition, we should also note that it is not true that the sum of the complexities for HeSHoTHoF (Figure 5c), HoSHeTHoF (Figure 5b), and HoSHoTHeF (Figure 6a) is equal to the complexity produced by HeSHeTHeF (Figure 6d). Again, heterogeneity is not numerically additive. (d) Figure 6. State transitions by iteration for functional heterogeneity. In the upper part of the subfigures, we indicate the distributions used for each case. "Poi" stands for Poisson, "Exp" for exponential, "Out" for out-degree, "Tri" for triangular, and "Non" for none. At the bottom of each subfigure, the average complexity (per iteration) for each case is shown. When only functional heterogeneity is added (a), greater complexity is observed with respect to HoSHeTHoF and HeSHoTHoF. Adding functional and structural heterogeneity (b) results in similar complexity to the HeSHeTHoF case, but having functional and temporal heterogeneity (c) keeps the complexity at an intermediate value. This is consistent with what we saw previously: HeS increases the complexity more than HeT. Finally, HeSHeTHeF (d) is the case that gives the major structure to the iterations. In fact, complexity is close to maximal for this case. To generate these images, we considered a network with 100 nodes, 100 steps each (time flows downwards), and a connectivity of K = 10.

Varying Functional and Temporal Parameters
Now that we know the effects of heterogeneous functionality in extending criticality, we can explore the parameter space associated with HeSHeTHeF and study which variations further extend criticality. Figure 7a shows the complexity curves for seven different functional distributions. Four of them are Gaussian, but with different standard deviations. The extreme case corresponds to a point-like distribution (PL), which is equivalent to the homogeneous functional case and has zero standard deviation. Approximately, for values greater than K = 6.5, the Gaussian and uniform distributions overcome the complexity of the extreme case. However, this is not the case for the triangular distribution, since for it, a much higher K is needed to extend the criticality. Since working with larger K values is computationally expensive, it is not possible at this time to know whether the sigmoid-like curves will continue to grow or at some point begin to decrease. It is clear that, for Gaussian distributions, there is a value of σ (close to 0.5) such that, after this number, the result no longer depends on the chosen standard deviation. This property of the standard deviations cannot be generalized, for example the uniform distribution has σ = 0.288, which is much smaller than the values of σ for the Gaussian distributions and still has a behavior very similar to them. Finally, we can state that, qualitatively, both the triangular and PL cases have a larger area under the curve than the other cases (which seems to extend criticality better); this is another indication that total heterogeneity is not always better. (d) Figure 7. Average complexity of RBNs as the average connectivity K is increased for variations in the functional parameters of the triple heterogeneity (3He) case. (a) shows the changes obtained by testing seven different functional distributions. Point-like indicates the extreme case when the standard deviation is zero; this is equivalent to the homogeneous functional case. When the distribution is Gaussian, convergence is observed after a certain value for the standard deviation. This behavior is not generalizable, because, with the uniform distribution, something similar is obtained, and its standard deviation is smaller. By varying only the mean of the functional distribution (b), we observe symmetric behavior, so that varying the mean only matters how far we are from µ = 0.5. (c,d) show what happens when the domain of the functional distribution is modified. In both cases, a similar order in the curves is observed, and also, as we make the interval shorter, the curves stick together. For all curves, a step of ∆K = 0.2 was used, and in every step, a total of 1000 iterations were averaged. Figure 7b shows the complexity curves for different means in the triangular distribution. The extreme values in this case are µ = 0 and µ = 1, and both again produce sigmoidal-like behavior. It is interesting that both cases produce such similar curves, and these are no exception; for µ = 0.25 and µ = 0.75, the curves also overlap, so we can say that there is a symmetry with respect to µ = 0.5. This is because of the complementarity of 0 and 1 given by the logical negation. Thus, choosing a different value for µ only matters how far away we are from µ = 0.5. Again, as we cannot reach higher values of K, we can only state that, for this interval, the smallest areas under the curve correspond to the extreme cases. This is due to the fact that they take longer to rise with respect to the other cases. Figure 7c shows the complexity curves taking different intervals as the domain for the triangular distribution. All intervals are symmetrical with respect to 0.5, and as the interval length is reduced, the associated standard deviation also decreases. Before K ≈ 6, all other curves exceed the one associated with the domain D = [0, 1]. However, after K ≈ 6, this curve is the one that surpasses the others. Because of this, it is difficult to know qualitatively which extends criticality more. However, it is clear that, as we make the interval smaller, the curves stick together. This effect is even more visible in Figure 7d, which is the same experiment, but using a uniform distribution. To make the analysis a little more quantitative, a simple integration algorithm was implemented to calculate the area under the curves. By doing the above, it was obtained that, for the triangular distribution, the interval that extends the criticality better is D = [0.1, 0.9], while, for the uniform distribution, the interval that extends the criticality better is D = [0.2, 0.8]. This implies that criticality depends not only on the length of the interval, but also on the distribution used.
So far, we have modified only the functional parameters. However, to further complement this work, we performed an experiment using two different temporal strategies. The first strategy is the out-degree that we have been using along the paper and consists of "the more connections a node has, the slower it is updated". The second strategy is called "Ceil", and this assigns to all the nodes the same maximum activation period. In order to have a more exhaustive trial, we decided to explore two types of functional distributions: the triangular and the uniform, both with a mean of 0.5 and domain [0, 1]. Figure 8 shows the results of the experiment described above. For both functional distributions, the Ceil strategy extends the criticality more compared to the out-degree strategy we have been using. Furthermore, it appears that, only for very large K, we could achieve an advantage by using out-degree updating. However, as previously discussed, this is not always possible. In general terms, we could argue that out-degree is a more heterogeneous temporal strategy compared to Ceil. However, it is the latter that gives us more criticality under a resource-constrained scenario, which again shows that maximum heterogeneity should not always be preferred.  . Average complexity of RBNs as the average connectivity K is increased for two temporal strategies using two functional distributions. For the triangular distribution, out-degree achieves a higher complexity for relatively large K. However, it is clear that, for most values of K, the Ceil strategy is dominant. For the uniform distribution, Ceil is always superior to out-degree. For very large values of K, a tie is observed, but out-degree never seems to completely outperform Ceil. For all curves, a step of ∆K = 0.2 was used, and in every step, a total of 1000 iterations were averaged.

Antifragility
So far, we have explored how heterogeneity extends criticality, which was partially studied in [34]. Another property that is fundamental to the evolution and adaptability of complex systems is antifragility, so it is worth studying whether heterogeneity also extends antifragility in some way. Figure 9 shows the variation of fragility with respect to X and O for five different connectivity values. For the 3Ho case, K = 1 represents ordered dynamics, K = 2 represents critical dynamics, and K = 3, 4, 5 represents chaotic dynamics. When the curves become negative, this implies antifragility. In the homogeneous case (Figure 9a,b), both the ordered network and the critical network show increased antifragility, which is absent in the chaotic cases. By adding a triple heterogeneity (Figure 9c,d), we observed an extension in the antifragility intervals for fragility vs. O and for fragility vs. X. Despite the extension of the antifragility intervals, we note that, for the homogeneous case, larger values of antifragility are reached, which again indicates a necessary balance between homogeneity and heterogeneity. On the other hand, with 3He for K > 2, other interesting things happen. First, we observe that the curves for K = 4 remain almost invariant to the change from homogeneity to heterogeneity. Moreover, there is a reversal of order in the curves for K = 3 and K = 5, since, for the homogeneous case, K = 3 is always the upper curve; however, for the heterogeneous case, K = 5 becomes predominant, and even regimes appear where K = 3 and K = 4 reach antifragility. These results are consistent with those found in [58]. Heterogeneity not only extends criticality, but also the region of ordered dynamics, where antifragility is highest. However, homogeneous networks are "more ordered" than heterogeneous ones. Thus, analogous to criticality, a balance between homogeneity and heterogeneity is also necessary to achieve "optimal" antifragility in complex systems. (d) Figure 9. Average antifragility of ordered, critical, and chaotic RBNs depending on X and O for the 3Ho and 3He cases. In the homogeneous case (a,b), only the ordered and critical networks achieve antifragility intervals, the former reaching much larger values. Adding triple heterogeneity (c,d) shows a reversal in the order of the previously chaotic networks, and in this case, two of the three curves for K > 2 reach antifragility intervals. Although, in the heterogeneous case, the number of curves attaining antifragility increases and the length of the intervals with antifragility is extended, in the homogeneous case, much larger antifragility values are obtained. To generate the fragility curves against O, a value of X = 40 was used. To generate the fragility curves against X, a value of O = 1 was used. For all curves, steps of ∆X = 1, ∆O = 1 were used, and in every step, a total of 1000 iterations were averaged. To obtain these results, 200 time steps were used.

Conclusions and Future Work
An open cosmological problem is the existence of a fine-tuned universe [59], where only certain values in the fundamental constants give rise to the complex structures that surround and sustain us. This has evoked a series of different explanations. One approach proposes that heterogeneity offers a natural way to extend criticality. As such, this perspective is consistent with the fact that diverse natural structures and social architectures are heterogeneous. Our work showed that, by adding some kind of heterogeneity to RBNs (structural, temporal, or functional), both criticality and antifragility are extended. The generality of this model suggests that the results obtained could be extended to different natural and social phenomena that can be modeled with networks. Although heterogeneity may not be the only way to increase the criticality or antifragility of a system, the results obtained show that, once heterogeneity is present, the dynamics will have certain evolutionary advantages. This implies that parameter fine-tuning is not necessary to obtain the benefits of criticality, making the evolution of complex systems faster and more reliable. It is also important to emphasize that, despite the benefits obtained, heterogeneity is not "the whole enchilada". Simulations show that a balance between homogeneity and heterogeneity is necessary for a complex system to thrive in a dynamic environment.
In our explorations, there are many parameters that could be varied. Certainly, we could perform a more thorough investigation of the parameter space, but there are also other interesting things to explore. The results point to the existence of a hierarchy in the heterogeneities. Each heterogeneity has a different weigh for extending criticality and antifragility. We do not have an explanation for this hierarchy. Why is HeF less effective at extending criticality than HeS? Why does HeS produce greater complexity than HeT, even though HeT seems to give more structure to node transitions? Implementing other types of heterogeneity (and their combinations) could provide answers to these questions. Could it be that the combinations between homogeneities and heterogeneities induce a hierarchy? What would be the rules governing this ranking [41]?
When varying the standard deviations in the Gaussian distribution, a clear tendency of the curves to converge after a certain value of σ was noticed. Although this behavior does not seem universal, this same phenomenon could be studied in other functional distributions. Likewise, when varying the mean, we found a symmetry with respect to µ = 0.5, and in fact, we could also investigate if this symmetry is present in other distributions and in non-Boolean networks [60]. If so, we could argue that there is a general symmetry, and as any continuous symmetry, this induces a conservation law (which can be highly non-trivial [61]). Then, there would be a general conserved quantity that governs the dynamics of the system and that can help us to better understand how complex systems work in general. Furthermore, the parameter that we investigated in a less furtive way was the temporal updating. Perhaps other temporal updating schemes (including some variations of the out-degree mechanism) could be devised to find out which is the best strategy in the temporal heterogeneous case. Moreover, when we study antifragility, we are focusing on a single form of disturbance, but it is possible to implement other types of disturbances and other ways of measuring the impact of such changes [62].
So far, we have shown that, although the presence of heterogeneity is important in extending criticality and antifragility, this rule is not entirely universal. Under certain conditions, decreasing heterogeneity and increasing homogeneity can also be beneficial in maximizing both criticality and antifragility. This is even more visible when adding functional heterogeneity, since, in this case, it can be argued that heterogeneity is not entirely additive, i.e., it is not true that the more heterogeneity, the higher the criticality threshold. Furthermore, we can appreciate a relationship between the number of nodes in the network (N) and the connectivity among them (K). For networks with few nodes, it takes a larger K to take advantage of this heterogeneity overlap. At the same time, if N is large enough, the connectivity from which we take advantage in the 3He case decreases. During the evolution of any complex system, there may be phases where conditions are not optimal for exploiting environmental resources, which reinforces the idea of balance described above. The functional and temporal variations reflect the high dependence between the variables chosen and the complexity curves obtained. For now, it seems difficult to state that a certain combination of parameters produces maximum criticality, and this is mostly due to the computational expense involved in working with large K values.
Despite their generality, it is not possible to model all natural and social phenomena using random Boolean networks. One of the major limitations of these models is that they only allow each node to have two states. However, it is not difficult to imagine dynamics in which the binary nature of the states is absent. Implementing networks with more than two states is a computational challenge, since, even in the Boolean case, the memory can be filled exponentially. Nevertheless, there are ecological models such as Roy et al.'s [63] that show the importance of heterogeneity for two-dimensional lattices, whose nodes can have three different states. With sufficient computational power, it will be possible to explore the effects of heterogeneity in much more general networks with non-binary states [64,65]. Another difficulty in sketching nature is the entanglement between its different scales. In general, networks only allow capturing the system at one of its scales. If one wishes to study a phenomenon at different scales (and the interactions between them) using networks, concepts such as causal emergence have been developed [66]. It is also possible to attack this problem by relating spatial structure and global densities;an example of this that also reflects the importance of heterogeneity was studied in ecology by Pascual et al. [67]. Finally, in many network models, the connectivity between nodes is fixed. This implies that elements can only interact with a fixed set of components throughout the entire dynamics. This dynamic rigidity can be abolished by using, for example, agent-based models and, more generally, with adaptive networks [68]. Although criticality and resilience have already been studied in these models [69], the effects of heterogeneity on them remain to be studied. to choose are the number of nodes in the network (N), the simulation time (T), and a probability p, which serves as the activation probability in the homogeneous functional case or as the mean of the functional distribution in the heterogeneous functional case. Depending on whether the fragility curve is being calculated against X or against O, a fixed value of O or X, respectively, is left as a parameter.
The program works as follows: First, a vector of five entries is introduced; the first entry corresponds to whether we are studying the heterogeneous temporal case or not (CRBN/DGARBN); the second is to give homogeneous structure to the network or not (in our case, all were considered "Complex" since this option allows us to study the functional heterogeneity in the network); the third is the type of temporal update we used (if we are in the homogeneous temporal case, we simply leave it empty); the fourth is to select the topology of the network (which can be determined by a Poisson or exponential distribution); the fifth input is the connectivity of the network; in our case, five values can be chosen: K = 1 is for an ordered network; K = 2 is for a critical network; K = 3, 4, 5 is for a chaotic network.
Once the type of network to be studied is known, the variable "number_of_iterations" is taken, and for each K (connectivity), an amount of "number_of_iterations" networks without perturbations are constructed; all of them are left to evolve, and at the end, the complexity of each of them is calculated. Then, the unperturbed networks are taken, and perturbations are added to them, resulting in the construction of "number_of_iterations" perturbed networks; all of them are left to evolve, and at the end, the complexity of each of them is calculated. Both the average complexity of the perturbed networks and the average complexity of the undisturbed networks are entered into the FRAGILITY function, then the fragilities for a given connectivity are calculated and averaged over the "num-ber_of_iterations" networks. At the end, for each connectivity, there will be an average fragility, and that is what we observe in Figure 9.