Lévy Walk in Swarm Models Based on Bayesian and Inverse Bayesian Inference

Graphical abstract


Introduction
Living systems are adapted to open dynamic environments and can overcome various problems arising in environments [1,2]. In other words, they can compute what the environments require and can resolve them to some extent. If the problem is well defined, optimization techniques can be applied to the problem. However, the problems living systems face can be ill-defined and can change over time. Notwithstanding these situations, living systems can reach quasi-optimal solutions [3,4]. While various attempts using bioinspired computing have proposed defining this adaptability [5][6][7][8][9][10], the issue of the essential property of adaptability has not yet been found.
A quasi-optimal solution in open environments could be achieved by balancing sticking to a specific well-defined problem with giving up that problem, since open environments could change the problems themselves over time. This can be replaced by balancing highly efficient computation with universal computation [1,11], balancing exploitation with exploration [12][13][14][15] and balancing specialist with generalist strategies in adaptation [16][17][18]. In terms of dynamical systems, balancing the exploitation at an attractor basin with the exploration of wandering various attractors is called the edge of chaos, since staying at an attractor shows oscillation or ordered patterns and wandering to various attractors shows chaotic patterns [19][20][21][22][23][24][25]. It is known that systems tuned at the edge of chaos or at the critical point can have computability, balancing universal computation with highly efficient computation [21][22][23][24][25]. The critical state is characterized by a power law distribution [26][27][28][29]. What can make a system be tuned to the critical state? Although the idea of self-organizing criticality is one of the candidates [30][31][32], no generalized method of tuning at the critical state has been found yet.
Swarming behavior can be a touchstone to determine the core of quasi-optimal solutions in open environments. Recently, many optimization techniques have been developed based on swarming behaviors [33][34][35][36]. A swarm is neither a machine-like order nor just disorder, and a swarm shows a critical state between order and disorder. Since the Self-Propelled Particle (SPP) model is one of the most powerful models for swarms, swarming behavior can be regarded as a critical phenomenon in the phase transition between order and disorder [37,38]. The SPP model is based only on velocity matching, although other models implement collision avoidance and flock centering [39][40][41][42]. Because the phase transition is controlled by the parameter, the critical state is not selforganizing, and it requires parameter tuning. Indeed, real animal swarms show various properties of the critical state, such as scale-free correlation and Lévy walk, which are characterized by a power law distribution [43][44][45][46][47][48][49]. The problem of what makes a system tuned to the critical state is an open problem.
Here, we show that a critical state balancing exploration (disorder) and exploitation (order) can be obtained through Bayesian inference [50][51][52][53] coupled with the specific interaction between data and hypotheses [54][55][56]. We take the SPP model as the basic swarm model and introduce an agent that can infer the state by using Bayesian inference. The question becomes how to obtain the critical behavior of the SPP. The critical behavior can be estimated by a power law distribution of the individual walk pattern, called the Lévy walk. Although there are many approaches to explain the Lévy walk, they are nonsystematic approaches [57][58][59]. Swarm approaches based on self-organized criticality [60][61][62] have universally failed to explain the Lévy walk [63][64].
We implement a swarm agent as the decision maker using Bayesian inference coupled with the data-hypothesis interaction. This interaction was previously proposed by us and is called inverse Bayesian inference [54][55][56]. Although hypotheses are not altered through Bayesian inference, both the probability and the likelihood of each hypothesis are constantly changed in our model. Normal Bayesian inference does not reveal critical behavior, but our inference system succeeds in explaining the critical behavior or Lévy walk. Real animal swarms, schools and/or flocks exhibit spatial and temporal changes in their collective behaviors, including translation, splash and tornado behaviors. Although detailed parameter tuning is required to reveal these special behaviors in previous swarm models [65][66][67][68], our model, aimed at showing this critical behavior, easily reveals the coexistence of translation, splash and tornado behaviors in a swarm. In the context of the edge of chaos, disturbed patterns and/or disorder are called ''chaos" independent of the definition of chaotic dynamics. In our paper, we use the term chaos in this sense, especially disorder. Since our inference system relying on data-hypothesis interaction easily and ubiquitously shows critical behaviors, the system exhibits universal criticality. Finally, we examine the significance of universal criticality in realizing the optimal design in an open environment for both natural and artificial design purposes.

SPP with Bayesian and inverse Bayesian inference
The collective behavior in swarms, flocks and schools could contain the various intrinsic dilemma between social norms and individual free decision [69,70]. One of the simplest models, the Self-Propelled Particle (SPP), implements the social norm by velocity matching of an individual's neighborhood and their freedom by fluctuation [37,38]. In our model, the kth agent at tth time step in a swarm calculates the average velocity in the form of an angle, h t k , for any jth agent satisfying the following condition: where x t k ; y t k À Á represents the location of the kth agent at the tth time step, and R represents the radius of the neighborhood. The (t +1)th location of the kth agent is determined by: where rnd e ð Þ is the perturbation randomly generated in 0.0 rnd e ð Þ e, and d t k 2 D ¼ f0; 1; 2; 3g is the angle against the velocity matching inferred through Bayesian and inverse Bayesian inference mentioned later. V is the unit velocity in the form of scalar. It is clear that if d t k ¼ 0, the kth agent accepts the velocity match under a small perturbation, and if d t k ¼ 2, the kth agent escapes from its flock mates in the reverse direction. While d=0 reveals ballistic behavior showing exploration, d=1, 2 and 3 reveal perturbed behavior showing exploration. Thus, these data reveal variation of exploitation and exploration, and the probability distribution of data constitutes a hypothesis. If d t k ¼ 0 for any d and t, then the swarm model is defined as simple SPP or simply SPP.

Bayes Only (BO)
The second Bayes-Only (BO) model is defined by the SPP model equipped with Bayesian inference. Each agent k makes d t k by using Bayesian inference. The inference process of the kth agent proceeds as follows. In regard to d t k in equation (3), the conditional probability of any hypothesis h 2 H ¼ f0; 1; 2; 3g under data d t k is expressed as: for any d 2 D can be chosen with the corresponding probability. The procedure based on the probability is the same as that in equations (7)(8)(9). The cumulative probability conforming to: is defined, and random variable 0:0 p 1:0 is then updated. Bayesian inference consists of the above procedure equations (4)- (11).
As mentioned before, in Bayesian inference, a set of hypotheses is not changed through the inference process. Since a hypothesis is defined by its likelihood such as PðdjhÞ, the distribution of PðdjhÞ is not altered in Bayesian inference. In contrast, here, we introduce the interaction between data and hypotheses below. If the swarm model is implemented by equations (1)- (11) and P t k ðdjhÞ is invariant over time, we call the model an SPP with BO model, or simply a BO model.
where #S for set S is the number of elements in S, and #Dat ¼ m The interaction between data and hypotheses or the inverse Bayesian inference is defined by: f tþ1 k is determined by the probability, for any given random variable 0:0 r 3:0, and f tþ1 k is obtained by: It is evident that P tþ1 (5) is symmetric to P tþ1 k djf tþ1 k ¼ P t k ðdÞ in equation (14), which is why the procedure of equation (14) is called inverse Bayesian inference. On the one hand, Bayesian inference contracts the condition of an event (hypothesis) by replacing the probability of P tþ1 On the other hand, inverse Bayesian inference extends the condition of an event (data) by replacing P tþ1 k djf tþ1 k with P t k ðdÞ. Thus, Bayesian inference balances contraction (exploitation) with extension (exploration). This is considered in a later section and compared to the mechanism of self-organized criticality.
In our simulation studies, the program for the SPP in which Bayesian and inverse Bayesian inference (BIB) are implemented is shown in Figure 1. In the program, if d t k ¼ 0 for any t, the program simulates the simple SPP under which each agent obeys the velocity matching expressed by equations (1)-(3) without any inference. If certain sourced codes implying equations (12)- (17) are commented out, the program simulates the SPP containing only Bayesian inference. This implies the likelihood of hypotheses P t k djh ð Þ ¼ P tÀ1 k djh ð Þ for any d and h. It is easy to compare the simple SPP, the SPP with only Bayesian inference (BO), and SPP with BIB.
The SPP part of our model contains various parameters, the radius of the neighborhood, R, the unit velocity of each agent, V, and the perturbation (noise) for the velocity matching, e (radian).
Throughout all simulation studies, we set R=20.0 and V =5.0 if there is no description. The perturbation e varies and corresponds to the phase parameter controlling the phase shift.

Step length distribution for the power law analysis
To estimate the critical behavior in the model swarm, the distribution of the step length is measured for the simulation studies. In research on animal foraging, the step length distribution for various animals has been measured and analyzed with respect to the power law distribution [71][72][73], since animal foraging always encounters exploitation and exploration dilemmas. If animals implement the exploitation strategy and consume food resources in a closed environment, the walking patterns reveal random walks. In contrast, if animals implement the exploration strategy to search for other food resources and leave their previous environments, the animals reveal ballistic walking trajectories, which could realize efficient searches for unknown resources. If the step length is defined by the distance between two bending points, the ballistic walk implies a very long step length, while the random walk implies a short step length normally distributed around the mean step length. Animals could balance exploitation with exploration, and the walk pattern would reveal a random walk pattern with a long tail, which is characterized by a power law distribution of the step length. If the exponent of the power law distribution ranges from 1.0 and 3.0, the walk pattern is called the Lévy walk. It is known that the foraging patterns of animals frequently indicate the Lévy walk [43][44][45][46][47].
Here, we define the step length by the following. Since the location of an agent is moved in a stepwise fashion, we define the bending angle a in the agent's walk as: is the location of the kth agent at the tth time step.
In the walking pattern of agents, if a t > a max , one step walk is then terminated, and the distance, D, between the previous bending point ðX k ; Y k Þ and the new bending point x tÀ1 k ; y tÀ1 k À Á is obtained by: This is consistent with the step length as the intermittent interval length [74]. After calculating the step length, the previous bending point is updated by: The power law distribution of the step length is estimated with respect to the frequency distribution of D. In our analysis in this paper, we define a max ¼ 2p=9.

Analysis of the translation, splash and tornado behaviors
As mentioned before, the critical behavior in artificial systems indicates universal and efficient computation. This behavior is characterized by the complex mixing of contraction and extension of information. In cellular automata, the critical behavior consists of locally stable oscillatory or fixed patterns and chaotic wave patterns propagating from one local site to another [22,23]. Natural and real animal groups also exhibit these behaviors. Locally fixed patterns and locally oscillating patterns are compared to translating movements (or schooling) and tornado (or massive tornado) patterns, respectively [65,67,68]. The chaotic propagating waves are compared to the splash patterns of animal groups. While the coexistence of translation, splash and tornado behaviors seems to be the attribute specific to the critical behavior, previous swarm models have never revealed the coexistence of these behaviors. Indeed, splash and tornado patterns require fine tuning of the parameters and initial conditions. Therefore, it is very important to estimate the coexistence of translation, splash and tornado patterns in swarms in terms of the critical behavior. Here, we define the index for the tornado pattern in a swarm by the number of agents satisfying the following: Inequality (21) implies that the average bending angle of the agent's walk is so large that the agent rotates around a point. Since the splash pattern implies that agents are radially dispersed, the index for the splash pattern is defined by the number of agents satisfying the following: where r t min denotes the distance between the kth agent and its nearest neighbor, and T represents a constant time interval. Thus, inequality (22) implies that the nearest neighbors come away. In contrast, since the translation pattern of a swarm implies that agents move without changing their moving direction, the index for the translation pattern is defined by the number of agents satisfying the following: Figure 2 (top) shows snapshots of the three patterns of swarm behaviors. The left snapshot depicts the tornado patterns in which all 500 agents rotate. The middle snapshot shows the splash patterns in which 500 agents are dispersed from the right to the left. The right snapshot shows the translation patterns from below to above. These snapshots are simulated by the SPP-based model implemented to reveal each pattern. The bottom graph shows the number of agents satisfying conditions (21), (22) and (23) over time for each swarm pattern. In the left graph, first, all agents match their velocity resulting in the translation pattern, after which all agents rotate. Thus, soon after starting to rotate, the number of agents satisfying condition (21) (i.e., the tornado index) increases. In the middle graph, the number of agents satisfying condition (22) (i.e., the splash index) gradually increases. In the right graph, the number of agents satisfying condition (23) (i.e., the translation index) increases, while there are some agents exhibiting the splash pattern due to the perturbation in the SPP.
Note that the three indexes are not complementary to each other, and one agent can therefore satisfy multiple conditions. Notwithstanding this ambiguity, Figure 2 shows that the three indexes are suitable factors to estimate the behavioral components of the splash, tornado and translation patterns.

Results
All the experiments were conducted on a PC with an Intel Core 17 processor running at 2.6GHz, and hard drive of 16 Gbytes. Our implementation was compiled using gcc (4.2.1).

Critical behavior consisting of splash, tornado and translation patterns
In this paper, we use the term criticality both in the phenomenological sense and in the strict sense. Criticality in the strict sense is defined here by the phenomena characterized by a power law distribution. In contrast, criticality in the phenomenological sense implies the midpoint of the transition between order and disorder in a broad sense. This phase transition is expressed as the phase transition between the ordered swarm showing translation and/or tornado behavior and the disordered swarm showing splash behavior. While tornado behavior is a typical exploitation behavior, splash and translation behaviors are typical exploration behaviors, in another sense. The coexistence of tornado, splash and translation behaviors implies criticality in the phenomenolog-ical sense. Therefore, the question arises whether criticality in the phenomenological sense entails criticality in the strict sense. Figure 3 shows the time development of the SPP model. Each diagram shows a snapshot of the swarm consisting of 1000 agents, where the kth agent is represented by the line connecting its location x tÀ1 k ; y tÀ1 k À Á with x t k ; y t k À Á . Although the SPP does not implement flock centering (the force to approach the denser swarm) but only velocity matching, the agents are gradually concentrated due to the periodic boundary condition (i.e., the left margin of the square space is connected to the right margin, and the top of it is con-  nected to its bottom). This results in a dense translation pattern, as shown in the right bottom snapshot. Typical time development is stored as Video_simple_SPP_1.
It is known that the SPP exhibits a phase transition expressed as polarization (degree of velocity matching) with respect to normalized perturbations. The original SPP, which implements only velocity matching, could reveal neither tornado nor splash patterns. If the SPP is coupled with both a self-propelled force and friction, a specific ratio of the force and friction could lead to selforganizing tornado patterns [67,68]. Not only the SPP but also other swarm models could entail tornado patterns under attraction and repulsion balancing [39,[65][66]. Tornado patterns require fine tuning of the parameter setting. Although the splash pattern is frequently observed in real bird flocks and fish schools, it has been examined with respect to the prey-predator scheme in simulation studies. Since the prey-predator scheme is implemented by balancing attraction and repulsion, not only the prey-predator scheme but also the general swarm model coupled with attraction and repulsion could reveal splash patterns. However, fine tuning of the parameters is required.
Our SPP coupled with Bayesian inference contains both decisions consistent with and contradictory to velocity matching. Thus, it could be similar to the balancing of the self-propelled force and friction and the balancing of attraction and repulsion. This suggests that the SPP coupled with only Bayesian inference (i.e., BO model) could reveal tornado and splash patterns. Figure 4 shows short trajectories of the SPP coupled only with Bayesian inference. Since the degree of obeying the velocity matching is determined by Bayesian inference, this behavior depends on the likelihood of the hypotheses, where P k ðdjhÞ is invariant throughout time for each kth agent. In our model, for any k, P k djh ð Þ ¼ 0:7 if d ¼ h; otherwise, P k djh ð Þ ¼ 0:1. Through time development, P k djh ð Þ is not altered, and initially, P k h ð Þ ¼ 0:25 for any h 2 H.
While each agent affects P k h ð Þ dependent on its previous data, d t k , no behaviors characterized by critical behavior are observed. Each agent frequently changes its decision to either obey the velocity matching or not, and various rotations against the mean velocity are then mutually canceled. This leads to a fluctuating moving swarm, which sometimes indicates a large swarm and is sometimes divided into small parts. Dispersing and gathering patterns are perpetually iterated and indicate a complex fluctuating behavior, although they never show explicit splash and/or tornado patterns. If P k ðdjhÞ is randomly given and remains invariant over time, while the swarm rarely shows a tornado. However, in these cases, the realized tornado is perpetually sustained and not broken. This case seems to be achieved by the random setting ofP k ðdjhÞ and could be compared to the fine parameter tuning by chance. Typical time development of SPP only with Bayesian inference is stored as Video_BO_1. Figure 5 shows a pair of snapshots of the SPP coupled with Bayesian and inverse Bayesian inference (i.e., BIB model). Typical time development of SPP with Bayesian and inverse Bayesian inference is stored as Video_BIB_1. It is evident that there are translation, splash and tornado patterns. Due to the coexistence of translation, splash and tornado behaviors, the swarm is perpetually dispersed and gathered. Since the initial likelihood ofP t k ðdjhÞ is rapidly replaced by another one and is perpetually altered, complex patterns consisting of translation, splash and tornado behaviors occur independent of the initial condition of the distribution of agents and the initial condition of P t k ðdjhÞ. Compared to the SPP with only Bayesian inference (i.e., BO model), there are distinct patterns of locally stable information processing manifested as a tornado pattern and of information transmission manifested as splash and translation patterns. In the BO, information processing and transmission are mixed and averaged, which involves perturbed and fuzzy information processing. Since the BO seems neither chaotic nor exhibits definite information processing, it cannot be used for universal and efficient information processing.
The indexes of the tornado, splash and translation patterns are estimated here for the behavior of the swarm. Figure 6 shows a comparison of the indexes in the SPP and the BIB. The indexes of the tornado, splash and translation patterns are plotted over time. In the SPP, there is no indication of the splash pattern, and the swarm shows a mixture of translation and tornado patterns due to the perturbation, while the major phenomenon is translation. After a very large swarm is generated, translation occurs due to velocity matching. Thereafter, due to the perturbation, the swarm is divided into various parts, and one swarm is again formed. In the dividing process, a small moving population is split to maintain a dense population. Therefore, swarm deformation does not influence the index of the splash pattern but affects that of the tornado pattern.
In contrast, the BIB is characterized by the coexistence of tornado, splash and translation behaviors. As mentioned above, since the indexes of the tornado, splash and translation patterns are not independent of each other, they are mixed, and the perturbed splash and/or perturbed tornado pattern can thus be estimated by both the tornado and splash indexes. As shown in Figure 6 (right), all indexes are high in the swarm of the BIB. Figure 7 shows a comparison of the BO and the BIB with respect to the tornado, splash and translation indexes. Both simulation results are obtained under low-fluctuation conditions. Although both simulations are performed for 1000 agents, the summation of the agents satisfying the three indexes is much smaller than 1000 for the BO. This implies that there are many agents satisfying neither tornado nor translation indexes whose average turn angles are over p/10. Since the BO allows agents to neglect velocity matching, agents can turn at a sharp angle. If this fluctuating behavior resulting from Bayesian inference is simply mixed with the agents obeying velocity matching, it is considered that there are many agents that can turn at an angle other than p/10. In contrast, the BIB reveals that both the tornado and translation indexes are satisfied by a high proportion of agents and indicates that the agents satisfying the three indexes are complementary with each other. This implies that splash, translation and tornado behaviors coexist in the swarm and that these three behaviors continuously connect with each other as if the information originating from the tornado pattern is effectively transmitted to other places by the splash and translation patterns.
The simulation results in Figure 7 are obtained under the condition of a small perturbation. The question arises whether the difference between the two kinds of SPPs, BO and BIB, results from this small perturbation. Therefore, we simulated under the condition of a large perturbation. Figure 8 shows the simulation results under a perturbation e of 0.2. The trend depicted in Figure 7 is also observed in Figure 8. This reveals that adequate combinations of local information processing and information transmission are realized by the BIB, independent of the perturbation magnitude.
We conducted a statistical test to assess the difference between the BIB and BO models with respect to the tornado-translation index and the splash-translation index (Figure 9). The normality was tested by the Shapiro-Wilk normality test, and normality was rejected (p < 0.001, all are smaller than 7.443e-14). The difference between the mean value of the ratio in the BIB model and that in the BO model was checked by the Wilcoxon rank sum test, and the results showed a significant difference between them under both conditions, small perturbation and large perturbation (p < 0.001, all are smaller than 2.2e-16). The statistical data for the test are shown in Table 1.
The tendencies found in Figure 9 are general trends in the BIB and BO models. While the behaviors in the BO model are basically perturbed translation, which is revealed by the tornado index, the behaviors in the BIB model are basically characterized by the coexistence of tornado, splash and translation behaviors. Thus, the tornado index normalized by the translation index and the splash index normalized by the translation index in the BIB model are smaller than those in the BO model. In addition, we simulated the special case of the SPP model, in which every individual has a different radius. In that case, individuals converge to a large swarm controlled by the smallest radius. Thus, it can be concluded that individuals with heterogeneous radii never generate the coexistence of translation, tornado and splash behaviors.

Power law distribution in the SPP with Bayesian and inverse Bayesian inference
Next, we determine whether the SPP with Bayesian and inverse Bayesian inference (BIB model) shows critical behavior in a term of power law distribution. In estimating the translation, splash and tornado behaviors, the periodic boundary condition is defined as to not disperse the swarm. In estimating the power law distribution of the step length, 1000 agents are initially located in a small   central area in an open space, and the walking pattern of any freely moving agents is estimated. Figure 10 shows the trajectories of the simple SPP (above) and the BIB swarm. Each column shows the conditions of e=0.001, 0.1, 0.2 and 0.4. The SPP model with e=0.001 (above left) shows ballistic trajectories in finite time since little perturbation influences the trajectories. Since the simulations were run for T=10000, most agents remained in a central area, and alignm ent interactions also contributed to the step length distribution under a large perturbation. The results of the power law distributions are the same as those under periodic boundary conditions. The more perturbations there are, the more rotations occur indicating plant-root-like networks. In contrast to the SPP, the BIB swarm reveals similar trajectory patterns independent of the extent of the perturbation. While there are no trajectories of the BO swarm in Figure 10, the apparent patterns are similar to those generated by the BIB, where there are no ballistic trajectories. First, we conducted a statistical test, the Kruskal-Wallis rank sum test, to determine whether there was a significant difference among the distributions of the step length in the SPP, BO and BIB models. Since the distribution of the step length in the BIB swarm is far from a normal distribution, the test is nonparametric, and we adopted the Kruskal-Wallis rank sum test. It was found that all differences among the distributions of the step length in the SPP, BO and BIB models were significant (p < 0.001, all are smaller than 2.2e-16). Figure 11 shows the comparison of the mean and variance among the BIB, BO and simple SPP models under various levels of noise. Noise levels 1, 2 and 4 imply that e is 0.1, 0.2 and 0.4, respectively.
For the trajectories from the central area in an open space, the distribution of the step length is estimated. As mentioned before, the walk bending at a smaller angle than 2p/9 is regarded as a straight walk. If the step length and its frequency are represented by D and f(D), respectively, the power law distribution is expressed as: and if 1 l 3, it is called the Lévy walk. Recent studies have demonstrated that animal walks are approximated as truncated power law distributions [46,47] such as: Regarding the step length distributions of the simple SPP, BO and BIB, we test whether the distribution is approximated by the truncated power law distribution or the exponential distribution with the Akaike Information Criterion (AIC) [75][76][77][78]. Truncated power law distribution is frequently used for data fitting in Lévy walk analysis, while there is a trade-off between efficiency of fitting and completeness in data fitting [79][80][81].
Two probability density functions are given, one for the truncated power law distribution, f x ð Þ ¼ ðl À 1Þ=ðx 1Àl min À x 1Àl max Þx Àl , and the other for the exponential distribution, f x ð Þ ¼ kexpðÀk x À x min ð Þ Þ , where x min is determined by using Kolmogorov-Smirnov statistics and x max is defined as the maximum value of the data [82]. After best-fit exponents for the truncated power law (l) and the exponential distribution (k) and loglikelihood are calculated, the AIC and Akaike weight for both models are calculated. Finally, the better-fitting model is determined based on the Akaike weight. Figure 12 shows the cumulative frequency distribution for the step length. We calculated the Akaike weight for the truncated power law distribution, 0.0 wpl 1.0, where the larger wpl is, the higher the likelihood of the truncated power law distribution is. In contrast, 1-wpl represents the likelihood for the exponential distribution. The left, middle and right columns represent the conditions with e equal to 0.1, 0.2, 0.4, respectively. In the case of the simple SPP (above panel in Figure 12), for e = 0.1, l=1.00, k=0.00, wpl=0.00, for e = 0.2, l=1.00, k=0.00, wpl=0.00, and for e = 0.4, l=3.00, k=0.03, wpl=0.00, respectively. This implies that the frequency distribution of the SPP best fits to an exponential distribution. In the case of the BO, for e = 0.1, l=3.00, k=0.09, wpl=0.00, for e = 0.2, l=3.00, k=0.08, wpl=0.00, and for e = 0.4, l=3.00, k=0.10, wpl=0.00. This also implies that the frequency distribution of the BO model best fits to an exponential distribution. In contrast, the BIB reveals different types of distributions. In the case of the BIB, for e = 0.1, l=2.65, k=0.02, wpl=1.00, for e = 0.2, l=2.59, k=0.02, wpl=1.00, and for e = 0.4, l=2.72, k=0.02, wpl=1.00. Thus, the step length distribution generated by the BIB model strictly indicates a power law distribution, especially Lévy walk, and critical behavior.  We checked whether the window size of the BIB inference (m) influenced the results of the step length distribution. We found that the window size did not influence the power law distribution, as shown in Figure 13. The exponent for the truncated (l) and exponential (k) power law for various windows were obtained by the following: for m=40, l =2.81 and k = 0.011; for m=20, l =2.85 and k = 0.010; for m=10, l =2.78 and k = 0.011; and for m=5, l =2.8 and 5k = 0.010. The Akaike weights, wpl, for all window sizes were 1.0, and then all step length distributions could be fit to a truncated power law distribution.
Finally, it is confirmed that the complex patterns consisting of tornado, splash and translation behaviors, generated by the BIB, imply critical behaviors. In other words, the interaction between data and hypotheses in Bayesian inference can self-organize its critical behavior. This kind of criticality is intrinsically different from the criticality of the phase transition. As mentioned before, a phase transition occurs in the SPP. With respect to the step length distribution, one cannot determine a power law distribution, and there are signs of critical phenomena only encountered at the critical point. In contrast, the criticality in the BIB is not observed in narrow critical regions. Over a wide range of perturbations (e) and window size (m), power law distribution can be found. In this sense, the criticality is ubiquitously found and could be called a universal criticality.

Discussion and Conclusion
First, we consider the problem of artificial and natural design to be realized in open environments. Since the process of genome editing was first proposed and developed, it has become possible to design an artificial genome specific to a particularly concrete function [83][84][85]. Artificial design involves the hard problem of which optimal design based on the above concrete function must be achieved in an open and real environment. There are two reasons why the artificial and natural design is challenging. The first is how to overcome the trade-off between the behavior under open conditions and the behavior specific to a certain function [1,2]. The second reason is, if one could assign a solution beyond the above trade-off, how could the solution be effectively obtained [3,4]. These two reasons are mutually related to each other. Herein, we proposed a kind of solution to this problem.
The above design faces the trade-off between open conditions and definite functions, and the trade-off can be replaced by the trade-off between universality and efficiency in computation and by the trade-off between exploration and exploitation. In this sense, the design to be achieved results in the coexistence of both exploration and exploitation beyond the trade-off. In a swarm model, called the SPP, the trade-off between exploration and exploitation is expressed as the phase transition between order Figure 13. Cumulative frequency distributions of the step length for various window sizes of the BIB inference. The data are approximated by two distributions, i.e., the exponential (blue) and truncated power law (green) distributions. The perturbation e is set to 0.001. and chaos. Thus, it is expected that universal and efficient design can be realized at the critical point or the edge of chaos in the SPP parameter (perturbation) space.
Although determining the critical point as the design solution requires parameter fine tuning, there have been attempts to autonomously identify the critical point, that is called metaheuristics. One of these attempts is the self-organizing criticality (SOC). While the SOC has been applied to determining the criticality in the SPP, previous attempts have failed to identify the critical behavior without ad hoc global knowledge such as a fitness function. The metaheuristics for the swarm model based on SOC require global knowledge such as a fitness function [63,64]. As such, the idea of the SOC is not autonomously easy beyond parameter tuning. While there are many swarm-based optimization techniques [33][34][35][36], they are not related to the critical state.
Our proposal based on the interaction between data and hypotheses in Bayesian inference easily achieves critical behavior. However, this is different from the idea of the SOC with respect to two points. First, while the SOC is a way to choose the optimal solution from among possible states, our proposal does not select a limited state from many possible states. Instead of choosing, our proposal turns most states into critical states. In particular, metaheuristics are not introduced as fitness is added to the phase transition but instead the phase transition itself is modified. Therefore, critical behavior is not achieved only at the edge of chaos but is achieved anywhere in the parameter space. Second, our proposal never requires fitness knowledge or metaheuristics. The agent of the SPP never sees all agents in a space, and each agent makes decisions based on Bayesian and inverse Bayesian (BIB) inference, which is a task based not on global knowledge but on local dynamic knowledge.
Bayesian inference implements information contraction since it replaces the probability with the conditional probability. The agent using Bayesian inference reduces the world to that experienced by itself. Thus, this occurs in the simple optimization framework, and Bayesian inference contributes to rapidly reaching the optimal solution. Inverse Bayesian inference implements information extension since it replaces the conditional probability with the probability. In other words, the agent perpetually perceives the real world (data) outside its own cognitive world (hypotheses). Because the combination of information contraction and extension is independent of the knowledge of fitness, it could be applied to various problems.
Although our proposal is based on agents implementing an inference system, it never assumes that the genome has the ability to make decisions even if our proposal is applied to, for instance, genome editing. The optimal design of a network or traveling salesman problem can be resolved by ant agents making decisions with the probability or pheromones. It is never implied that the network itself has the ability to make decisions. Bayesian and inverse Bayesian inference denote that the probability is not globally given but is temporally and locally defined. Balancing information contraction and extension in Bayesian and inverse Bayesian inference leads to the universal criticality.