A Multilayer Structure Facilitates the Production of Antifragile Systems in Boolean Network Models

Antifragility is a property to not only resist stress and but also to benefit from it. Even though antifragile dynamics are found in various real-world complex systems where multiple subsystems interact with each other, the attribute has not been quantitatively explored yet in those complex systems which can be regarded as multilayer networks. Here we study how the multilayer structure affects the antifragility of the whole system. By comparing single-layer and multilayer Boolean networks based on our recently proposed antifragility measure, we found that the multilayer structure facilitated the production of antifragile systems. Our measure and findings can be utilized for many applications from understanding properties of biological systems with multilayer structures to designing more antifragile engineered systems.

We measured antifragility of BNs based on the change of complexity before and after adding perturbations, in which the BNs were all single-layer networks. However, numerous realworld complex systems are composed of multiple subsystems that interact with each other, which can be regarded as multilayer networks [17,18]. Here we aim at investigating how the multilayer structure affects the antifragility of the whole system by assessing the antifragility of single-layer and multilayer RBNs and comparing them.
The rest of this paper is organized as follows. In the section of "Measurement of Antifragility in single-layer and multilayer RBNs", we explain single-layer and multilayer network models, how to calculate their complexity, perturbations to networks, and how to assess the antifragility. In the section "Experiments", specific experimental designs are described. In the section of "Results and Discussion", the results about the antifragility of singlelayer/multilayer RBNs and a biological BN are mentioned. The last section summarizes and concludes the paper.

Single-layer and Multilayer Random Boolean Networks
RBNs were suggested as models of gene regulatory networks (GRNs) in cells that are present in all known living organisms [28][29][30]. Although RBNs are highly simplified models, they can greatly explain relevant properties of life and its possibilities. Accordingly, they have been actively used in many fields such as systems biology and artificial life [31][32][33][34][35][36]. In this study, a single-layer RBN represents a GRN at a single cell level, a multilayer RBN indicates coupled GRNs at a multicellular level.
A RBN is also called Boolean network, where is the number of nodes, and is the number of input links per node. Here self-links are allowed. In a RBN, the links are randomly arranged, and Boolean functions are randomly assigned to each node as well. Once the topology and Boolean logic rules are determined, they remain fixed. Each node represents a gene. The state of a node can have either 0 (off, inhibited) or 1 (on, activated), and it is updated by the states of input nodes and corresponding Boolean functions.
A state space of a RBN is the set of all possible configurations (2 N ) of a system including the transitions among them. In the state space, stationary configurations are attractors (point or cyclic), and the others converging into attractors are their basin of attraction. Depending on the structure of the state space, RBNs have ordered, chaotic, or critical regimes. The ordered and chaotic regimes indicate phases, while the critical regime refers to the phase transition boundary between them. The dynamics of RBNs can be systematically varied by K: ordered for = 1, critical for = 2, and chaotic for 3, on average under internal homogeneity (i.e., probability of being activated or inhibited) = 0.5 [37].  Our multilayer RBN model is composed of two layers: intercellular and intracellular [38,39].
In an intercellular layer, cells interact with each other and in an intracellular layer, genes interact with each other (Figure 1). All the cells have the same RBNs, and cellular topologies representing interactions between cells keep changing in each simulation. The assumption on such dynamical cellular topology is based on research showing that interacting cells continue to change by cell movements or cell growth [40].
In the multilayer RBN model, the states of all the nodes are simultaneously updated as a whole system. The specific update rules are as follows:  Communicating genes: In each RBN, communicating nodes are assigned for cell-cell interactions, which follows cell signaling in Flann et al.'s model [31]. The state of a communicating node is determined by communicating nodes of neighboring cells. If even one is activated among the communicating nodes of its neighbors, it is activated.
Here the neighbors mean source nodes from which the links originate.  The other genes: The states of the other nodes are determined in the same way as the update of the node states in a single-layer RBN mentioned above.

Complexity of RBNs
We measure complexity of RBNs based on our previous approach [41,42] as follows: (1) where is the "emergence" of node , is the probability that the state of the node is ( = 0, 1), (0 1) is the complexity of the network, and (0 1) is the average of the emergence values for all the nodes. ( ) is calculated by counting the number of 0s (1s) expressed in node during time steps. Specifically, ( ) is computed from simulation time +1 to 2 not from 1 to , which is to obtain ( ) in more stable state transitions (i.e., closer to attractors). Figure 2 shows an example calculating complexity of a RBN.
Emergence here is understood as novel information being produced, so it can be measured precisely with Shannon's information (1). Complexity indicates a balance between regularity ( and change ( [29] and for regular RBNs it is maximized at the phase transition, i.e. in the critical regime. From an information viewpoint, the regularity enables information to be preserved, and the change allows new information to be explored [42]. In the context of RBNs used as GRN models, keeping and changing the node states which point out genetic information can be connected with stability to maintain existing functions and adaptability to flexibly adapt to a new environment. Based on equation (2), an optimal balance between regularity and change reaches maximum when the emergence is 0.5 ( = 0.5 = 1). It is when the expression of state 0 or 1 is highly probable, that is, or 0.89. [41,43]. On the contrary, the complexity becomes minimum when the emergence is 0 or 1 ( = 0 or 1 = 0). It is when only one state is expressed ( or = 1; = 0) or the two states are evenly distributed ( = = 0.5; = 1). The coefficient 4 is added to normalize to the [0,1] interval.

Network Perturbations to RBNs
We perturb RBNs by flipping node states. For a RBN consisting of nodes, we randomly choose nodes and perturb them with frequency during 2 time steps. The perturbations are added whenever the time step is divisible by ( = 0). For example, = 3, = 5, and = 10 mean that we randomly choose three nodes of a network at each time step, and flip the node states every five time steps from the initial to 2 10 time steps. The perturbed three nodes are different every five time steps. Then, we calculate fragility based on the states transitions during ten time steps from = 11 to = 20.
To obtain normalized fragility values, we define the degree of perturbations as follows: where 0 1.

Antifragility of RBNs
where is the degree of "satisfaction" from perturbations, and is the degree of perturbations added to a system. As mentioned in Introduction, antifragile systems not only resist stress but also benefit from it. We describe the benefit as satisfaction. The satisfaction can be differently measured depending on systems. In this study, we calculate the satisfaction based on complexity. In other words, networks will have a higher satisfaction the closer they are to criticality. However, one can measure the satisfaction using other criteria such as performance and fitness.
If a system does not get the satisfaction from perturbations and rather is damaged, it means the system is fragile. If the system does not change against perturbations, then it is robust. If the system increases its satisfaction with perturbations, it is antifragile. Antifragility is simply represented as 1-. and are normalized to the interval [-1,1] and [0,1] respectively.
In RBNs, the satisfaction is computed based on the difference of complexity before and after adding perturbations. is computed by the following equation: where is complexity of a network before adding perturbations, and is complexity of the network after adding perturbations. The same initial states are used at = 0. Because the value of complexity is between 0 and 1, has values between -1 and 1 (-1 1). for RBNs was defined in the previous section.
If the of a RBN has a negative value, the RBN is considered antifragile. If is a positive value, the RBN is fragile. If is close to zero, the RBN is robust.
Based on equation (4), has negative values when is larger than , which means that the complexity of the system is improved by perturbations. On the other hand, has positive values when is greater than , which indicates that the complexity is reduced by the perturbations. becomes zero when is equal to . It represents that the complexity does not change with perturbations.

Experiments
We measured fragility in single-layer RBNs, multilayer RBNs, a single-layer biological BN, and a multilayer biological BN.
In the case of single-layer RBNs, we produced ordered ( = 1), critical ( = 2), and chaotic ( = 3, 4) RBNs. For each RBN, we randomly chose 10 different initial states and then investigated their state transitions during 2 400 = 800 time steps, respectively. Using the same initial states, we also looked into the state transitions of the perturbed RBNs. Comparing them during the last 400 time steps, we calculated mean of for 10 initial states. The values shown in the plots are averages calculated from 100 different RBNs per .
For multilayer RBNs, we generated three regimes of multilayer networks taking ordered, critical, chaotic RBNs. The topology at an intercellular layer was randomly determined based on the number of links randomly chosen between 1 and 81 (because the number of cells was set to nine, the intercellular network can have a maximum of 81 links.). For an individual RBN, the genes were set to eighteen and the communicating genes were set to six, which is based on the numbers of genes and cell signaling molecules in the real biological system we used (i.e., CD4+ T-cell). In the same manner as the fragility of single-layer RBNs, of multilayer RBNs was calculated.  For a biological BN, we made use of a CD4+ T cell network consisting of 18 nodes. The network for CD4+ T cell differentiation and plasticity is modeled to study how CD4+ T cells control immune responses depending on environmental signals and immunological challenges [44]. In the network, there are six cytokines which are cell signaling molecules related to cell-cell communications. We considered the six cytokines communicating nodes in our multilayer network model. For both a single-layer CD4+ T cell network and a multilayer CD4+ T-cell network, we randomly chose 1000 different initial states and then examined their state transitions during 2 400 = 800 time steps. Changing the parameters and , we computed .
For the simulation, the following parameters were used: The specific values of parameters follow Table 1. Our simulator was implemented in Java.

Antifragility in Multilayer RBNs
Figure 3(a) shows average fragility of multilayer RBNs for = 1, 2, 3, 4 depending on perturbation frequency when perturbed node size = 80. As grew (i.e., the period of adding perturbations became longer and longer), the antifragile or fragile dynamics of the ordered, critical and chaotic networks changed into robust beyond = 30 although about half nodes were perturbed ( = 80). This result means that how often perturbations are added has a greater impact on antifragility than how many nodes are perturbed. Figure 3(b) represents average fragility of multilayer RBNs depending on when = 1. For all , there were certain ranges of for the networks to be antifragile. We also found that the certain ranges of decreased and antifragility declined as increased. These findings indicate that a moderate level of perturbations can induce "optimal" antifragility, and multilayer RBNs take bigger benefits from perturbations in the order of ordered, critical, and chaotic networks. From Figure 3(a) and 3(b), we can see that maximal antifragility is obtained with a moderate level of perturbations. It is worth noting that the maximum antifragility varies for different values. Figure 4 explains the reason why multilayer RBNs obtained more improved antifragility in the order of ordered, critical, chaotic networks. In Figure 4(a), the complexity before adding perturbations gradually increased as got bigger, while in Figure 4(b) the complexity after adding perturbations was increased as got smaller except for the early range of (1 20). Figure 4(c) shows the difference of complexity before and after perturbations. As seen in the figure, the smaller was, the larger the difference was. It means that the complexity representing the balance between regularity and change can be improved by perturbations, and the degree of the improvement is much larger for smaller .  Regarding the complexity, one interesting finding is that the complexity before perturbations (Figure 4(a)) is not consistent with that of previous studies [41,45,46]. The existing studies demonstrated that critical RBNs have higher complexity, while our result revealed that chaotic multilayer RBNs have larger complexity. The difference is due to a cellular level. The previous studies were performed at a single cell level, and our study was conducted at a multicellular level. The result that multilayer RBNs and single-layer RBNs exhibit different dynamics emphasizes the need for research in the context of multicellular settings. Figure 5 shows average fragility of multilayer RBNs depending on two parameters related to the multilayer structure: the number of communicating genes and the number of links of an intercellular network. In Figure 5(a), as the number of communicating genes increased, multilayer RBNs in all the regimes became antifragile. In Figure 5(b), the fragility values did not change significantly except for the early range as the number of links of an intercellular network increased. These results indicate that the communicating genes have a larger effect on antifragility at a multicellular level. Also, it suggests the possibility that the number of communicating genes representing the degree of interactions between cells might be able to be used as an indicator to estimate the effect of multilayer structure on antifragility.

Comparison of Antifragility Between Multilayer and Single-layer RBNs
To study how the multilayer structure has an effect on the production of antifragile networks, we calculated probability of how many antifragile networks were generated in multilayer and single-layer RBNs, respectively and then compared them. Figure 6 shows heat maps representing the probability values in a diverse range of and for multilayer and singlelayer networks.  For all the cases of (a), (b), and (c) in Figure 6, we found that they had similar trends that the smaller was, the more frequently antifragile networks were produced. It means that the trends of antifragile dynamics at a single-cell level are still maintained at a multicellular level. Also, as the perturbed node size increased and the period of adding perturbations became shorter, the probability of generating antifragile networks decreased overall. It is worth noticing that, especially for large values, there is not much difference between the singlelayered RBNs independently of their size (Figures 6(b) and 6(c)). On the other hand, multilayered RBNs are clearly more antifragile (Figure 6(a)). Figure 7 shows the comparison of the probability acquired from Figure 6 based on statistical tests. Firstly, to investigate the difference of antifragility between multicellular and singlecell systems with the same system size, we compared multilayer RBNs with 162 nodes and single-layer RBNs with 162 nodes. As seen in the figure, the multilayer networks produced antifragile networks more frequently at = 2, 3, 4. In the case of = 1, the single-layer networks had a higher probability, but the values were not so different. Secondly, to examine the difference of antifragility between a multicellular system and its component, we compared multilayer RBNs with 162 nodes and single-layer RBNs with 18 nodes. We found that the multilayer networks generated antifragile networks with higher probability for all . The findings from Figure 7 indicate that the multilayer structure helps to produce the greater number of antifragile networks, especially for larger values.

Comparison of Antifragility Between Multilayer and Single-layer CD4+ T-cell networks
Using a CD4+ T-cell network related to the immune system as an example of biological systems, we calculated average fragility and the probability of generating antifragile networks for an individual CD4+ T-cell network and coupled CD4+ T-cell networks. As seen in Figure 8(a) and 8(c), the general tendency of antifragile dynamics in the CD4+ T-cell network at a multicellular level was practically the same as the tendency at a single-cell level. We compared multilayer and single-layer CD4+T-cell networks based on the probability values in Figure 8(b) and 8(d). We found that multilayer CD4+ T-cell networks generated antifragile networks with a little higher probability on statistical tests (i.e., multicellular level: 99.26% single-cell level: 97.74% under p-value 0.05). From a biological viewpoint, these results suggest that the properties of biological systems might be enhanced in the structure of interacting multiple subsystems.
When compared to multilayer and single-layer RBNs, the CD4+ T-cell network showed similar antifragile dynamics to the dynamics of multilayer and single-layer networks at = 1. From this, we can infer that the CD4+ T-cell network may be ordered. This can be understood because immune cells probably not only have a variable environment, but actually have evolved to thrive on it. The finding on the ordered dynamics of the CD4+ Tcell network is consistent with the results of many studies demonstrating living organisms are ordered or critical [47][48][49][50].

Conclusions
In this study, applying our (anti)fragility measure to multilayer and single-layer BNs, we studied how the dynamics of the networks are varied depending on relevant parameters, and how the multilayer structure affects the antifragility of the whole system. We found that systems showed different dynamics depending on the degree of perturbations and the degree of interaction between system components: fragile, robust or antifragile. Also, we found that the multilayer structure facilitated the production of antifragile systems. Probably this is related to the modular structure of multilayer RBNs [33], although further studies should be made.
The findings can be utilized for various applications such as systems biology and bio-inspired engineering. We can control system properties from fragile through robust to antifragile dynamics by adjusting the size and frequency of perturbations, and the number of communicating nodes. Our results may also help to better understand dynamical behaviors of multicellular organisms. Besides, using the multilayer structure, we could create engineered systems with an increased antifragility.
Our study has a few limitations. Firstly, our multilayer RBN model is the one where identical RBNs are randomly coupled. However, there are many other systems where different subsystems are connected to each other and they communicate in a certain way. Secondly, we used only one biological example in explaining the dynamical behaviors of biological systems. To obtain more generalized findings on antifragility, we plan to develop different kinds of multilayer network models, explore them, and use various biological systems.

Data Availability
Our source code and data are available at https://github.com/bin20005/antifragility

Conflicts of Interest
The authors declare that there is no conflict of interest regarding the publication of this paper.

Funding Statement
This research was partially supported by CONACYT and DGAPA, UNAM.