Possibility of Controlling Self-Organized Patterns with Totalistic Cellular Automata Consisting of Both Rules like Game of Life and Rules Producing Turing Patterns

The basic rules of self-organization using a totalistic cellular automaton (CA) were investigated, for which the cell state was determined by summing the states of neighboring cells, like in Conway’s Game of Life. This study used a short-range and long-range summation of the cell states around the focal cell. These resemble reaction-diffusion (RD) equations, in which self-organizing behavior emerges from interactions between an activating factor and an inhibiting factor. In addition, Game-of-Life-type rules, in which a cell cannot survive when adjoined by too many or too few living cells, were applied. Our model was able to mimic patterns characteristic of biological cells, including movement, growth, and reproduction. This result suggests the possibility of controlling self-organized patterns. Our model can also be applied to the control of engineering systems, such as multirobot swarms and self-assembling microrobots.


Introduction
Self-organization phenomena, in which global structures are produced from purely local interactions, are found in fields ranging from biology to human societies. If these emergent processes were to be controlled, various applications would be possible in engineering fields. For this purpose, the conditions under which they emerge must be elucidated. This is one of the goals of the study of complex systems. Such control methods may allow the automatic construction of machines. They could also be applied to the control of engineering systems, such as controlling robot swarms or self-assembling microrobots. In addition, the methods allow for the growth of artificial organs or the support of ecosystem conservation, as well as clarifying the emergence of the first life on ancient earth.
Mathematical modeling of self-organization phenomena has two main branches: the mathematical analysis of reaction-diffusion (RD) equations, and discrete modeling using cellular automata (CA). The Turing pattern model is one type of RD model. This was introduced by Turing in 1952 [1], where he treated morphogenesis as the interaction between activating and inhibiting factors. Typically, this model achieves self-organization through the different diffusion coefficients for two morphogens, equivalent to an activating and an inhibiting factor. The general RD equations can be written as follows:  CA models reduce the computational load by removing the need to solve differential equations numerically. However, it is necessary to translate the partial differential equations of the RD equations into the transition rules of the CA model. In a study of up to present, no general method currently exists for doing this, with the exception of the ultra-discrete systems [24] used in certain soliton equations. Normally, the transition rules driving the emergent phenomena must be found by trial and error. Studies that attempted to derive the transition rules using genetic algorithms [11,25] demonstrated the difficulty of deriving generalized rules. The Ishida cell division model [12] is a hybrid of the Wolfram-and Conway-type CAs, in which part of the rules are totalistic transition rules, and the others are head-to-head-type transition rules, and these rules were discovered by trial and error. The Young model [13] is one of the 2D totalistic models that bridges the RD equations and CA model; this model is used to produce Turing patterns. Some other examples to produce Turing patterns are below. Adamatzky [26] studied a binary-cell-state eight-cell neighborhood (2) a cell with four or more living neighbors dies; (3) a cell with three living neighbors is born; and (4) a cell with two or three living neighbors survives. The typical patterns that emerge can be categorized into three groups: still, moving, or oscillation.
Micromachines 2018, 9, x FOR PEER REVIEW 3 of 19 give rise to the diversity of patterns-edge of chaos. However, this could not derive the rule to construct favorite patterns.  CA models reduce the computational load by removing the need to solve differential equations numerically. However, it is necessary to translate the partial differential equations of the RD equations into the transition rules of the CA model. In a study of up to present, no general method currently exists for doing this, with the exception of the ultra-discrete systems [24] used in certain soliton equations. Normally, the transition rules driving the emergent phenomena must be found by trial and error. Studies that attempted to derive the transition rules using genetic algorithms [11,25] demonstrated the difficulty of deriving generalized rules. The Ishida cell division model [12] is a hybrid of the Wolfram-and Conway-type CAs, in which part of the rules are totalistic transition rules, and the others are head-to-head-type transition rules, and these rules were discovered by trial and error. The Young model [13] is one of the 2D totalistic models that bridges the RD equations and CA model; this model is used to produce Turing patterns. Some other examples to produce Turing patterns are below. Adamatzky [26] studied a binary-cell-state eight-cell neighborhood CA models reduce the computational load by removing the need to solve differential equations numerically. However, it is necessary to translate the partial differential equations of the RD equations into the transition rules of the CA model. In a study of up to present, no general method currently exists for doing this, with the exception of the ultra-discrete systems [24] used in certain soliton equations. Normally, the transition rules driving the emergent phenomena must be found by trial and error. Studies that attempted to derive the transition rules using genetic algorithms [11,25] demonstrated the difficulty of deriving generalized rules. The Ishida cell division model [12] is a hybrid of the Wolframand Conway-type CAs, in which part of the rules are totalistic transition rules, and the others are head-to-head-type transition rules, and these rules were discovered by trial and error. The Young model [13] is one of the 2D totalistic models that bridges the RD equations and CA model; this model is used to produce Turing patterns. Some other examples to produce Turing patterns are below. Adamatzky [26] studied a binary-cell-state eight-cell neighborhood two-dimensional cellular automaton model with semitotalistic transitions rules. Dormann [27] also used a 2D outer-totalistic model with threes states to produce a Turing-like pattern. Tsai [28] analyzed a self-replicating mechanism in Turing patterns with a minimal autocatalytic monomer-dimer system.
Young's CA model uses a real number function to derive the distance effects, with two constant values within a grid cell: u 1 (positive) and u 2 (negative). The calculation begins by randomly distributing black cells on a rectangular grid. Then, for each cell at position R 0 in 2D fields, the next cell state of R 0 , due to all nearby black cells at position R i , are added up. R i is assumed to be a black cell within radius r 2 from R 0 cell. The function v(r) is a continuous step function, as shown in Figure 3b. The activation area, indicated by u 1 , has a radius of r 1 and the inhibition area, indicated by u 2 , has a radius of r 2 (r 2 > r 1 ). At position R 0 , when R i is assumed to be a grid within r 2 , function v(|R 0 − R i |) value is added up. Function |R 0 − R i | indicates the distance between R 0 and R i . If ∑ i v (|R 0 − R i |) > 0, the grid cell at point R 0 is marked as a black cell. If ∑ i v (|R 0 − R i |) < 0, the grid cell becomes a white cell. If ∑ i v (|R 0 − R i |) = 0, the grid cell does not change state [13]. Using these functions, Young described that a Turing pattern can be generated. Spot patterns or striped patterns can be created with relative changes between u 1 and u 2 . two-dimensional cellular automaton model with semitotalistic transitions rules. Dormann [27] also used a 2D outer-totalistic model with threes states to produce a Turing-like pattern. Tsai [28] analyzed a self-replicating mechanism in Turing patterns with a minimal autocatalytic monomer-dimer system. Young's CA model uses a real number function to derive the distance effects, with two constant values within a grid cell: u1 (positive) and u2 (negative). The calculation begins by randomly distributing black cells on a rectangular grid. Then, for each cell at position R0 in 2D fields, the next cell state of R0, due to all nearby black cells at position Ri, are added up. Ri is assumed to be a black cell within radius r2 from R0 cell. The function v(r) is a continuous step function, as shown in Figure 3b. The activation area, indicated by u1, has a radius of r1 and the inhibition area, indicated by u2, has a radius of r2 (r2 > r1). At position R0, when Ri is assumed to be a grid within r2, function v(|R0 − Ri|) value is added up. Function |R0 − Ri| indicates the distance between R0 and Ri. If ∑i v (|R0 − Ri|) > 0, the grid cell at point R0 is marked as a black cell. If ∑i v (|R0 − Ri|) < 0, the grid cell becomes a white cell. If ∑i v (|R0 − Ri|) = 0, the grid cell does not change state [13]. Using these functions, Young described that a Turing pattern can be generated. Spot patterns or striped patterns can be created with relative changes between u1 and u2. In this Young model, let u1 = 1, u2 = w (here 0 < w < 1), and further, if the state of the cell is set to 0 (white) and 1 (black), this model can be arranged as follows. The state of cell i is expressed as ( ) ( ( ) = [0, 1]) at time t. The following state ( + 1) at time + 1 is determined by the states of the neighboring cells. Here, N1 is the sum of the states of the domain within S meshes of the focal cell. Similarly, N2 is the sum of the states of the domain within t meshes of the focal cell, assuming that S < t.
Here, S (T) is the number of cells within s (t) meshes from focal cell. Figure 4 shows the schematic of the total sum of states N1 and N2. The next time state of the focal cell is determined by the following Expression (1): In this Young model, let u 1 = 1, u 2 = w (here 0 < w < 1), and further, if the state of the cell is set to 0 (white) and 1 (black), this model can be arranged as follows. The state of cell i is expressed as c i (t) (c i (t) = [0, 1]) at time t. The following state c i (t + 1) at time t + 1 is determined by the states of the neighboring cells. Here, N 1 is the sum of the states of the domain within S meshes of the focal cell. Similarly, N 2 is the sum of the states of the domain within t meshes of the focal cell, assuming that S < t.
Here, S (T) is the number of cells within s (t) meshes from focal cell. Figure 4 shows the schematic of the total sum of states N 1 and N 2 . The next time state of the focal cell is determined by the following Expression (1): Micromachines 2018, 9, 339 5 of 18 Cell state at the next time step= (1) Figure 4. Schematic of the summing of states N1 and N2. Each grid cell has state 0 (white) or 1 (black).
The inner area has a domain within s grids of the focal cell and the outer area a domain within t grids.
The essence of the Turing pattern model is that the short-and long-range spatial scales are each affected by separate factors [29], and pattern formation emerges from nonlinear interactions between the two factors. Turing used two chemicals with different diffusion coefficients to create these short-and long-range spatial effects. However, as long as there exists a difference between long-and short-range effects, other implementation methods can be applied. This model used two ranges of s mesh and t mesh to create a difference. It is therefore thought that this model is essentially the same as an RD model. Figure 5 shows the results for the Expression (1) model using a square grid. As the initial conditions, 300 state 1 cells were set randomly in the calculation field 80 × 80. Turing-like patterns emerged as parameter w changed. When w was greater than 0.4, spot patterns were observed, whereas at w = 0.3, stripe patterns emerged. When w was 0.20 or lower, all cells in the field had a value of 1 (shown in black). Changing the parameters N1 and N2 merely changed the scale of the patterns.
On the other hand, the GoL is one of the 2D totalistic CAs to emerge patterns of diversity. In this model, the state of the focal cell is activated when the total number of surrounding states is within a certain range. Therefore, we considered the extended model of GoL. Whereas the GoL uses the states of the eight surrounding cells (the Moore neighborhood), an extended neighborhood is possible that considers two or more meshes from the focal cell like Young model. Such models are expected to generate a similar variety of patterns to the GoL.
In natural phenomena, short-distance influences play a stronger role than more distant influences. Under this circumstance (1) The influence decreases as the distance from the focal cell increases, which is N1 + N2 × w type model (w: weight parameter, 0 < w < 1). However, reaction diffusion phenomena are created with long-distance influences that play a negative role than more near influences. (2) The influence reverses as the distance from the focal cell increases, which is N1 − N2 × w type model (w: weight parameter, 0 < w < 1). The essence of the Turing pattern model is that the short-and long-range spatial scales are each affected by separate factors [29], and pattern formation emerges from nonlinear interactions between the two factors. Turing used two chemicals with different diffusion coefficients to create these short-and long-range spatial effects. However, as long as there exists a difference between long-and short-range effects, other implementation methods can be applied. This model used two ranges of s mesh and t mesh to create a difference. It is therefore thought that this model is essentially the same as an RD model. Figure 5 shows the results for the Expression (1) model using a square grid. As the initial conditions, 300 state 1 cells were set randomly in the calculation field 80 × 80. Turing-like patterns emerged as parameter w changed. When w was greater than 0.4, spot patterns were observed, whereas at w = 0.3, stripe patterns emerged. When w was 0.20 or lower, all cells in the field had a value of 1 (shown in black). Changing the parameters N 1 and N 2 merely changed the scale of the patterns.
On the other hand, the GoL is one of the 2D totalistic CAs to emerge patterns of diversity. In this model, the state of the focal cell is activated when the total number of surrounding states is within a certain range. Therefore, we considered the extended model of GoL. Whereas the GoL uses the states of the eight surrounding cells (the Moore neighborhood), an extended neighborhood is possible that considers two or more meshes from the focal cell like Young model. Such models are expected to generate a similar variety of patterns to the GoL.
In natural phenomena, short-distance influences play a stronger role than more distant influences. Under this circumstance (1) The influence decreases as the distance from the focal cell increases, which is N 1 + N 2 × w type model (w: weight parameter, 0 < w < 1). However, reaction diffusion phenomena are created with long-distance influences that play a negative role than more near influences. (2) The influence reverses as the distance from the focal cell increases, which is N 1 − N 2 × w type model (w: weight parameter, 0 < w < 1). Micromachines 2018, 9, x FOR PEER REVIEW 6 of 19 Figure 5. Results for the model of Expression (1) with s = 5 and t = 9 in a square grid. As the initial condition, 300 state 1 cells were set randomly in the calculation field. Turing-like patterns emerged as parameter w changed. When w was greater than 0.4, spot patterns were observed, whereas at w = 0.3, stripe patterns emerged. When w was 0.20 or lower, all cells in the field had a value of 1 (shown in black).
The first type of the model is expected to produce results similar to those of other totalistic CAs. Under specific conditions, GoL-like patterns will be found. The second model is expected to produce results similar to those of a Young model, which similarly takes account of the nearness of the activating and inhibiting factors.
This study adopted the second type of the model, and the emergence of a variety of patterns, such as the Turing pattern, were confirmed. In addition, by applying GoL-type rules to the second model type, more complex patterns are expected, such as self-reproducing and growth patterns. Our study model was able to produce not only a Turing-like pattern while using a simple transition expression, but also the model produces various patterns such as still, moving, and oscillation. Furthermore, this model can control these patterns with two parameters, which means the possibility of controlling self-organized patterns with a totalistic CA model.

Model
We investigated the results of an N1 − N2 × w type totalistic CA model in a 2D field that was implemented in Java. The GoL-like rules assumed in this model which limit the range of survival as state 1. To address this, by extending Expression (1), we introduced Expression (2). This expression was obtained by adding Rule 4 to Expression (1). When N2 = 0, it is assumed that Rule 4 takes precedence.
Cell state of the next time step Figure 5. Results for the model of Expression (1) with s = 5 and t = 9 in a square grid. As the initial condition, 300 state 1 cells were set randomly in the calculation field. Turing-like patterns emerged as parameter w changed. When w was greater than 0.4, spot patterns were observed, whereas at w = 0.3, stripe patterns emerged. When w was 0.20 or lower, all cells in the field had a value of 1 (shown in black).
The first type of the model is expected to produce results similar to those of other totalistic CAs. Under specific conditions, GoL-like patterns will be found. The second model is expected to produce results similar to those of a Young model, which similarly takes account of the nearness of the activating and inhibiting factors.
This study adopted the second type of the model, and the emergence of a variety of patterns, such as the Turing pattern, were confirmed. In addition, by applying GoL-type rules to the second model type, more complex patterns are expected, such as self-reproducing and growth patterns. Our study model was able to produce not only a Turing-like pattern while using a simple transition expression, but also the model produces various patterns such as still, moving, and oscillation. Furthermore, this model can control these patterns with two parameters, which means the possibility of controlling self-organized patterns with a totalistic CA model.

Model
We investigated the results of an N 1 − N 2 × w type totalistic CA model in a 2D field that was implemented in Java. The GoL-like rules assumed in this model which limit the range of survival as state 1. To address this, by extending Expression (1), we introduced Expression (2). This expression was obtained by adding Rule 4 to Expression (1). When N 2 = 0, it is assumed that Rule 4 takes precedence.
Cell state of the next time step= Here, 1 > w 2 > w 1 > 0, and N 1 is the sum of the states of the domain within s meshes of the focal cell. Similarly, N 2 is the sum of the states of the domain within t meshes of the focal cell, assuming that s < t. S (T) is the number of cells within s (t) meshes from focal cell. As with the Young model, we used "Unchange" as the equality in Expression (2) to prevent all cells becoming state 1 or 0 in the next time step when N 1 = 0 and N 2 = 0.
The model used 2D hexagonal grids ( Figure 6), in which the transition rules were simple to apply. Although square grids are generally used in 2D CA modeling, we also used hexagonal grids for the reason that a hexagonal grid is isotropic as opposed to a square grid. The models were implemented under the following conditions: -Calculation field: 80 × 80 cells in hexagonal grids -Periodic boundary condition -Initial conditions were assumed following two types of conditions: Condition (a): some state 1s were set randomly in the calculation field; Condition (b): some state 1s were set in the center of the calculation fields ( Figure 7).
Here, 1 > w2 > w1 > 0, and N1 is the sum of the states of the domain within s meshes of the focal cell. Similarly, N2 is the sum of the states of the domain within t meshes of the focal cell, assuming that s < t. S (T) is the number of cells within s (t) meshes from focal cell. As with the Young model, we used "Unchange" as the equality in Expression (2) to prevent all cells becoming state 1 or 0 in the next time step when N1 = 0 and N2 = 0.
The model used 2D hexagonal grids ( Figure 6), in which the transition rules were simple to apply. Although square grids are generally used in 2D CA modeling, we also used hexagonal grids for the reason that a hexagonal grid is isotropic as opposed to a square grid. The models were implemented under the following conditions: -Calculation field: 80 × 80 cells in hexagonal grids -Periodic boundary condition -Initial conditions were assumed following two types of conditions: Condition (a): some state 1s were set randomly in the calculation field; Condition (b): some state 1s were set in the center of the calculation fields ( Figure 7).   Here, 1 > w2 > w1 > 0, and N1 is the sum of the states of the domain within s meshes of the focal cell. Similarly, N2 is the sum of the states of the domain within t meshes of the focal cell, assuming that s < t. S (T) is the number of cells within s (t) meshes from focal cell. As with the Young model, we used "Unchange" as the equality in Expression (2) to prevent all cells becoming state 1 or 0 in the next time step when N1 = 0 and N2 = 0.
The model used 2D hexagonal grids (Figure 6), in which the transition rules were simple to apply. Although square grids are generally used in 2D CA modeling, we also used hexagonal grids for the reason that a hexagonal grid is isotropic as opposed to a square grid. The models were implemented under the following conditions: -Calculation field: 80 × 80 cells in hexagonal grids -Periodic boundary condition -Initial conditions were assumed following two types of conditions: Condition (a): some state 1s were set randomly in the calculation field; Condition (b): some state 1s were set in the center of the calculation fields ( Figure 7).

Results
Figure 8 maps the model results against parameters w 1 and w 2 . Initial condition is based on 300 randomly arranged state 1 cells. Different patterns were observed at different values of w 1 and w 2 . As in the Young-type model, the parameter w 1 determined the pattern type, whether spotted or striped. As w 2 became smaller, pulsing and unsteady white spots were born within the black pattern. The vibration of these white spots became more intense as the value of w 2 was reduced. Figure 9 shows the results when varying the number of states 1 (black cell) to place them randomly in the initial state. The calculation conditions were s = 5, t = 9, w 1 = 0.40, and w 2 = 0.80 in a hexagonal grid. Based on the results of this figure, if there are more than a certain number of states 1, it is evident that a similar pattern formation develops. Furthermore, Figure 10 shows the transition of the ratio of states 1 and 2, and the ratio of the applied cells of Rule 4 at each time step. This result was calculated when w 2 was changed (s = 5, t = 9, and w 1 = 0.35). When w 2 was 1, Rule 4 was not applied, and the result was equivalent to the Young model. As w 2 was lowered, the application ratio of Rule 4 increased. At this time, white regions were generated in the black spots, and the instability increased and became pulsating. This phenomenon is attributed to Rule 4.

Results
Figure 8 maps the model results against parameters w1 and w2. Initial condition is based on 300 randomly arranged state 1 cells. Different patterns were observed at different values of w1 and w2. As in the Young-type model, the parameter w1 determined the pattern type, whether spotted or striped. As w2 became smaller, pulsing and unsteady white spots were born within the black pattern. The vibration of these white spots became more intense as the value of w2 was reduced. Figure 9 shows the results when varying the number of states 1 (black cell) to place them randomly in the initial state. The calculation conditions were s = 5, t = 9, w1 = 0.40, and w2 = 0.80 in a hexagonal grid. Based on the results of this figure, if there are more than a certain number of states 1, it is evident that a similar pattern formation develops. Furthermore, Figure 10 shows the transition of the ratio of states 1 and 2, and the ratio of the applied cells of Rule 4 at each time step. This result was calculated when w2 was changed (s = 5, t = 9, and w1 = 0.35). When w2 was 1, Rule 4 was not applied, and the result was equivalent to the Young model. As w2 was lowered, the application ratio of Rule 4 increased. At this time, white regions were generated in the black spots, and the instability increased and became pulsating. This phenomenon is attributed to Rule 4. Figure 11 shows the effect of the initial condition which is centrally arranged on four black cells. The instability became more pronounced as w2 decreased, and the central spotted pattern grew as an unstable wave in the calculation field. By contrast, when w2 was relatively large, the central spot did not spread.   Figure 11 shows the effect of the initial condition which is centrally arranged on four black cells. The instability became more pronounced as w 2 decreased, and the central spotted pattern grew as an unstable wave in the calculation field. By contrast, when w 2 was relatively large, the central spot did not spread.  This result was calculated when w2 was changed in s = 5, t = 9, and w1 = 0.35. When w2 was 1, Rule 4 was not applied, and it was equivalent to the Young model. As w2 was lowered, the application ratio of Rule 4 increased. In comparison with Figure 8, white regions were generated in the black spots, and instability increased and began pulsating. This phenomenon is attributed to Rule 4.   This result was calculated when w2 was changed in s = 5, t = 9, and w1 = 0.35. When w2 was 1, Rule 4 was not applied, and it was equivalent to the Young model. As w2 was lowered, the application ratio of Rule 4 increased. In comparison with Figure 8, white regions were generated in the black spots, and instability increased and began pulsating. This phenomenon is attributed to Rule 4. This result was calculated when w 2 was changed in s = 5, t = 9, and w 1 = 0.35. When w 2 was 1, Rule 4 was not applied, and it was equivalent to the Young model. As w 2 was lowered, the application ratio of Rule 4 increased. In comparison with Figure 8, white regions were generated in the black spots, and instability increased and began pulsating. This phenomenon is attributed to Rule 4. Figure 11. Results for the model with s = 5 and t = 9 in a hexagonal grid. The image of each w1 and w2 value is based on centrally arranged four black cells as shown in Figure 7. Instability became stronger as w2 decreased. The initial central spotted pattern grew as an unstable wave in the field. When w2 was relatively large, the central spot did not spread.  Figure 12 can be viewed in the Video S1 of Supplementary Materials. As the initial condition, four state-1 cells (black) were set in the center of the field shown in Figure 7. The central spotted pattern was divided into two spots, and these divisions spread until the field was filled with spotted patterns. Figure 13 shows the results when the initial condition of state 1 was subjected to the same conditions as in Figure 12 (s = 4, t = 8, w1 = 0.40, and w2 = 0.75 in a hexagonal grid). When there was only one state 1 in the initial field, Rule 4 was applied around this state 1, which disappeared at the next step. If there are two or more states 1 as the initial state, it was observed that the state 1 survived, and pattern formation similar to Figure 12 occurred. In the case of several state 1s, due to the balance between Rules 3 and 4, complicated patterns were observed in two or three steps from the beginning of the calculation, and thereafter, a self-replicating pattern was observed.
In addition, Figure 14 shows the calculation results when increasing the black region as the initial condition, under the same conditions as in Figure 12 (s = 4, t = 8, w1 = 0.40, and w2 = 0.75 in a hexagonal grid). The point symmetrical initial shape tended to be a fixed pattern due to the balance of the rules. On the other hand, if the black area was larger than a certain size, Rule 1 was mostly applied, and disappeared at the second step. To produce a self-replicating pattern, an asymmetric shape must initially occur. Figure 15 shows the transition of the ratios of states 1 and 2, and the ratio of applied cells of Rule 4 per time step under the same conditions as in Figure 12 (s = 5, t = 9, w1 = 0.40, and w2 = 0.75). In 90 to 100 steps, the frequency of pattern division became high, and the application ratio of Rule 4 increased. Rule 4 generated a white area in the black spot patterns, and instability increased; hence, the symmetrical shape collapsed and became a trigger for division. Figures 16 and 17 show the results with an initial pattern of a symmetrical ring with state 1s. In the case where the initial shape was symmetrical, the still pattern or the oscillation pattern often appeared in a wide range of w1 and w2. When w1 = 0.45 and w2 = 0.55, the ring pattern spread while dividing spots. In an opposite way, as w2 was approached as w1, the instability increased and the pattern disappeared. Figure 11. Results for the model with s = 5 and t = 9 in a hexagonal grid. The image of each w 1 and w 2 value is based on centrally arranged four black cells as shown in Figure 7. Instability became stronger as w 2 decreased. The initial central spotted pattern grew as an unstable wave in the field. When w 2 was relatively large, the central spot did not spread. Figure 12 shows the time series results of the model with s = 4, t = 8, w 1 = 0.40, and w 2 = 0.75. In this figure, i shows the iteration number of the calculation. The video of Figure 12 can be viewed in the Video S1 of Supplementary Materials. As the initial condition, four state-1 cells (black) were set in the center of the field shown in Figure 7. The central spotted pattern was divided into two spots, and these divisions spread until the field was filled with spotted patterns. Figure 13 shows the results when the initial condition of state 1 was subjected to the same conditions as in Figure 12 (s = 4, t = 8, w 1 = 0.40, and w 2 = 0.75 in a hexagonal grid). When there was only one state 1 in the initial field, Rule 4 was applied around this state 1, which disappeared at the next step. If there are two or more states 1 as the initial state, it was observed that the state 1 survived, and pattern formation similar to Figure 12 occurred. In the case of several state 1s, due to the balance between Rules 3 and 4, complicated patterns were observed in two or three steps from the beginning of the calculation, and thereafter, a self-replicating pattern was observed.
In addition, Figure 14 shows the calculation results when increasing the black region as the initial condition, under the same conditions as in Figure 12 (s = 4, t = 8, w 1 = 0.40, and w 2 = 0.75 in a hexagonal grid). The point symmetrical initial shape tended to be a fixed pattern due to the balance of the rules. On the other hand, if the black area was larger than a certain size, Rule 1 was mostly applied, and disappeared at the second step. To produce a self-replicating pattern, an asymmetric shape must initially occur. Figure 15 shows the transition of the ratios of states 1 and 2, and the ratio of applied cells of Rule 4 per time step under the same conditions as in Figure 12 (s = 5, t = 9, w 1 = 0.40, and w 2 = 0.75). In 90 to 100 steps, the frequency of pattern division became high, and the application ratio of Rule 4 increased. Rule 4 generated a white area in the black spot patterns, and instability increased; hence, the symmetrical shape collapsed and became a trigger for division. Figure 12. Results for the model with s = 4, t = 8, w1 = 0.40, and w2 = 0.75 in a hexagonal grid. As the initial condition, four state 1 cells (black) were set in the center of the field as shown in Figure 7. The central spotted pattern was divided into two spots, which then spread until the field was filled with spotted patterns. Figure 13. Results with the initial condition of state 1 under the same conditions as in Figure 12 (s = 4, t = 8, w1 = 0.40, and w2 = 0.75 in a hexagonal grid). If there is only one state 1 in the initial field, Rule 4 was applied at that location, which disappeared at the next step. If there are two or more state 1s as the initial state, state 1 was observed to survive, and a pattern formation similar to Figure 12 occurred. Figure 12. Results for the model with s = 4, t = 8, w 1 = 0.40, and w 2 = 0.75 in a hexagonal grid. As the initial condition, four state 1 cells (black) were set in the center of the field as shown in Figure 7. The central spotted pattern was divided into two spots, which then spread until the field was filled with spotted patterns.
Micromachines 2018, 9, x FOR PEER REVIEW 11 of 19 Figure 12. Results for the model with s = 4, t = 8, w1 = 0.40, and w2 = 0.75 in a hexagonal grid. As the initial condition, four state 1 cells (black) were set in the center of the field as shown in Figure 7. The central spotted pattern was divided into two spots, which then spread until the field was filled with spotted patterns. If there is only one state 1 in the initial field, Rule 4 was applied at that location, which disappeared at the next step. If there are two or more state 1s as the initial state, state 1 was observed to survive, and a pattern formation similar to Figure 12 occurred. Figure 13. Results with the initial condition of state 1 under the same conditions as in Figure 12 (s = 4, t = 8, w 1 = 0.40, and w 2 = 0.75 in a hexagonal grid). If there is only one state 1 in the initial field, Rule 4 was applied at that location, which disappeared at the next step. If there are two or more state 1s as the initial state, state 1 was observed to survive, and a pattern formation similar to Figure 12 occurred. Figures 16 and 17 show the results with an initial pattern of a symmetrical ring with state 1 s. In the case where the initial shape was symmetrical, the still pattern or the oscillation pattern often appeared in a wide range of w 1 and w 2 . When w 1 = 0.45 and w 2 = 0.55, the ring pattern spread while dividing spots. In an opposite way, as w 2 was approached as w 1 , the instability increased and the pattern disappeared.  Figure 12 (s = 4, t = 8, w1 = 0.40, and w2 = 0.75 in a hexagonal grid). The point symmetrical initial shape tended to be a fixed pattern, due to the balance of each of the rules. On the other hand, if the black area is larger than a certain size, Rule 1 is mostly applied, and disappeared at the second step. To produce a self-replicating pattern, an asymmetric shape must initially occur.  Figure 12 (s = 5, t = 9, w1 = 0.40, and w2 = 0.75). From 90 to 100 steps, the frequency of pattern division became high, and the application ratio of Rule 4 increased. Figure 18 shows the initial pattern of a deformed ring with a thick upper boundary. When w1 = 0.45 and w2 = 0.65, the ring pattern was divided while moving. When w1 = 0.45 and w2 = 0.72, a constantly moving ring-like pattern formed. When w1 = 0.45 and w2 = 0.75, a single fixed symmetrical ring pattern formed. Thus, when the initial shape is asymmetric, a moving pattern or dividing pattern tends to be generated from the breaking of symmetry. As w2 approaches 1, it approaches fixed patterns, and results in a pattern close to the Turing pattern. Figure 14. Calculation results, when increasing the black region as the initial condition, under the same conditions as in Figure 12 (s = 4, t = 8, w 1 = 0.40, and w 2 = 0.75 in a hexagonal grid). The point symmetrical initial shape tended to be a fixed pattern, due to the balance of each of the rules. On the other hand, if the black area is larger than a certain size, Rule 1 is mostly applied, and disappeared at the second step. To produce a self-replicating pattern, an asymmetric shape must initially occur.
Micromachines 2018, 9, x FOR PEER REVIEW 12 of 19 Figure 14. Calculation results, when increasing the black region as the initial condition, under the same conditions as in Figure 12 (s = 4, t = 8, w1 = 0.40, and w2 = 0.75 in a hexagonal grid). The point symmetrical initial shape tended to be a fixed pattern, due to the balance of each of the rules. On the other hand, if the black area is larger than a certain size, Rule 1 is mostly applied, and disappeared at the second step. To produce a self-replicating pattern, an asymmetric shape must initially occur.  Figure 12 (s = 5, t = 9, w1 = 0.40, and w2 = 0.75). From 90 to 100 steps, the frequency of pattern division became high, and the application ratio of Rule 4 increased. Figure 18 shows the initial pattern of a deformed ring with a thick upper boundary. When w1 = 0.45 and w2 = 0.65, the ring pattern was divided while moving. When w1 = 0.45 and w2 = 0.72, a constantly moving ring-like pattern formed. When w1 = 0.45 and w2 = 0.75, a single fixed symmetrical ring pattern formed. Thus, when the initial shape is asymmetric, a moving pattern or dividing pattern tends to be generated from the breaking of symmetry. As w2 approaches 1, it approaches fixed patterns, and results in a pattern close to the Turing pattern.   . Results for the model with s = 5 and t = 8 in a hexagonal grid. The initial pattern was a symmetrical ring. When w1 = 0.45 and w2 = 0.65, the ring pattern became an oscillation pattern with period 2. When w1 = 0.45 and w2 = 0.72, the ring pattern likewise became an oscillation pattern with period 2. When w1= 0.45 and w2 = 0.75, a still pattern was formed.   Figure 18 shows the initial pattern of a deformed ring with a thick upper boundary. When w 1 = 0.45 and w 2 = 0.65, the ring pattern was divided while moving. When w 1 = 0.45 and w 2 = 0.72, a constantly moving ring-like pattern formed. When w 1 = 0.45 and w 2 = 0.75, a single fixed symmetrical ring pattern formed. Thus, when the initial shape is asymmetric, a moving pattern or dividing pattern tends to be generated from the breaking of symmetry. As w 2 approaches 1, it approaches fixed patterns, and results in a pattern close to the Turing pattern.

Discussion
In this study, we constructed a model in which the focal state 1 survives only when a moderate rate of the number of state 1 cells surrounds the focal cell. This model created a central white part (state 0) within the black domain by changing parameter w2. The patterns then became increasingly unstable over the calculation field, and this instability increased the rate of pattern emergence, including self-reproducing spot patterns and stripe patterns.
The self-reproducing spot pattern of Figure 12 seems to have some association with the Gray-Scott model, which is a type of RD model. Generally, the Gray-Scott model produces a Figure 18. Results for the model with s = 5 and t = 8 in a hexagonal grid. The initial pattern was a deformed ring with a thick upper boundary. When w 1 = 0.45 and w 2 = 0.65, the ring pattern was divided while moving. When w 1 = 0.45 and w 2 = 0.72, a moving ring-like pattern formed. When w 1 = 0.45 and w 2 = 0.75, a single fixed symmetrical ring pattern formed.

Discussion
In this study, we constructed a model in which the focal state 1 survives only when a moderate rate of the number of state 1 cells surrounds the focal cell. This model created a central white part (state 0) within the black domain by changing parameter w 2 . The patterns then became increasingly unstable over the calculation field, and this instability increased the rate of pattern emergence, including self-reproducing spot patterns and stripe patterns.
The self-reproducing spot pattern of Figure 12 seems to have some association with the Gray-Scott model, which is a type of RD model. Generally, the Gray-Scott model produces a self-reproducing spot pattern when the parameters are adjusted. In future work, we will investigate the relationship between the current model and the Gray-Scott model.
Particularly, in the case of the model with s = 5 and t = 8 shown in Figures 16-18, the result suggests the possibility of controlling self-organized patterns. When w 1 is increased, a spot shape can be generated, and by decreasing w 1 , a stripe pattern can be grown from spotted seeds.
In the symmetric shape, when w 2 is large, it is stable in many cases. However, by making the value of w 2 small and making it close to w 1 , the spotted pattern becomes unstable, and asymmetric shapes are born. Furthermore, these spots are moving or divided by the controlling w 2 parameter, as shown in Figure 19.
In a future study, if CA rules can translate to a chemical reaction process, a control technology for nanomachines becomes feasible. Thus, the simple totalistic CA presented in this paper allows the emergence of a wide range of self-organization patterns to be observed.
Based on these results, if the parameters w1 and w2 are changed at each time step, it is possible to generate a continuously changing pattern. Figure 20 shows a calculation case in which w1 and w2 are changed in time over time. The video of Figure 20 can be viewed in the Video S2 of Supplementary Materials. As in this result, it is possible to create changes by changing w1 and w2; birth of annulus -> growth -> division -> move -> decline -> growth -> annihilation. Figure 21 shows the transition of the ratio of states 1 and 2, and the ratio of the applied cells of Rule 4 in each time step during the calculation of Figure 20. Rule 4 was applied more frequently when spot patterns replicated, as compared to when the spot moved or vibrated. Furthermore, when the pattern disappeared, the application frequency of Rule 4 increased. Thus, it is suggested that pattern stability can be controlled by changing the application ratio of Rule 4.
Since this research is still in its first stage, discussions are limited mainly to qualitative considerations, but in future work, it will be necessary to investigate further quantitative assessments. There is a possibility of finding an index like Langton's λ parameter with the evaluation of many calculation cases. Furthermore, since it has rules similar to the Game of Life, basically like the Game of Life, many trials are necessary to find the desired pattern. Slight changes in the initial value may develop into different patterns. In the future, it is necessary to make a catalog of pattern formation on the relation between the initial condition and patterns based on many calculation results.   Figure 13. The ring pattern can be controlled with parameter w 2 . When w 2 = 0.65, the ring pattern was divided while moving. When w 2 = 0.72, a moving ring pattern formed. When w 2 = 0.75, a fixed ring pattern formed.
The robustness of the initial conditions was found to have the following tendencies. (1) If state 1 is one, it disappears; (2) When state 1 is randomly arranged, if two or more states 1 are distributed in the range for counting the number of surrounding states, the development of the pattern is recognized. Next, depending on the degree of instability due to the w 2 parameter, it is decided whether the pattern develops spatially. Due to random initial placement, an identical pattern does not occur, but qualitatively similar patterns can be stably generated; (3) When several states 1 are gathered, complicated pattern changes occur in the first two or three steps due to the balance between Rules 3 and 4, and if a shape that is not point symmetric is subsequently generated, a pattern spreading in the space is generated. When a point symmetric shape is born, a circular or ring-fixed pattern is generated; (4) When the black (state 1) area is large, the black area disappears due to Rule 1.
In a future study, if CA rules can translate to a chemical reaction process, a control technology for nanomachines becomes feasible. Thus, the simple totalistic CA presented in this paper allows the emergence of a wide range of self-organization patterns to be observed.  Based on these results, if the parameters w 1 and w 2 are changed at each time step, it is possible to generate a continuously changing pattern. Figure 20 shows a calculation case in which w 1 and w 2 are changed in time over time. The video of Figure 20 can be viewed in the Video S2 of Supplementary Materials. As in this result, it is possible to create changes by changing w 1 and w 2 ; birth of annulus -> growth -> division -> move -> decline -> growth -> annihilation. Figure 21 shows the transition of the ratio of states 1 and 2, and the ratio of the applied cells of Rule 4 in each time step during the calculation of Figure 20. Rule 4 was applied more frequently when spot patterns replicated, as compared to when the spot moved or vibrated. Furthermore, when the pattern disappeared, the application frequency of Rule 4 increased. Thus, it is suggested that pattern stability can be controlled by changing the application ratio of Rule 4.
Since this research is still in its first stage, discussions are limited mainly to qualitative considerations, but in future work, it will be necessary to investigate further quantitative assessments. There is a possibility of finding an index like Langton's λ parameter with the evaluation of many calculation cases. Furthermore, since it has rules similar to the Game of Life, basically like the Game of Life, many trials are necessary to find the desired pattern. Slight changes in the initial value may develop into different patterns. In the future, it is necessary to make a catalog of pattern formation on the relation between the initial condition and patterns based on many calculation results.  Figure 20. Rule 4 was applied more frequently when spot patterns replicated than when the spot moved or vibrated. Furthermore, when the pattern disappeared, the application frequency of Rule 4 increased. Thus, it is suggested that pattern stability can be controlled by changing the application ratio of Rule 4.

Conclusions
To observe special patterns of emergent behavior that mimic biological systems, we developed a totalistic CA model with an activating inner area and inhibiting outer area. This is a CA model, yet it works in a way that is equivalent to that of RD approaches like the Turing pattern model. If a multiple CA model can be designed that allows the number of states to be increased, we may witness the emergence of more complex patterns, like those of living things.
Our model also has potential applications in the engineering field. For example, it may be possible to develop swarming robots that can work as a single unit or divide into groups. As shown in the result of Figure 13, spot-like regions can be formed in the space by this algorithm. If many robots are dispersedly arranged in this space, it becomes possible to gather by grouping the robots within this spot-like rounded area. Furthermore, using this algorithm, it is also possible to divide one spot area into two. Therefore, robots can be divided into two groups simply by paying attention so that the robot that was in the first spot area does not come out of the area.
The model can further be applied to information networks. If the transition rules can be applied to Internet of Things (IoT) devices through the Web, hierarchical structures of IoT devices may emerge. For example, if the IoT devices know their own position, this algorithm can be applied by transmitting tokens to the nearby devices. Then, the spotted structure can be formed in the network  Figure 20. Rule 4 was applied more frequently when spot patterns replicated than when the spot moved or vibrated. Furthermore, when the pattern disappeared, the application frequency of Rule 4 increased. Thus, it is suggested that pattern stability can be controlled by changing the application ratio of Rule 4.

Conclusions
To observe special patterns of emergent behavior that mimic biological systems, we developed a totalistic CA model with an activating inner area and inhibiting outer area. This is a CA model, yet it works in a way that is equivalent to that of RD approaches like the Turing pattern model. If a multiple CA model can be designed that allows the number of states to be increased, we may witness the emergence of more complex patterns, like those of living things.
Our model also has potential applications in the engineering field. For example, it may be possible to develop swarming robots that can work as a single unit or divide into groups. As shown in the result of Figure 13, spot-like regions can be formed in the space by this algorithm. If many robots are dispersedly arranged in this space, it becomes possible to gather by grouping the robots within this spot-like rounded area. Furthermore, using this algorithm, it is also possible to divide one spot area into two. Therefore, robots can be divided into two groups simply by paying attention so that the robot that was in the first spot area does not come out of the area.
The model can further be applied to information networks. If the transition rules can be applied to Internet of Things (IoT) devices through the Web, hierarchical structures of IoT devices may emerge. For example, if the IoT devices know their own position, this algorithm can be applied by transmitting tokens to the nearby devices. Then, the spotted structure can be formed in the network space of devices, and the devices can be divided into plural groups. Furthermore, if our method can be applied to microsystems, it may become possible to make a production process of self-organizing a microrobot from nanoparts automatically.
The model also produces Turing-like patterns using a simple algorithm and without the need to solve complex simultaneous differential equations. However, a weak point is the necessity to add up the state quantities within a certain area. While these calculations are simple to perform using electronic devices, it is difficult to see how an individual living cell in an organism would be able to calculate the states of distant cells. To apply this model to morphogenesis, therefore, further modification is needed.
Other outstanding questions include the mathematical relationship between this model and the RD equations, and the relationship with a Turing machine. It would also be interesting to consider the entropy of the model.
Funding: This research was supported by grants from the Project of the NARO (National Agriculture and Food Research Organization) Bio-oriented Technology Research Advancement Institution (the special scheme project on regional developing strategy).

Conflicts of Interest:
The authors declare no conflict of interest.