Enhanced robustness of evolving open systems by the bidirectionality of interactions between elements

Living organisms, ecosystems, and social systems are examples of complex systems in which robustness against inclusion of new elements is an essential feature. A recently proposed simple model has revealed a general mechanism by which such systems can become robust against inclusion of elements with totally random interactions when the elements have a moderate number of links. The interaction is, however, in many systems often intrinsically bidirectional like for mutual symbiosis and competition in ecology. This study reports the strong reinforcement effect of the bidirectionality of the interactions on the robustness of evolving systems. We show that the system with purely bidirectional interactions can grow with twofold average degree, in comparison with the purely unidirectional system. This drastic shift of the transition point comes from the reinforcement of each node, not from a change in structure of the emergent system. For systems with partially bidirectional interactions we find that the regime of the growing phase gets expanded. In the dense interaction regime, there exists an optimum proportion of bidirectional interactions for the growth rate at around 1/3. In the sparsely connected systems, small but finite fraction of bidirectional links can change the system’s behaviour from non-growing to growing.

each extinction event will also modify the fitness of the other species, the fitness is re-calculated and the least-fit species is re-identified. The deletion of species is continued until the minimum fitness value of a species becomes positive, meaning that nothing more will happen in terms of the above process. After finding such a state, we proceed to the next time step by adding a new species into the system. The m interacting species are chosen randomly from the resident species with equal probability and the directions are also determined randomly. The link weights are again assigned randomly using the standard normal distribution. Then, we re-calculate the fitness of each species to find the species that will become extinct.
The results of the model have shown that the repetition of this process gives rise either to continuous growth or stagnation for the system size, depending on the model's unique parameter, namely the number of weighted unidirectional links per a new node, m. Furthermore it turns out that the system can grow only if the connections are moderately sparse, i.e. the average degree is within a rather narrow interval for 5 ≤ m ≤ 18. If the network becomes denser, there is a transition from growth to decay. It was shown that this transition originates from a balance between two effects, namely while on one hand the inclusion of species with more interactions makes each node more robust, on the other hand it also increases the impact of the loss of a node 14,15 . This relation might be the origin of moderately sparse network structures ubiquitously found in real world cases, with the average degree of up to few tens [16][17][18][19] .
We should note that the approach of simple models we adopt here, for obtaining a new way to determine the system's robustness, is quite different from the more standard and often used models of population dynamics. In our simple model, all the intrinsic dynamics of the systems is indirectly represented by the survival or fitness condition, i.e. f i > 0, which does not depend on the population of the species, as we have not introduced the population size as the degree of freedom. With the focus on robustness of the system while keeping the model general and "as simple as possible but not simpler", we have ignored such details as the intensity of the flux of energy or nutrition (mass) in ecosystems, rates of chemical reactions, and trade volume in a company network, etc. Our aim with this general and simplified modelling approach is to study in addition to ecological systems even social and economic systems, in which it is hard to find precise equation of motion or evolution rule.
As a further support for this kind of universality seeking approach we mention that in several cases such modelling has been very successful. For example, one can find good agreement between a simple model and a real ecological system in the statistics such as the distribution of the lifetime of species 20 . Further examples of the power of a simplified modelling approach comes from looking neuronal and gene regulatory networks that have quite different types of elements considered as nodes, yet the model dynamics is the same for both of them with elements getting a shunt or off state going extinct and never getting excited. As a last justification worth mentioning is to regard the survival condition in our model as the condition for persistence in an original complex dynamical model. An interesting finding is that we can show a relatively natural example that the extinction condition in a population dynamics model with ratio-dependent predation term, proposed and studied in ecology to improve the model's correspondence with reality 21 , reduces to the same population-independent form as is the case in our simple model 22 .
In the previous yet even simpler model, the interactions between species of the system were generated randomly with certain probability as unidirectional links between a pair of species i and j. This means that a random unidirectional link from j to i with a unidirectional link from i to j existing already, could be generated as a second order effect proportionally to its link density, as ~ρ 2 . Hence the existence of influence from species i to j is independent of the existence of an influence in the opposite direction, i.e. from j to i. However, the simulation results of this model of self-emergence have shown that the accumulation of bidirectional influence is negligible, especially for large systems their density turned out to be vanishingly small.
However, in real systems, the interactions can on one hand be considered unidirectional as in case of reaction dynamics or signaling, and on the other hand bidirectional reflecting various types of bidirectional influence, like symbiosis or competition between species, predation-prey ecology, and even action-reaction law of Newtonian mechanics, etc. Therefore, in real systems we could distinguish five main types of interactions, of which two are unidirectional in nature, and three reflect bidirectional influence in terms of being cooperative, predatory, and competitive, as shown in Table 1 with the characterization and some real life examples. Our aim in this study is to gain better understanding of the robustness and growth of real systems, with computational modeling. In this study we will focus on studying systems in which the interactions around the species newly introduced into a system can be unidirectional, bidirectional, or a mixture of both in certain proportions. In the Model section we revisit the previous model and extend it to include bidirectional influence between species. In the following Results section we describe the simulation results for the cases purely and partially bidirectional interactions. Next in the Discussion section we summarize our findings. This is then followed with the extensive Methods section.

Model
As in this study we are interested in the effects of bidirectional influences on the robustness of the evolving open system, we take the previous model as the base describing it first in detail and then extend it by introducing bi-directionality in a controlled way. In the previously model the evolution of an open system of species was set up to take place in a graph or network of randomly connected species or links describing the influences between species. This dynamical system evolves through interaction dependent survival of each species and a successive addition of new species to the system, as depicted in Fig. 1. This network consists of nodes (species) that are sparsely connected to each other with directed and (positively or negatively) weighted links, each describing the strength of interaction between a pair of species. In this system the survival of each species is determined by its fitness, which is calculated by the sum of the weights of its incoming links such that if the fitness of a given species is positive, it will survive. On the other hand if the fitness of a species is non-positive that species and all the interactions from and to it are deleted. If two or more species have non-positive fitness values, the species with the minimum fitness are deleted and the fitnesses of other species are recalculate. This so called extinction process is repeated until all the remaining species have positive fitness and once such a network structure is reached, nothing further will happen until a new node is born or introduced into the system. According to the terminology of ecology this is called a persistent state rather than a stable state.
When the system has relaxed to such a persistent state, a new species (node) is added to the system, and the above described process continues. Next M unidirectional links between the new node and the nodes randomly selected from the resident nodes are introduced.
Till this point our extended model is the same as the previous model, but instead of introducing unidirectional links a mixture in total of M unidirectional and bidirectional links is introduced by generating the bidirectional links with probability P and unidirectional links with probability 1 − P. The weight for each link is chosen again from a normal distribution. Thus our extended model has only two control parameters, M and P, instead of one in the previous model. In order to run the dynamics of both these models the above described addition and survival check process is repeated. A pseudo-code like description of the entire procedure of the present model is shown in the Methods section.
Counting each bidirectional link as two unidirectional links, one can say that m = (1 + P)M unidirectional links to the new species is introduced, which on average constitute PM bidirectional links. The extreme cases of P = 0 and P = 1 correspond to the original unidirectional model (with the number of unidirectional links per new species, m = M) and the pure bidirectional model (with m = 2M), respectively. We denote the density and the amount of bidirectionality of interactions in the emergent system in the same way, by the average number of interactions per species, i.e. K = ∑ i M i /N, where M i and N represent the number of links of i th resident and total number of resident species, respectively. We denote the ratio of bidirectional links by Q. Then the resident species The introduction of a new species (red). In this case the inclusion leads to the extinctions of a resident species (yellow) and that extinction triggers another extinction of a species, which is not directly connected to the newly introduced species (magenta). (c) Finally, the system reaches to a new persistent state. The numbers in the square brackets demonstrate the control parameters; M: the number of interacting species of the newly introduced species, P: the ratio of the bidirectional relation in the initially assigned interactions, and m = (1 + P)M: unidirectional link degree of the newly introduced species, and observables of the emerging network; N: the number of species, K: average number of interacting pairs per species in the system, Q: the ratio of the bidirectionally interacting pairs to the all interacting pairs in the system, and k = (1 + Q)K: average unidirectional link degree of the system. has k = (1 + Q)K unidirectional links such that they form QK bidirectional links. The above definitions are summarized in Table 2.

Results
Purely bidirectional interaction (P = 1). Transition in growth behavior at twofold critical point. We first examine the system with purely bidirectional interaction, i.e. we choose P = 1. Similarly with what was obtained for the model with unidirectional interactions, the purely bidirectional model yields two distinct phases characterized by its growth behavior. This is dependent on the unique parameter that specifies the number of interactions, M, we introduce to each new species ( Fig. 2(a)). While the system keeps growing without any limitation to its size if M is smaller than a certain value We call these phases diverging phase and finite phase, respectively. From the systematic simulations, we have estimated the critical number of bidirectional interactions to be at between M = 19 and 20. As in this case all the interactions are bidirectional, the critical point, measured by the number of unidirectional links, is ). This value is almost two times that of the pure unidirectional model (i.e. = .
= m 18 5 ). Another clear difference from the unidirectional model is the absence of a finite phase in the very sparse regime. In the unidirectional case, the systems with m ≤ 4 are in finite phase, due to the fragility of the emerging network. In contrast, in the bidirectional case, the system even with M = 1 turns out to exhibit diverging behavior. This is basically caused by the easiness of forming mutually supporting "dimers", which is enough to let the nodes survive even when they are not globally connected. This regime will be discussed later in more detail. Our primary aim in what follows is to understand the twofold transition point in the well-connected regime.
To illustrate why we measure the interaction density by the number of unidirectional links and call the transition point "twofold", we will briefly review the mean field argument on the basis of which the understanding and the rough prediction of the transition in the unidirectional model was obtained 14,23 . In a system the introduction of a new species can lead to the extinctions of the resident species. We denote the extinction probabilities during this event as E i . It should be noted that the extinction of resident species can also cause a cascading extinction through the deletion of the links from it, with a certain probability E e . Assuming a structure-less random network topology, which is supported by the direct observation of the emergent network, the average number of extinctions N E per a newly introduced species can be calculated as a simple infinite geometric series as follows for the unidirectional model, and, Reinforcement of elements. In order to clarify the crucial factor for the shift of the transition point, we evaluate the transition point using Eq (1) and Eq (2) with empirically obtained extinction probabilities (see the Method section). The average number of total extinctions per a newly introduced species increases monotonically as a function of the degree of a new species both for the unidirectional model and for the bidirectional model. The quantities = N E P 0 and = N E P 1 cross the critical value 1 at m * ~ 19.5 for the unidirectional model and 2M * ~ 39 for the bidirectional model, respectively, as depicted in Fig. 3(a). These critical points m * and 2M * agree well with the real transition points, = .
= m 18 5 . This result suggests that the shift of the transition point stems from the decrease in the extinction probability of the emergent system. In Fig. 3(b) we see that the quantities E i and E e decrease monotonically as a function of the degree of a new species both for the unidirectional and bidirectional models. Comparing the bidirectional model with the unidirectional model, however, the bidirectional model gives a smaller E i and E e if the new Further investigation of the emerging system reveals that the larger decrease in the extinction probability of the bidirectional model is caused by the slight increase in the average degree, together with the small change in the link weight correlation (see Method section for more details).

Model with partially bidirectional interaction (P ∈ [0, 1]). Bidirectionality dependence of the transition
of the growth behavior. Next, we examine the systems with partially bidirectional interactions to investigate the influence of the proportion of bidirectionality P on the growth behavior. We find that P ∈ (0, 1) expands the area of the diverging phase, as seen in Fig. 4(a). Interestingly or even surprisingly, the system with a certain level of bidirectionality steadily grows even if the system with purely unidirectional interactions (P = 0) or purely bidirectional interactions (P = 1) is in the finite phase. There are two phase boundaries, one found at around 18 < M < 22 and another at around 1 ≤ M < 5, as can be seen in Fig. 4(a). Here the finite phase remains not only in the dense interaction regime and but also in the sparse interaction regime.
One may expect a strong effect of birectionality, because strong accumulation and cooperative interaction is possible. The most prominent example of this is the case of M = 1, in which each cooperative dimer can survive independently letting the system to consist of many independent clusters with broad size distribution. However, in the dense regime, one cannot find such a strong structuring effect. For example, in the dense regime, the proportion of bidirectional links in the emerging system Q remains almost the same as the initial proportion P (see Methods). Growth behavior in the dense interaction regime (19 ≤ M ≤ 22). Figure 5(a) shows the speed of divergence in the dense interaction regime. Here the speed of divergence is asymmetrically distributed as a function of a proportion of the bidirectional links. It is seen that there is an optimum bidirectionality for the growth rate P ~ 1/3. Thus, the upper phase boundary in the denser interaction regime (18 < M < 22) is asymmetric for the proportion of the bidirectional links ( Fig. 4(b)). We also find that there exists an optimum bidirectionality for the growth rate at around P ~ 1/3. The system with a moderate amount of bidirectionality grows for M = 19, 20, and 21, even if the system with purely unidirectional interactions or purely bidirectional interactions (P = 0 or P = 1) is in the finite phase. In addition, a small amount of bidirectionality makes the partially bidirectional system grow. On the other hand, the system with M = 22 remains in the finite phase only for all P ∈ [0, 1].
Growth behavior in the sparse interaction regime (1 ≤ M ≤ 4). In Fig. 5(b) it is seen that the speed of divergence increases monotonically from zero with the increasing proportion of the bidirectional links. Thus, a small but finite amount of bidirectionality brings the system from the original finite phase to the diverging phase, as is evident in Fig. 4(c). The area of the finite phase gradually decreases as a function of the amount of bidirectionality. Then if the proportion P of the bidirectionality is set larger than a certain value, the sparsely connected system involves only a diverging phase. We can also observe that, in this sparse regime, the emergent systems tend to have  Table 3 in the Method section). This indicates relatively strong self-organization that is regarded as the main reason for the change in the system's growth behavior such as the disappearance of the finite phase.

Discussion
We have studied the effects of introducing new species with variable degree of bidirectionality in the interactions between species, focusing especially on the dynamics of the robustness of the network of species.
First we found that the system with purely bidirectional interactions is more robust against the inclusion of new species, which is characterized by the twofold transition point from the growing phase to the finite phase. The twofold transition point also means that one can roughly identify or even mix the unidirectional and bidirectional  interactions to apply our framework to real systems, in which the detailed information of interactions and their nature is often not available. The upper limit of the average number of interacting species per species for letting the system to grow its diversity is around 20 irrespective of the bidirectionality of interaction. This is again confirmed to be significantly looser than the classical diversity-stability relation known for random dynamical systems 1 and for generalized Lotka-Volterra system with random inclusion of new species 24,25 . This limiting condition seems to be more consistent with the persistence condition of ecological models with realistic nonlinear interactions 8 and of a model with spatial structure 26 , but most importantly with empirical data 19 .
In the model with partial bidirectionality, we find that there exists an optimum ratio of bidirectional interactions at around 1/3, for the growth of the system to take place. Because the bidirectionality of the emergent system remains near to its input proportion, this optimum in the middle seems to be consistent with the enhancement of the system by the diversity of interaction types in the system, found in a more complex models of ecosystems 27 . We also find that only a small, though finite, portion of bidirectional interactions let the very sparse (M ≤ 4) system to grow to a large and connected structure.
Except for the very sparse regime, all those drastic changes in the entire system's robustness are caused by a little change in the degree and also in the emergent correlations between the link weights. The network characteristics such as clustering and assortativity are also found to be not much deviated from that of Erdös-Rényi random graph (see Methods). This is quite noteworthy because the system in our model self-emerges without any apparent constraints and hence more drastic self-organization could be possible. Indeed, stronger self-organization and its effect on the growth behavior of the system is observed in the sparse regime (M ≤ 4).
A study of an extended model in which we force the system to have stronger correlations in the link and weight topologies could be the next step 6,26,28,29 , to further relate our framework to the real systems. (iii) Go back to the step 1 (to re-evaluate the fitnesses). 3. A new species is added to the system ( Fig. 1(b)).  Table 3. Network characteristics of the emergent systems. For the case in the diverging phase, we observe the system when it first reaches a certain size N obs = 20,000. The numbers in the parenthesis for clustering coefficient and nestedness are those ratios to the ones of Erdös-Rényi random graph with the same size and degree, i.e. N = 〈N〉, K = 〈K〉, Q = 〈Q〉.
(ii) The new species forms bidirectional links, with probability P.
(iii) With probability 1 − P, we assign a unidirectional link. Its direction is randomly chosen with probability 1/2. (iv) The link weights are independently drawn from the standard normal distribution. The link weights between incoming and outgoing unidirectional links which consist up one bidirectional link are also independent: 〈a ij a ji 〉 = 0. (v) Go to 1.
The mean-field estimation of the transition point. In the original unidirectional model (P = 0), the system keeps growing only if the interaction is moderately sparse (5 ≤ m ≤ 18) and cannot grow outside of this regime (m ≤ 4 and m ≥ 19). The basic mechanism of the transition from the diverging phase to the finite phase at m c = 18.5 is known to be explained by the mean-field argument 14,15 , which assumes a correlation-less random network structure for the interaction network of self-emerging system. The validity of this assumption is confirmed in the simulation, from the direct observation of the system. With such an assumption, the system's growth behaviour against the inclusions of new species is characterized by two average extinction probabilities, E i and E e . E i is the extinction probability of the resident species after obtaining a direct incoming link from the newly introduced species. The other extinction probability, E e , is the probability of a species to get extinct after loosing one incoming link which is brought by the extinction of the neighboring species. Because of the assumption of random network topology, the total number of extinctions of the resident species per newly introduced species 〈N e 〉 is calculated from a simple geometric sum as Let us perform the same calculation in the system which consists of fully bidirectional interactions (P = 1). In this case, the expected number of extinctions in the layer of direct connection to the newly introduced species is ME i . These extinctions cause the following link-deletion events, subsequently giving rise to KE e extinctions in the next-nearest neighbors. After this layer, the contribution of the cascading extinction in each deeper layer is (K − 1)E e . Therefore, the total number of extinctions per a newly introduced species in the bidirectional model is calculated as, For the simple evaluation of the transition point, we further assume that the average degree of the emergent system is same as the initial degree, i.e. m = k and M = K. Because the correlation-less network assumption it also means that the extinction probability is a function of the unique system parameter, i.e. average degree, meaning that it can be written as a function of initial degree m = 2M. Based on the empirical findings, we also assume that the extinction probabilities are similar and much smaller than 1 near the transition point . Then the condition for the unidirectional case (Eq. 3) becomes the following very simple form 14 c c and the corresponding simple form for the bidirectional system (Eq. 4) is These results gives us the basic estimation for the transition point for the bidirectional systems, which says that the transition occurs at the same critical unidirectional degree: This is why we regard the observed transition point, Calculation of the extinction probabilities in the emerging systems. The extinction probabilities E i and E e in an emerging system are calculated according to those definitions as, where ≤ # a f ij i and # link represent the number of species with smaller fitness than the link weight receiving from each link and the total number of links of the emergent system, respectively.
The network structure and its effect in the bidirectional model (P = 1). Degree distribution. As shown in the Fig. 6(a), the average unidirectional link degree of the emergent system k increases almost linearly as a function of the degree of a new species m, both for the unidirectional model and for the bidirectional model. For the unidirectional model and the bidirectional model with the same degree of a newly incoming species m P=0 = 2M P=1 , the average degree of the bidirectional model is slightly larger than that of the unidirectional model. This difference comes mainly from the broader tail of degree distribution in the bidirectional case. The peak position of the degree distribution at k = m and the tail on the large degree side is fatter than that on the other side ( Fig. 6(b)).
One can find a scaling relation in the fitness distributions of the unidirectional model and the bidirectional model ( Fig. 7(a)). The fitness distributions of the unidirectional model and the bidirectional model match well when the control parameter of those satisfies the relation: m P=0 = 2(M P=1 + 1). This scaling relation holds also for the distributions of the link weights ( Fig. 7(b)). All the results above can be explained by the slight increase in the actual average degree in the emergent system in the bidirectional case: k P=1 (M) = k P=0 (2(M + 1)).
Weak correlation in link weights. According to the scaling relation, the unidirectional model and the bidirectional model give the same average degree of the emergent systems k P=0 = 2K P=1 at m P=0 = 2(M P=1 + 1). Therefore, the unidirectional model and the bidirectional model are expected to have the same extinction probabilities at k P=0 = 2K P=1 (Fig. 8(a)). However, while E i is found to obey this relation, E e of the bidirectional model is slightly smaller than that of the unidirectional model ( Fig. 8(b)). It suggests that the system of the bidirectional model has stronger weight-fitness correlation than the unidirectional model.
Here we conclude that the shift of the transition point stems from the small decrease in the extinction probability of the bidirectional model due to the weak shift in the average degree and the correlation in link weights, not from strong structuring or emergence of a certain motif in the emergent system.
The network structure of the partially bidirectional model. Dense regime. The self-emerging network of the partially bidirectional model in the dense regime is found to have only weak structuring, similar to the fully bidirectional case. The proportion of bidirectional links in the system also remains the same as the newly introduced interactions (Q = P).
Sparse regime. In the sparse regime, the emergent systems tend to have larger degree and larger proportion of bidirectional links than input values (〈K〉 > M and 〈Q〉 > P, see Table 3). This indicates the relatively strong self-organization effect in this regime, which is regarded as the main reason for the disappearance of the finite phase.
In order to investigate the network structure in this regime, we analyze the network size distribution. As shown in Fig. 9(a), the systems are composed of multiple clusters in the diverging phase, whereas all resident species construct a single big cluster in the finite phase. In the diverging phase, in contrast, the systems with M = 2, 3 and 4 are composed of a single big cluster with some small fragments such as dimers, trimers, etc. (Fig. 9(b)). The cluster size distribution of the system with M = 1 is continuous up to its largest cluster size, meaning that there is no dominant giant component. The basic mechanism of the self-organization seems to be a formation of cooperative dimers. Such a dimer can survive independently and hence can be a nuclei of non-connecting cluster. We can see this effect in the cluster size distribution clearly in the case of M = 1. Although it becomes less visible in the cluster size distribution in the case of M > 1, because of the random reconnection (or, bridging) process due to the introduction of new nodes, the dimer nucleus still let the network, which consists of sparsely connected almost independent clusters, survive. Assortativity, Nestedness, and Clustering. Some network characteristics of the emergent systems are shown in Table 3. For this we treat unidirectional and bidirectional links equally as directionless links. Therefore the assortativity coefficient is calculated from the directionless link degrees of the nodes in the both ends of link i, κ α i and κ β i , as follows where a ij and K i represent the adjacency matrix element and the degree of species i, defined by the directionless links. The clustering coefficient is also calculated from the number of triples and triangles formed by the directionless links as,  Extinction probability E e . Emergent systems of the unidirectional and bidirectional models with same average unidirectional link degree show the same extinction probability against the direct disturbance event from the new species, E i . On the other hand, the extinction probability against the loss of a link, E e , is slightly smaller for the bidirectional model than that of the unidirectional model.