Abstract
The Turing pattern model is one of the theories used to describe organism formation patterns. Using this model, self-organized patterns emerge due to differences in the concentrations of activators and inhibitors. Here a cellular automata (CA)-like model was constructed wherein the Turing patterns emerged via the exchange of integer values between adjacent cells. In this simple hexagonal grid model, each cell state changed according to information exchanged from the six adjacent cells. The distinguishing characteristic of this model is that it presents a different pattern formation mechanism using only one kind of token, such as a chemical agent that ages via spatial diffusion. Using this CA-like model, various Turing-like patterns (spots or stripes) emerge when changing two of four parameters. This model has the ability to support Turing instability that propagates in the neighborhood space; global patterns are observed to spread from locally limited patterns. This model is not a substitute for a conventional Turing model but rather is a simplified Turing model. Using this model, it is possible to control the formation of multiple robots into such forms as circle groups or dividing a circle group into two groups, for example. In the field of information networks, the presented model could be applied to groups of Internet-of-Things devices to create macroscopic spatial structures to control data traffic.
1. Introduction
The Turing pattern model was introduced by Alan Turing in 1952 [1] to account for morphogenesis (a developmental process of patterns found in nature) via interactions between activating and inhibiting factors. This typical self-organization model is described by reaction–diffusion (RD) equations that use different diffusion coefficients for two morphogens that are equivalent to the activating and inhibiting factors. The general RD equations can be written as
where and are the morphogen concentrations, the functions and represent the reaction kinetics, and and are the diffusion coefficients. Various functions of and have been previously studied and several models, such as the linear model, the Gierer–Meinhardt model [2], and the Gray–Scott model [3], produce typical Turing patterns.
The first observation of Turing pattern formation in nature was conducted by Kondo and Asai in a study of stripes on the marine angelfish Pomacanthus [4]. The two factors responsible for this skin pattern, however, remain to be identified. Kondo reported that cells transfer signals through the pseudopodium, rather than using chemical agents [5]. Numerous studies have proposed candidate substances related to Turing patterns, however, there is no scientific proof that these patterns are related to the Turing pattern model [6]. Some biological patterns have been described as containing features of the Turing pattern model. However, it is doubtful that all of an organism’s patterns can be described by the diffusion of two chemicals.
The Turing pattern model states that short- and long-range spatial scales are affected by separate factors [6]. Self-organized pattern formation is then caused by nonlinear interactions between these two factors. Turing considered two chemicals with different diffusion coefficients. However, as long as there is a difference between long- and short-range effects, other methods of implementation can be considered, such as chemical aging or protein modification.
It is unrealistic to believe that the entire shape of an organism can be determined by differences between diffusion coefficients, because factors comprised of similarly sized chemicals diffuse at similar speeds, and differences in diffusive concentrations are small. Therefore, sufficiently far from a morphogen source, the concentration gradient is effectively zero, at which point it becomes difficult to profit from concentration differences between chemicals elsewhere in a body.
In the case of fish stripes, however, the epidermal cells of a fish originate from nerve cells and these cells are connected to each other through the pseudopodium. These epidermal cells are believed to measure distance by means of stimuli through the pseudopodium [6]. However, it is difficult to apply such a mechanism to other morphogenetic mechanisms inside the body of other organisms.
A morphogen is defined as a signaling molecule that acts directly on cells to produce specific cellular patterns dependent on the concentration of the morphogen. The Wolpert model [7, 8] is a standard conceptual theory of a biological morphogen. However, the drawback of the Wolpert model is that the morphogen source is fixed. Therefore, this model cannot describe the revival phenomenon wherein part of an entire body, which has been divided into multiple parts, regrows into an entire body, such as that observed in a planarian. The Turing pattern model can be used as a method to compensate for this drawback. Green and Sharp stated that the Wolpert model was a special case of the Turing pattern model [9].
Here, a new cellular automata (CA)-like model was constructed to create Turing patterns via the exchange of integer values between adjacent cells. I believe that this is a basic model of short-range activating and long-range inhibiting factors. In general, CA models are discrete models in both time and space wherein the new state of each cell is determined by the current states of its neighbors. A state transition is based on previously specified transition rules. Here, the grid cells on the simulation space are referred to as cells and real cells in the organism are referred to as living cells.
Many studies have applied CA models to biological phenomena, such as the generation of seashell-like structures by Wolfram [11], the development of an RD model in a CA setting that reproduced the patterns of animal coloration by Young [12], the application of the Belousov–Zhabotinsky reaction CA model by Madore and Freeman [13], the catalytic modeling of proteins by Gerhardt and Schuster [14], the modeling of the immune system by De Boer and Hogeweg [15] and Celada and Seiden [16], the tumor growth model of Moreira and Deutsch [17], the genetic disorder model of Moore and Hahn [18], the hippocampal modeling by Pytte et al. [19], and the dynamic modeling of cardiac conduction by Kaplan et al. [20]. These historic studies demonstrated the effectiveness of CA models and exhibited various in vivo phenomena in terms of local interactions.
The following two CA models related to the RD equation have been studied. Markus [21] and Schepers [22] demonstrated that a CA model could produce the same output as RD equations. The Markus model is comprised of complex rules that model the RD equations. Conversely, Young [12] created a simple method to construct Turing patterns on a cellular grid using a simplified distribution pattern of two morphogens.
The Young model is a two-dimensional (2D) outer-totalistic model. In an outer-totalistic CA, each cell state is determined by summing the states of the neighboring cells. Adamatzky [23] studied a binary-cell state, eight-cell neighborhood, 2D cellular automaton outer-totalistic model with semi-totalistic transitions rules to produce Turing-like patterns. Dormann [24] also used a 2D outer-totalistic model with three states to produce a Turing-like pattern. Tsai [25] analyzed a self-replicating mechanism in Turing patterns using a minimal autocatalytic monomer–dimer system. Adamatzky [26] made a Glider-based model to compute reaction–diffusion patterns using hexagonal CA.
The preceding examples are representative CA models related to the RD phenomenon. The studies [12, 23–26] all used outer-totalistic CA models. These are considered models that will succeed in reproducing the essential phenomena of life, such as self-replication.
Young’s CA model [12] uses a real number function to derive distance effects using two constants within a grid cell, (a positive value) and (a negative value). The function is a continuous step function dependent upon and , as shown in Figure 1(b). The variable is the distance from the focal cell. This function assigns the activation area to and the inhibition area to . The activation area has a radius of and the inhibition area has a radius of ( > ), as shown in Figure 1(a). There are two states, black and white, on a rectangular grid. The calculation begins by distributing the black cells randomly on a rectangular grid. Then, for each cell (at position ), the next cell state is calculated by summing the values at positions . The values of the function, are summed for each black cell within a distance from on the grid. Here, is the distance from to , and is the number of black cells within a radius from . If , the grid cell at point is marked as a black cell. If , the grid cell is marked as a white cell and if , the grid cell does not change state [12]. Young showed that a Turing pattern could be generated using these functions. In the Young model, the parameters , , , and are chosen independently to find the combination of parameters that produce Turing patterns.
(a)
(b)
These 2D outer-totalistic models, including Young’s model, support Turing patterns via a simple algorithm, without the need to solve complex simultaneous differential equations. However, the drawback of this model is the necessity to sum the state quantities within a certain area. These calculations are simple to implement with computer processing. However, an individual living cell in an organism cannot integrate the states of distant living cells by itself.
To overcome this drawback, a CA-like model was constructed to support Turing patterns via the exchange of integer values between adjacent cells. In this model, each cell state changes according to information from the adjacent cells only; information of more distant cells is not used directly. In addition, the distinguishing characteristic of this model is that it presents a different pattern formation mechanism that uses only one kind of token, like a chemical agent that ages via spatial diffusion.
The proposed model was based on the exchange of integer signals between adjacent cells. Tokens labeled with integer values were moved between adjacent cells, as shown in Figure 2. The integer values were increased with each token movement. Each cell can contain multiple tokens simultaneously. Using these integer increments, the effects of both long and short distances could be implemented: a smaller integer indicates a shorter range, and vice versa. A detailed description of the model is given in the Methods.
In this model, the movement destinations of some tokens are determined, not by the local cell states but by the global simulation algorithm through stochastic token diffusion. This model is, therefore, not a strict CA model, in which the cell states are determined only by the cells, so the term CA-like model is used in this paper.
This CA-like model supports different Turing patterns (e.g., spotted or striped) for different parameter settings. Moreover, this model has the ability to support Turing instability that propagates in the neighborhood space; patterns are seen to spread from locally limited initial patterns. This model is not a substitute for the conventional Turing model but rather a simplified Turing model.
A distinguishing characteristic of this model is that it uses only one kind of token, like a chemical agent that ages via spatial diffusion. As indicated in Figure 3, there are various ways to use this feature. If the integer labels on the tokens are used to represent different concentrations, this becomes similar to the conventional model that uses two diffusive chemicals. If instead, they are used to represent signal levels, this is equivalent to the pseudopodium–nerve model for modeling the skin patterns of fish. Therefore, the presented model is a basic model that can be applied to various processes, such as concentration differences, signaling, chemical aging, and protein modification. The applicability of this model to protein modification is considered in the discussion.
These models help explain findings related to morphogens in living organisms and can be used to innovate the construction of engineering systems, such as multi-robot control and self-forming robots.
Another characteristic of this model is that each cell state changes based on the information from adjacent cells only. Using this model, it is possible to describe the process by which gap–junction intercellular communication in living cells contributes to morphogenesis. This is one form of intercellular communication between animal cells that directly connects the cytoplasm of two cells, allowing various molecules, ions, and electrical impulses to pass directly between the cells. Embryos with blocked gap junctions fail to develop normally. The links between antibodies blocking the gap junctions and morphogenesis are unclear but systematic studies have been undertaken to elucidate the mechanism [27, 28]. In this model, as well as in gap junctions, there is a common point of information transfer with neighboring cells. Furthermore, the algorithm of this model may lead to the elucidation of the mechanism of the morphogenesis of gap junctions. However, further research is required to understand the detailed relationships.
2. Methods
This model used 2D hexagonal grids wherein the transition rules are simple to apply. Square grids are commonly used in 2D CA, however, hexagonal grids were used here for the following reasons.
(a) In a square grid, the subsequent cell state is based on the states of the considered cell and its eight neighbors. There are only six neighbors in a hexagonal grid, which simplifies the number of transition rules to be considered.
(b) A hexagonal grid is isotropic, whereas a square grid is not. Since this model includes the process of distributing tokens to adjacent cells, it is simpler to apply when the distances between adjacent cells are equal. This model can also be applied to a square grid, but the pattern that is created is not isotropic. Many previous studies, such as those by Adamatzky [26] and Schepers [22] used isotropic grids.
The specifications of the model are as follows:
(a) Two-state cells are assumed, namely, black and white .
(b) Each cell can contain multiple tokens simultaneously. These tokens are labeled with integer values.
(c) The black cells are the source of the tokens. This is where the first parameter is set: The token’s production number is the initial number of tokens on a black cell at each time-step. This number is fixed for all cells and all time-steps.
(d) Some of the tokens are distributed to adjacent cells. The tokens are created only in black cells but they can move on either black cells or white cells. The residual token ratio is the second parameter, which is the fraction of unchanged (unmoved) tokens in each cell. The parameter is fixed for all cells and all time-steps. As shown in Figure 4, the following calculations are performed, where is the number of tokens in cell at time : (1) A proportion, of the tokens diffuse toward the six adjacent cells evenly and their labels are incremented by 1 upon arrival. (2) A proportion, , of the tokens remain in the original cell and their labels are incremented by 1. This is the process whereby the token leaves the original cell, then returns to the original cell again in one time-step. (3) The residual, tokens remain in the original cell without modification of their integer labels. If the number of tokens is not an integer multiple of 7, the remainder of the tokens are distributed between adjacent cells with equal probability.
As each token moves, its label value increases by 1 in each time-step. These modified tokens move to other adjacent cells and are modified again. In this process, some tokens remain in their original cells without modification. The incremental labeling in each time-step is repeated until the maximum label number reaches (where is a fixed, even integer). The above mentioned transfers are conducted on all grid cells. As the calculations progress, various numbered tokens are transferred over the grid. The integer value on the token indicates the number of movements; small integers indicate short travel and large integers indicate long travel. These values correspond to the short- and long-range spatial scales of the Turing model. If , all tokens travel long distances counting from the first time-step and there are no short travelers. Therefore, it is necessary that .
(e) The third parameter is the token dissipation rate . If there are fewer than tokens in a cell, those tokens are removed from the cell. At this time, the state of the cell does not change. This removal process controls the effective distance of the long-range signals.
In this model, if , spatial patterns do not emerge. This is because, if the tokens in cells with less than tokens are not eliminated from the simulated space, some tokens travel cells ahead. This distance is beyond the scale of a Turing wave and it becomes the disturbance factor in the pattern formation. Therefore, is an important parameter in this model.
(f) Processes (c), (d), and (e) are applied to all cells in the simulated space and the placement of tokens at the next time-step is determined.
(g) After the movement and modification of the tokens, the subsequent state of each cell is determined by the fraction of “younger” tokens that it contains. The number of tokens with labels less than indicates the short-range factor and the number of tokens with labels larger than indicates the long-range factor. The ratio of the number of tokens labeled equal to or lower than the total number of tokens in the cell are used in formula (3), from which the subsequent cell state (black or white) is determined. Figure 5 shows this process of cell state change.
Here,
= number of tokens labeled from 1 to .
= number of tokens labeled from + 1 to .
Parameter is equivalent to / in the Young model and controls patterns, such as spotted patterns or striped patterns. The results were calculated by varying the values of the parameters , , , and . Figure 6 shows the calculation flow.
It is impossible for self-organized patterns to emerge from the distribution of tokens alone. The tokens are distributed in a diffusion-like manner because they do not have an intention to move in a specific direction. Therefore, there is no spatial structure within the tokens’ configuration of the hexagonal grid. The essence of the present model is the transition function, formula (3), which uses the ratio of the numbers of tokens that have traveled short and long distances. In formula (3), the fourth parameter is identified. The next state is determined by the ratio of these two short and long factors.
For this model, a two-dimensional 100 × 100 hexagonal grid was used with periodic boundary conditions in which the border grids continued to the border of the other side without any gaps. The model was implemented in Java. Each grid cell had two associated variables: (1) the cell state (black or white) and (2) the number of tokens counted for each label.
3. Results
When the maximum number of modifications () is relatively small, the difference between the short-range effect and the long-range effect becomes small and spatial patterns do not emerge. When is approximately 24, spatial patterns are observed. Changing merely changes the pattern scale. As increases, the relevant pattern simply becomes larger. The same effect can be achieved by enlarging the activator and inhibitor areas in Young’s model. Therefore, (the maximum number of modifications) was set to 24 for all the results below.
In this study, four parameters, , , , and , were examined. There are various combinations of these four parameters. After a large number of trial simulations, the parameter ranges that caused an emergence of spatial patterns were found empirically. Within the parameters ranges, the parameters which emerged as the Turing-like spotted pattern were:
The number of tokens produced : 10,000; The residual ratio of the morphogen : 0.05; The dissipation rate of the morphogen : 0.008; and The weight coefficient : 0.60.Four black cells were placed in the central simulation field as an initial condition (Figure 7). Figure 7 shows the formation process of a spotted pattern under the conditions , , , and . The initial four black cells became a more larger single spot at the center ( in Figure 7) and the spot divided in two ( in Figure 7), as did each subsequent spot. Finally, the simulated space was covered with a spotted pattern. This formation process was more similar to the Gray–Scott model than the Turing pattern. However, this process was similar to typical Turing instability which generates heterogeneous spatial patterns. The model had the ability to propagate central small instabilities in the neighborhood space. Tsai [26] reported similar results.
Figure 8 shows the results of the model for various values of parameters and . Four black cells were placed in the simulation field as the initial condition. For , self-organized patterns emerged over a wide range of . Spotted patterns appeared at and striped patterns appeared at For various spatial patterns appeared. For unique patterns could only emerge over a limited range of , indicating that the dissipation rate was a sensitive parameter. Figure 8 only shows the results of discrete parameter values of and . However, there were threshold parameter values for shifting from striped patterns to whole black color, and threshold parameter values at which a fixed spotted pattern started to spread in the field. For and , the reaction wave stopped after a certain number of steps and a steady circle emerged. Whether this phenomenon was observed was dependent on the value of , which suggests the possibility of pattern control.
As described in the Methods section, is the weighting parameter of the short-distance effect and the long-distance effect. Because parameter determines the number of long travel tokens eliminated beyond a certain distance, can partially control the long-distance effect. Therefore, partially controls the patterns, such as the spotted or striped patterns, under the same value.
Figure 9 shows the effect of varying the residual morphogen ratio r and weight factor for and . Spatial patterns appeared for The residual rate had the same effect as , i.e., it could control the transition from spotted patterns to striped ones. The parameter controls the amount of tokens within the short-range effect. Therefore, the parameter can change the patterns, as can the parameter w.
Figure 10 shows the results of varying the token production number for , , and . The number of the source tokens in each black cell does not affect the simulation results for This is because the transition function depends on the token ratio and not the absolute number. For example, if is increased by a factor of 10, the values on both sides of formula (3) increased by 10.
In the case of , this suggests a shortage of tokens to count the long and short effects. This result also shows the possibility of creating patterns near .
These results were based on initial conditions whereby four black cells were placed in the simulation field. If only one black cell was used instead, Turing instability did not occur and the one-spot steady state remained. For more than two initial black cells, Turing instability spread within the grid field. The reasoning here is that one black cell causes perfectly symmetric token configuration, which prevents the formation of Turing instability.
4. Discussion
The present model used only the time variation of integer values and not the difference between the diffusion coefficients of two chemicals. This model can control spatial patterns using one parameter or under the condition whereby the other two parameters are fixed. However, the present striped patterns obtained are thought to differ from patterns that emerge from the typical Turing model. This model has the ability to propagate Turing instabilities in the neighboring spaces, which are observed by the spread of the spotted pattern from the central limited area to the entire field by the dividing processes. Conversely, the Young model requires the placement of multiple state 1 cells randomly on the simulated space.
When numerically solving the RD equations, various patterns are produced according to the reaction terms. The present model is not an alternative to the Turing pattern model but rather is one of its sub-models. Further study is required to find RD equations equivalent to the present model, the main feature of which is the highly self-propagating ability of unstable waves. If the initial pattern is comprised of more than two black cells, self-organized spatial patterns spread in the simulated space, as in the Gray–Scott model. The Gray–Scott model can create a process in which spotted patterns are continuously replicated under certain conditions. The present model created the Gray–Scott like patterns under certain conditions. However, further studies are needed as to how the present model relate to the Gray–Scott model mathematically.
It is possible to control the pattern easily using the present model, such as creating a spotted one. Applying this algorithm to robotic communication could control the formation of multiple robots, such as forming circle groups or dividing a circular group into two groups. As shown in Figure 7, this algorithm can form spot-like regions in the space. If multiple robots are dispersedly arranged in this space, it is possible to group the robots within this spot-like rounded area. Furthermore, using this algorithm, it is possible to divide a spot area into two halves. Therefore, robots can be classified into two groups simply by maintaining their position such that the robots that are in the first spot area do not leave that area. In the information network field, the current model could be applied to groups of Internet-of-Things devices to create macroscopic spatial structures or groups to control data traffic. For example, if the Internet-of-Things devices are aware of their own positions, this algorithm can be employed by transmitting tokens to the neighboring devices. Then, a spotted structure can be formed in the network space of the devices and these devices can be divided into several plural groups.
One of the potential processes that the algorithms can be applied to in living cells is proteins that gradually degrade and change their structure. That is, it is a process in which the next reaction is determined from the ratio of the numbers of proteins with a certain functional structure. Considering the modification of proteins (post-translation modification) during inter-cell processes, the present model could produce Turing patterns using one chemical aging agent instead of two diffusive substances. No previous CA model has addressed the aging effect that produces spatial patterns. It is envisaged that one protein could act as a signaling factor that is released into the intercellular space and gradually modified on the surface of neighboring cells. As this factor becomes more diffused, it could change the shape of the protein. The difference in the shape of the protein can be correlated to the number of tokens labeled less than and the number of tokens labeled larger than . This modifying protein is observed as an extracellular protease in real organisms. If this factor could perform more than one modification, the protein that is modified more than a certain number of times would change the signaling message inside the cell. These processes produce two types of factors virtually. The diffuse morphogens are modified on neighboring cell membranes in a manner similar to the mechanism of proteolysis. Once the protein has spread sufficiently, each cell membrane has various amounts of the two factors, which depend on the distance from the original source cell. The effect of distance can be measured using the ratio of these two modified protein factors. This mechanism can produce a form of spatial pattern. However, further studies are needed into how the proposed algorithm can model such protein-based communication processes.
Therefore, the algorithm presented by the model can mimic similar mechanisms (e.g., the cleavage of peptide chains) seen in living cells. It could, therefore, contribute to identifying the real morphogen in the morphological formation process. Concrete candidate substances were not specified here, however, and should be researched in the future. Application of a proteolysis model could contribute to the control of morphogenesis. Proteins are modified and used as signals between cells and future studies are needed to address how this model can be applied to those biochemical processes.
As indicated in Figure 3, the study model can be used to address various processes, such as chemical concentrations, signals, and aging effects. Furthermore, extended models using several types of tokens can be considered. This model may be able to explain various chemical reactions by the conversions of tokens. In addition, an extended model could potentially express both the chemical reaction and the self-organization phenomena demonstrated in this study.
Research evaluating the stability of this model will be necessary in the future. As in Wolfram's one-dimensional CA model, the model presented here may be able to form steady patterns, periodic patterns, or chaotic patterns depending on the choice of parameters and initial values. Regarding this point, much further calculation and consideration is needed and the results will be summarized in a future paper.
5. Conclusions
Herein, a model based on a simple mechanism was introduced, which showed the possibility of pattern formation evoked by cell division or tissue growth. The present model was developed to combine the Turing pattern model with a pathway model for biochemical reactions in a living cell. In the future, a rigorous mathematical description and a concrete in vivo example are needed. If BioBrick parts, which are DNA sequences that conform to a restriction-enzyme assembly standard using iGEM, are combined with this token-based CA-like model, various bio-robot systems may be designed in a manner similar to the CAD systems of mechanical design. It is hoped that the present model will lead to a deeper understanding of the morphogenesis process and contribute to regenerative medicine. Future plans for this token-based model include constructing a 3D model, a circulation model (considering the circulation of tokens), and an energy model (considering the energy budget of the transition rules).
Data Availability
No data were used to support this study.
Conflicts of Interest
The author declares no conflict of interest regarding the publication of this paper.
Funding
This research was supported by Grants from Japan Society for the Promotion of Science, KAKENHI Grant Number 19K04896.
Acknowledgments
The authors would like to thank Enago (https://www.enago.jp/) for the English language review.