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

Abstract

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.

1. 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. An example of this is Griffiths phases [20], which arise when a certain degree of spatial disorder is applied. Under this condition, the critical temperature, $T_c$, at which finite-size scaling is observed, becomes a wide range of T [21]. Griffiths phases and other rare region effects have been studied in more general models, such as complex networks [22,23,24]. Quenched disorder can be replaced by topological heterogeneity in networks with finite topological dimension [25]. These results have a wide range of implications for propagation phenomena and other dynamical processes in networks. Much of network science [26] studies the effect of non-trivial (i.e., heterogeneous) topologies in different properties of complex systems.

Heterogeneity is pervasive in many phenomena, from public transportation systems [27,28] to social and biological networks [29,30]. Twenty years ago, scientists started exploring alternative topologies in RBNs [31]. In 2003, Aldana showed that RBNs with a scale-free topology expand the advantages of the critical regime into the ordered phase [32], i.e., a scale-free topology expands the range of the critical regime. Later, Gershenson and colleagues studied the role of redundancy in the robustness of RBNs. They concluded that the redundancy of nodes (a particular type of heterogeneity) prevents mutations from propagating in RBNs [33]. Heterogeneity is common, but how relevant is it for the dynamics of complex systems? It seems the answer is more than what was previously thought [34], also highlighting the importance of nonlinearity. This principle (heterogeneity extends criticality) seems to apply to a broad range of situations.

One example of structural heterogeneity is modularity [35,36], as it provides a level of organization that promotes at the same time robustness and evolvability. Since changes do not propagate as easily between modules in the case of RBNs, modules broaden the range of the critical regime towards the chaotic phase [37]. Another example of the importance of structural heterogeneity is found in the evolution of biological networks. By modifying the physical interactions in them and the way in which constraints are imposed on the states of their nodes, interesting dynamic changes are observed [38]. It is also possible to take advantage of heterogeneity in human-built structures. For instance, power grid networks can produce heavy-tailed distributions without a self-organizing mechanism [39]. This demeanor is due to the fact that power grids are weighted, hierarchical modular networks and particular heterogeneous topologies can be exploited to avoid global failures.

Spatial disorder generates Griffiths phases, which, in analogy to critical points, are characterized by a slow relaxation of the order parameter and divergences of quantities such as the susceptibility. However, for dimensions greater than or equal to 2, it is possible to induce Griffiths phases through time heterogeneity using absorbing states [40]. Temporal Griffiths phases are characterized by generic power-law scaling of some magnitudes and generic divergences of the susceptibility. Temporal heterogeneity is also prevalent in rank dynamics [41], where a few more relevant elements change slower than most less relevant elements in a variety of systems, possibly facilitating systems at the same time with robustness and adaptability.

The above is a non-exhaustive list of relevant examples of heterogeneity. Some years ago, Wang and Zhou found that neural networks with a modular structure can extend the parameter region of coupling strength over which critical states are reached compared to non-modular networks [42], suggesting that the brain may take advantage of the hierarchical modular structure to attain criticality and enhanced function. One year after that, Moretti and Muñoz showed the existence of Griffiths phases in synthetic hierarchical networks and also in empirical brain networks such as the human connectome [43]. According to them, stretched critical regions produce higher functionality in a generic way, facilitating the task of self-organizing, adaptive, and evolutionary mechanisms that select criticality. Recently, Ratnayake and colleagues pointed out that complex networks with heterogeneous nodes have a greater robustness as compared to networks with homogeneous nodes [44]. Furthermore, Sormunen, et al. showed adaptive networks with critical manifolds that can be navigated as changes are made to the parameters [45], i.e., criticality can be associated with a manifold in a multidimensional system, instead of being restricted to a single value.

2. Methods

2.1. 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,\dots$ 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\ge 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 ${\left(f_a\right)}_{a\in 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.

A Boolean network is called random if the assignment $a\mapsto 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 $\varphi :V⟶\{0,1\}$, $a\mapsto {\varphi }_a$, is called a state of the random Boolean network on D. The value ${\varphi }_a$ is called the state of a. The updating function $F\left(\varphi \right)$ of a state $\varphi$ is the function $F:V⟶\{0,1\}$, $a\mapsto {\varphi }_a^{\prime }$, defined as

$${\varphi }_a^{\prime }=f_a\left({\varphi }_{a_1},\dots ,{\varphi }_{a_k}\right)$$

(1)

Given the complexity exhibited by real genetic networks, Kauffman assumed the following four conditions in his model:

• Each gene is regulated by exactly K other genes.
• The K genes that regulate each node are chosen randomly using a uniform probability.
• All genes are updated synchronously, i.e., at the same timescale [47].
• 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\to \infty$. In fact, for each value of p, there exists a critical connectivity $K_c={\left[2p(p-1)\right]}^{-1}$ such that all perturbations in the initial state vanish if $K<{\left[2p(p-1)\right]}^{-1}$ and a small perturbation in the initial state propagates throughout the system for $K>{\left[2p(p-1)\right]}^{-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.

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.

A generalized Boolean network on a directed graph D consists of a family ${\left(f_a\right)}_{a\in V}$ of functions $f_a:{\{0,1\}}^{k_a^{-}}⟶\{0,1\}$ with $k_a^{-}\ge 1$ the in-degree of a. 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 $\varphi$ of a random heterogeneous Boolean network on D is the function $\tilde{F}\left(\varphi \right):V\times \mathbb{N}⟶\{0,1\}$, defined by

$$\left(a,t\right)\mapsto \left\{\begin{array}{cc}{\varphi }_a^{\prime }\hfill & \mathrm{if}t\mathrm{is}\mathrm{a}\mathrm{multiple}\mathrm{of}{K}_a^{+}\hfill \\ {\varphi }_a\hfill & \mathrm{otherwise}\hfill \end{array}\right.$$

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.

2.2. 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:

$$I=-K\sum _{i=1}^b{p}_{i}log{p}_i.$$

(2)

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:

$$C=4I(1-I),$$

(3)

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.

2.3. 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” $\mathsf{\Delta }\mathsf{ס}$ and network perturbation $\mathsf{\Delta }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 ($tmodO=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:

$$\mathsf{\Delta }x=\frac{X·\left(T/O\right)}{N·T},$$

(4)

where $0\le \mathsf{\Delta }x\le 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:

$$\mathsf{\Delta }\mathsf{ס}=C-C_0,$$

(5)

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\le \mathsf{\Delta }\mathsf{ס}$ ≤ 1. Using (4) and (5), we define the fragility of the network as

$$\mathcal{A}=-\mathsf{\Delta }\mathsf{ס}·\mathsf{\Delta }x.$$

Since $\mathsf{\Delta }\mathsf{ס}$ ∈ [−1, 1] and $\mathsf{\Delta }x\in \left[0,1\right]$, then $\mathcal{A}\in \left[-1,1\right]$. Note that, if the complexity gain decreases with perturbations, then the system is fragile, i.e., positive values of $\mathcal{A}$ indicate fragility. On the other hand, if the complexity gain remains constant with perturbations, we say that the system is robust, from which it follows that values close to zero of $\mathcal{A}$ indicate robustness. Similarly, if the complexity gain increases with perturbations, we say that the system is antifragile; therefore, negative values of $\mathcal{A}$ indicate antifragility. Since $\mathcal{A}$ depends on $\mathsf{\Delta }\mathsf{ס}$, then $\mathcal{A}$ can have different values depending on the state of the nodes. Thus, using multiple final states, we compute the average antifragility of an RBN.

3. 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.

3.1. Homogeneity vs. Heterogeneity

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).

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 Figure 5 and Figure 6 the transitions of the states of the nodes in the network for a fixed K.

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).

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.

3.2. 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 $\sigma$ (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 $\sigma =0.288$, which is much smaller than the values of $\sigma$ 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.

Figure 7b shows the complexity curves for different means in the triangular distribution. The extreme values in this case are $\mu =0$ and $\mu =1$, and both again produce sigmoidal-like behavior. It is interesting that both cases produce such similar curves, and these are no exception; for $\mu =0.25$ and $\mu =0.75$, the curves also overlap, so we can say that there is a symmetry with respect to $\mu =0.5$. This is because of the complementarity of 0 and 1 given by the logical negation. Thus, choosing a different value for $\mu$ only matters how far away we are from $\mu =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\approx 6$, all other curves exceed the one associated with the domain $D=[0,1]$. However, after $K\approx 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.

3.3. 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.

4. 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 $\sigma$ 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 $\mu =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.