Formation of morphogenetic patterns in cellular automata
FFormation of morphogenetic patterns in cellular automata
Manan’Iarivo Rasolonjanahary and Bakhtier Vasiev ∗ Department of Mathematical Sciences, University of Liverpool, Liverpool, UK ∗ e-mail: [email protected] Abstract.
One of the most important problems in contemporary science, and especially in biology,is to reveal mechanisms of pattern formation. On the level of biological tissues, patterns form due tointeractions between cells. These interactions can be long-range if mediated by diffusive molecules orshort-range when associated with cell-to-cell contact sites. Mathematical studies of long-range interac-tions involve models based on differential equations while short-range interactions are modelled usingdiscrete type models. In this paper, we use cellular automata (CA) technique to study formation ofpatterns due to short-range interactions. Namely, we use von Neumann cellular automata representedby a finite set of lattices whose states evolve according to transition rules. Lattices can be considered asrepresenting biological cells (which, in the simplest case, can only be in one of the two different states)while the transition rules define changes in their states due to the cell-to-cell contact interactions. Inthis model, we identify rules resulting in the formation of stationary periodic patterns. In our analysis,we distinguish rules which do not destroy preset patterns and those which cause pattern formation fromrandom initial conditions. Also, we check whether the forming patterns are resistant to noise and analysethe time frame for their formation. Transition rules which allow formation of stationary periodic patternsare then discussed in terms of pattern formation in biology.Keywords: pattern formation, mathematical modelling, cellular automata.
Biological pattern formation is a complex phenomenon which is studied experimentally in a numberof model organisms [1] and theoretically by means of various mathematical techniques [2, 3]. One ofthe model organisms commonly used in biological studies of pattern formation is the fly (Drosophila)embryo [4]. Patterning in the fly embryo takes place in two directions: along the head-to-tail (antero-posterior) axis where the pattern occurs as a repeated structure formed by segments and along thedorso-ventral axis where patterning results into formation of internal morphological structures. It isknown that the formation of segments is preconditioned by the formation of spatially periodic patternsof gene expressions which takes place at early stages of embryonic development [5]. There are many otherexamples of formation of periodic patterns in biology including formation of stripes on skin of animals(zebra) and fish.In this work, we focus on formation of periodic patterns which allow discrete representation. Discretemodels often complement continuous models by giving description of patterning processes on differentspatial or time scales. For example, the continuous (reaction-diffusion) model of anterio-posterior pat-terning in fly embryo presented in [6] allows to conclude that the modelled units (nuclei) quickly reachtheir stable equilibrium states and there are only four such states in normally evolving embryo. This,in turn, allows to simplify the model by reducing its quantitative details and to describe the consideredpatterning as evolution of a chain of interacting nuclei on discrete level. Another example is given bycontinuous model of pigmentation patterning on skin of growing reptilian [7]. Due to slow diffusion be-tween skin scales this model also reduces to discrete model for interaction between differently pigmented1 a r X i v : . [ n li n . C G ] J u l cales, although the pigmentation of each scale is described by differential equations. Generally, discretemodels, as compared to continuous models, often offer a simple explanation for basic properties (such asspatial scaling) of morphogenetic patterns.Since we are primarily interested in formation of periodic patterns, we can reformulate the problemand ask the general question: what kind of local interactions can result in formation of a stationaryperiodic pattern? In our study, we use von Neumann’s cellular automata (CA) to address this question.This model has been invented, more than fifty years ago, by Stanislaw Ulam and John von Neumann [8,9].The model is represented by a collection of cells forming a regular grid which evolves over discrete timesteps according to a set of rules based on the state of the cells [10]. A detailed description of Neumann’sCA can be found in a number of sources [11, 12].Extensive study of various modifications (involving different numbers of dimensions, allowed statesper site and neighbours affecting transitions) of this model have been performed by Wolfram [12–14].He has found that patterns forming in cellular automata fall into one of the following four classes: 1.Homogenous, 2. Periodic structures, 3. Chaotic structures and 4. Complex. Correspondingly, all automataalso fall into four classes, although the class in which the given automata falls depends on the imposedinitial conditions.In this work, we will explore rules resulting in the formation of periodic patterns in the CA whereeach cell can only be in two distinct states. This is an extension of the research reported in [12,13] in twoways: firstly, we consider all 256 rules (while the classification made by Wolfram was restricted by 32 socalled “legal” rules) and, secondly, we are interested in periodic patterns forming in a domain of finitesize, and therefore the found rules appear to fall into classes 2 and 4 identified by Wolfram. Although thismodel is very simple, the results obtained allow conclusions to be drawn on the mechanisms of periodicpattern formation in biology, for example the anterio-posterior patterning in fly embryo. The simplest version of von Neumann’s cellular automata was designed by Wolfram [12, 13] and called“elementary cellular automaton” (ECA). It consists of a regular lattice of cells which form a chainand can only be in two states. The position of a cell in the chain will be denoted by the symbol i ,the total number of cells - by n and the state of cell by s i ( s i = 0 or s i = 1). The state s t +1 i of acell i at time t + 1 is determined by the states of cells i − i and i + 1 at time t according to atransition rule, defined as s t +1 i = (cid:0) s ti − , s ti , s ti +1 (cid:1) . A cell interacting with its two closest neighbours leadsto consideration of 2 = 8 possible configurations, which are 000, 001, 010, 011, 100, 101, 110 and 111.For each of these 8 configurations, the resulting states can be represented as: 0 s
0, 0 s
1, 0 s
0, 0 s s
0, 1 s
1, 1 s s
1. The complete set of middle elements of all resulting states forms a binarynumber s s s s s s s s which varies between 0 and 255. In other words, the configurations aboveallow 2 = 256 possible rules or cellular automata. The binary number, translated into decimal, is takento be the rule number.Different rules applied to different initial conditions can result into formation of various nontrivialspatio-temporal patterns including oscillations and propagating waves. Among these 256 rules, a fewsymmetry groups can be distinguished so that the rules belonging to one group result in the formation ofsimilar patterns. These groups contain up to four rules and their symmetry-based relationships termedcomplement, mirror image and mirror complement [14]. The definitions of these relationships are basedon certain features of the binary representation of rules, namely:2he complement of a given rule is a rule whose binary representation is obtained from the binaryrepresentation of the given rule where (a) the digits are taken in the reverse order and (b) each digit isinterchanged with the opposite digit.The mirror image of a given rule is the transition rule such that if the given rule contains an elementaryrule f ( b, c, d ) = a then the mirror image contains an elementary rule f ( d, c, b ) = a , that is, the sameoutput is allocated to the rule with mirrored states of neighbours.The mirror complement is a complement to a mirror image.Our study is performed using the one-dimensional CA described above. Note that this model can alsobe seen as a chain of logical elements. It is represented by a onedimensional regular lattice of elementswhich can be in one of two states. The state S i of the element i , can be either “false” (commonly denotedby “0” and illustrated in black colour) or “true” (denoted by “1” and drawn in white). The state of eachelement changes with each time increment (except for the border cells which remain unaltered) accordingto the applied rules.The goal of this study is to find out what rules can cause the formation of stationary periodic patternsin a chain of finite size where the states of border elements are fixed. We also check whether the periodicpatterns are stable, i.e. noise sensitive. To do so, we introduce a noise by changing the state of eachelement randomly with a certain (small) probability. In application to biology, each logical elementcorresponds to a biological cell and the cell’s state reflects the expression of certain genes. The ECArepresents a two-state model and can only be used for study of patterns formed by expression of a singlegene along cells forming a chain. Our study of periodic patterns forming in the ECA has been performed in a few steps. First, we havefocused on two-periodic patterns and identified the rules which (1) do not destroy a preset two-periodicpattern, (2) allow recovery of two-periodic patterns perturbed by noise and (3) allow formation of two-periodic patterns from random initial conditions. Then, we have repeated the above study for the caseof three- and more- periodic patterns.
Our first set of simulations aimed to find all transition rules which do not destroy pre-existing periodicpatterns. These simulations have been started with preset two-periodic patterns (1 white + 1 black) likethe one shown in Fig 1 but with more (usually 60 or more) elements. Starting with two-periodic initialconditions, we have found that 64 rules out of 256 (25%) do not destroy it. Their values, in decimal form,are as follows:4, 5, 6, 7, 12, 13, 14, 15, 20, 21, 22, 23, 28, 29, 30, 31, 68, 69, 70, 71, 76, 77, 78, 79, 84, 85, 86, 87, 92,93, 94, 95, 132, 133, 134, 135, 140, 141, 142, 143, 148, 149, 150, 151, 156, 157, 158, 159, 196, 197, 198,199, 204, 205, 206, 207, 212, 213, 214, 215, 220, 221, 222 and 223.3 igure 1:
Two-periodic pattern in the ECA composed by 20 cells. The binary expression for all these rules hasthe form: xx xx xx , where x can be either 0 or 1. That is, the pattern consisting of a repetition of black andwhite remains the same when the transition rule leaves the configurations 3 and 6 unaltered: (010) → (0
0) and(101) → (1 The second set of simulations was aimed at finding all rules which are such that, starting with periodicinitial conditions (see Fig 1) and applying a noise, the perturbed pattern is able to recover. The noiselevel of 0.1% has been applied meaning that the state of each cell can be altered with probability of0.001 (i.e. one out of 1000 cells is altered each time step). Simulations show that 33 out of the above 64rules (52%) would allow resistance to the perturbation caused by the noise. These rules are given by thefollowing numbers:6, 13, 14, 15, 20, 28, 30, 69, 70, 77, 78, 79, 84, 85, 86, 92, 93, 134, 135, 141, 143, 148, 149, 156, 157, 158,159, 197, 198, 199, 213, 214 and 215.The recovery of a perturbed pattern can happen in two ways: locally, when the perturbation dis-appears without affecting surrounding cells, or globally, when the perturbation propagates along themedium to one of its borders and then vanishes. Also, for both scenarios there are exceptional rules (28,70, 78, 92, 141, 156, 157, 197, 198 and 199) such that if the perturbed cell is located next to the bordercell then it cannot recover (unless it is hit by another perturbation).
The following rules allow the local extinction of perturbations:13, 28, 69, 70, 77, 78, 79, 92, 93, 141, 156, 157, 197, 198 and 199.For the analysis of the recovery processes, we point out that if the probability of perturbation issmall, then the probability that two neighbouring cells are perturbed simultaneously is much smallerand, therefore, can be neglected. With this in mind, let us consider the recovery of the periodic patternin a chain described, for example, by the rule 77, whose binary number is 01001101. Consider the periodicsequence . . . . . . and the scenario when, for example, the element in the 4 th position changes to due to the perturbation caused by noise. Then, the above sequence will become . . . . . . . Wewant to know how this sequence will evolve over time and how many time steps it will take for thepattern to recover. The changed element can only affect its two closest neighbours, i.e. the zeros in 3 rd and 5 th positions. The 3 rd zero is in the middle of 10 and will not change following the transition(10 ) → (10 ). The perturbed cell is in the middle of 0 followingthe transition (0 → (0 th zero is in the middle of
01 and will remain unaltered followingthe transition ( → ( . . . . . . ) in a medium of arbitrary (finite) size. It is easy to show that if the noisechanges any of the 1s to 0 in the periodic sequence it will take only one time step to recover. On theother hand, it takes two time steps to restore when the noise changes any of 0s to 1. To illustrate this,let us consider the finite periodic sequence . . . . . . . where the noise changes the digit 0 atthe 3 rd position into 1, so that the sequence becomes . . . . . . . The 1 in the second positionwill stay unaltered following the transition: (011) → (011). The 1 in the third position will change to0 following the transition (111) → (101). The 1 in the fourth position will change to 0 following thetransition (110) → (100). So, after the first step we get the sequence . . . . . . . This is identicalto the case (considered above) when the state of the (cid:0) th (cid:1) cell in periodic pattern is changed from “1”to “0”. Its recovery takes only one time step and thus it takes two time steps for the pattern to recoverwhen, owing to the noise, any of 0s change to 1. The stated result is exactly the same for the rule 69,which is the mirror image of the rule 13, and it is slightly amended for the rules 79 and 93, which arecomplements of rules 13 and 69 respectively. Namely, in the cases of rules 79 and 93 it takes one timestep to recover when any “0” is changed to “1” in the preset periodic pattern and two time steps whenany “1” is changed to “0”. Rule Complement Mirror image Mirror complement77 → →
13 79 69 93 → → → → → → → →
70 157 28 199 → → → → → → → →
78 141 92 197 → → → → → → → →
156 198 156 → → → → Table 1:
Summary of recovery dynamics for rules allowing local recovery of perturbed two-periodic pattern. Thenotation “ a → b ( n TS)” is used to state that the perturbed by noise cell recovers from state “ a ” to its originalstate “ b ” in n time steps. The exceptional cases, when the recovery doesn’t take place, (and when the noisehits the cell next to the border cell) are described in the extra (last) line. For example, the notation “(. . . 000)”indicates that the rule doesn’t allow recovery in the case when the perturbed chain ends with three 0s. → (100), which allows for the 3 rd cell from the right border remaining unaltered, and(000) → (000), for which the state of the second (perturbed) cell remains unaltered. The only way forthe periodic pattern to recover is that the noise strikes the same cell once again. The results for rule 70can, with small modifications, be extended to rules 157, 28 and 199 which are respectively complement,mirror image and mirror complement of rule 70. The properties of all rules which allow local recovery ofperturbations on preset two-periodic patterns are given in Table 1.The rules listed in Table 1 allow local recovery from any noise, namely from 0 to 1, as well as from 1to 0 perturbations. However, there are rules which allow the recovery of periodic patterns only from onetype of perturbation but NOT from the other and which have not been included in Table 1. These arerules 5, 76, 94, 95, 133 and 205. Rules 76, 94 and 95 allow the recovery of periodic pattern when, due tothe noise, 0 has changed to 1. Rules 5, 133 and 205 allow the recovery in case when, due to the noise, 1has changed to 0. All six rules, when they allow the recovery of perturbed pattern, do it in only 1 timestep. In this section, we consider a set of rules which also allow the recovery of perturbed periodic patterns butin a different manner: namely, the perturbation is driven to one of the two edges of the chain where itcommonly disappears. In the most common scenario, the perturbation moves like a wave with a constantspeed and the time required for the perturbation to get to the chain edge is proportional to the initialdistance from the perturbed cell to the corresponding edge of the chain. However, for some rules thewaves of perturbation which form are not regular, so that their speed and the size of perturbed areachange over time. The rules exhibiting propagating waves of perturbation are:6, 14, 15, 20, 30, 84, 85, 86, 134, 135, 143, 148, 149, 158, 159, 213, 214 and 215.The perturbation waves propagate from LEFT to RIGHT in the case of rules 6, 14, 15, 30, 134,135, 143, 158 and 159 and from RIGHT to LEFT for the remaining rules. The most important questionabout the propagating waves of perturbations is at what speed they propagate and, consequently, howlong it takes for the perturbation to reach the chain’s border where it usually vanishes. Our study showsthat for rules 15 and 85 (which form a symmetry group, see Table 2), the perturbed cell recovers inone time step and at the same time the perturbation shifts one cell (to the right/left in the case of rule15/85 respectively). Thus, the size of the wave is 1 cell and the speed of the perturbation wave is 1cell/time-step (see Fig 2) and, therefore, it takes no more than n time steps (where n is a size of thechain) for the perturbations to reach the chain edge and disappear.Properties of the perturbation waves forming under rules 14, 143, 84 and 213 (form a symmetrygroup, see Table 2) are similar to what was observed for rules 15 and 85. However, now it takes twotime steps for each perturbed cell to recover (see Fig 2B). The recovery (that is, the contracting borderof the perturbed area) as well as the perturbation (that is, the expanding border of the perturbed area)6hifts one cell per timestep and the number of perturbed cells (i.e. the size of the wave) is 2 cells. Whenthe perturbation reaches the edge of the chain, it doesn’t always disappear: for example, in the case ofrule 14 the state of the second rightmost cell remains perturbed (at “0”) if the state of rightmost cell inthe initial periodic pattern is 0 (the chain in the stationary state has 000’ on its right border). Similarexceptions occur for the other three rules in the symmetry group (see Table 2). Figure 2:
Illustrations of propagating perturbation waves in recovering patterns. Each row shows a fragment(perturbed) of initially periodic pattern. Rows (top to bottom) show the states of the fragment after each consec-utive time step. Perturbed cells are highlighted in green/blue to distinguish between the initial perturbationsfrom state “1” to state “0” or from “0” to “1” respectively. A: Each perturbed cell recovers in one time step. B: Each perturbed cell recovers in two time steps.
C: Green case - perturbed “0” recovers in one time stepand perturbed “1” - in two.
Blue case - perturbed “0” recovers in five and perturbed “1” in six time steps.
D:Green case - same as on Panel C.
Blue case - perturbed “0” recovers in nine and perturbed “1” in ten timesteps. E: Recovery is slow and completely irregular.
The recovery of the periodic pattern under rule 6 is more complicated and takes place in one of twoscenarios (Fig 2C). The first scenario takes place when the originally perturbed cell was in the state “1”and the perturbation changes its state to “0”. In this case, the perturbation expands with the speed of 1cell/time-step and the recovery takes place in a way that two cells recover simultaneously every secondtime step so that the size of the wave alternates between 1 and 2 cells. This is because it takes two timesteps to recover from “0” to “1” and one time step to recover from “1” to “0”. However, if the originallyperturbed cell was in state “0” then the recovery is much slower. After the initial transition process, theperturbation wave stabilizes, so that each perturbed 0 recovers in five time steps while each “1” recoversin six. The speed of the perturbation expansion is still 1 cell per time-step, but its size alternates between5 and 6 cells. In similarity to the case of rule 15 the perturbation after reaching the edge of the chain,doesn’t always disappear but the second leftmost cell remains perturbed (at “0”) if the state of leftmostcell in the initial periodic pattern is 0 (the chain in the stationary state ends with “000” on the right).Similar scenarios are observed for the three other rules (159, 20 and 215) in the symmetry group (seeTable 2). The recovery of the periodic pattern under rule 134 is identical to that under rule 6 when theinitially perturbed cell was in state “1” (Fig 2D). A difference is observed when the state of the initiallyperturbed cell changes from “0” to “1”. After the initial transition process when the perturbation wavestabilizes, the recovery turns to be even slower: each perturbed “0” recovers in 9 time steps while each“1” in 10. The speed of the wave is still 1 cell per time-step but its size alternates between 9 and 107ells. This observation naturally extends to rules 158, 148 and 214 forming the symmetry group withrule 134.
Rule Complement Mirror image Mirror complementLEFT-to-RIGHT RIGHT-to-LEFT15 85
Speed of the wave - 1 cell/TS Same as 15.Size of perturbed area - 1 cell
14 143 84 213
Speed - 1 Speed - 1 Speed - 1 Speed - 1Size - 2 Size - 2 Size - 2 Size - 2(. . . 000) (. . . 111) (000. . . ) (111. . . ) →
0: Speed - 1 0 →
1: Speed - 1 1 →
0: Speed - 1 0 →
1: Speed - 1Size 1/2 Size 1/2 Size 1/2 Size 1/2(. . . 000) (. . . 111) (000. . . ) (111. . . )0 →
1: Speed - 1 1 →
0: Speed - 1 0 →
1: Speed - 1 1 →
0: Speed - 1Size - 5/6 Size - 5/6 Size - 5/6 Size - 5/6(. . . 000) (. . . 111) (000. . . ) (111. . . )
134 158 148 214 → → → → →
1: Speed - 1 1 →
0: Speed - 1 0 →
1: Speed - 1 1 →
0: Speed - 1Size - 9/10 Size - 9/10 Size - 9/10 Size - 9/10(. . . 000) (. . . 111) (000. . . ) (111. . . )
30 135 86 149
Same as 149. Same as 149. Same as 149. Speed of expansion 1.Recovery is slow and irregular.Size is unlimited.
Table 2:
Summary on the recovery dynamics for rules allowing formation of propagating perturbation waves.The direction, speed and size of the waves formed are given. If the properties of perturbation waves depend onthe form of original perturbation, they are stated separately: the notation “ a → b ” is used to indicate that theoriginally perturbed cell changed from the state “ a ” into the state “ b ”. The exceptional cases, when the recoveryis not complete (the perturbation wave hits the border where the perturbation remains), are described in theextra (last) line. For example, the notation “(. . . 000)” indicates that the recovery is not complete when theperturbed chain ends with three 0s. The recovery of the periodic pattern under rules 30, 135, 86 and 149 (forming a symmetry group) issignificantly different from what we have seen so far (see Fig 2E). The perturbed area under these rulesstill expands with the speed of 1 cell/time-step (i.e. under rule 30, the right border of the disturbed areashifts 1 cell/time-step) but the recovery is much slower: the contracting border (i.e. the left border underrule 30) of the perturbed area is slow, is affected by oscillations and its motion is completely irregular.As a result, the perturbed area quickly expands until it reaches the border of the chain, i.e. its size is onlylimited by the size of unrecovered area in the chain (i.e. the area to the right of the contracting borderunder rule 30). Although the contracting border moves much more slowly than the expanding one, itcan be shown that its speed can’t be less than 2 cells per 9 time-steps. The properties of propagatingperturbation waves under all the rules allowing their formation are summarized in Table 2.8 .4 Formation of two-periodic patterns from random initial condition
In the last two sections, we dealt with the rules which allow the recovery of preset periodic patternsperturbed by a rare and random noise. In this section, we shall focus on a subset of those rules whichnot only allow the recovery of perturbed periodic pattern but also the generation of periodic patternsfrom any initial conditions. These are the following six rules: 15, 30, 85, 86, 135 and 149.These six rules form two symmetry groups: (15, 85) and (30, 86, 135 and 149) (see Table 2). Theprocess of periodic pattern formation is significantly different for these two groups: under rules 15 and85 periodic patterns are generated considerably faster than under the rules forming the second group.
Figure 3:
Formation of periodic patterns from random initial conditions under rule 15.
Three examples of short simulations illustrating the formation of stationary periodic patterns underthe rule 15 are shown in Fig 3. One can see that periodic patterns form in the course of time, first arisingat the right border and then expanding to the left with a speed of 1 cell/time-step. For rules 15 and 85,one can formulate propositions like the one below:
Proposition for rule 15:
Starting from any initial conditions, the state of the chain evolves towardstwo-periodic pattern (010101. . . or 101010. . . ) developing from the LEFT to the RIGHT with a speed ofat least 1 cell/time-step. The final stationary pattern may end with “11” or “00” on its RIGHT edge.
Proof of this proposition is omitted. Proposition for rule 85 differs only in the directionality of theprocess: under the rule 85 the periodic pattern evolves from right to left. The formation of patterns underthe other four rules (30, 86, 135 and 149) is considerably different and in a line with the informationprovided in Table 2. For example, in the case of rule 30 the formation of a periodic pattern starts atthe right border and develops to the left. The speed at which the periodic pattern expands is irregularand alters from 4 cells per time-step to 2 cells per 9 time-steps. Furthermore, depending on the initialconditions, the final pattern may remain non-stationary and exhibit oscillations at the right border of thechain. The oscillating spot can contain one element, ( . . . ↔ ( . . . . . . ↔ ( . . . . . . → ( . . . → ( . . . → ( . . . → ( . . . So far, we have focused on rules allowing existence/formation of two-periodic patterns. In this section, weaddress the same questions with respect to three-periodic patterns. Two types of three-periodic patternscan be set in ECA: two black cells followed by a white cell or two white cells followed by a black (seeFig 4). 9 AB Figure 4:
Three-periodic patterns in ECA composed of 20 cells. A: Repetitions of 2 black and 1 white cells. B: Repetitions of 2 white and 1 black cells.
Consider the periodic pattern formed by the repetition of one white and two black cells. There are32 rules which do not destroy preset periodic pattern of this type. These rules are:4, 5, 12, 13, 36, 37, 44, 45, 68, 69, 76, 77, 100, 101, 108, 109, 132, 133, 140, 141, 164, 165, 172, 173, 196,197, 204, 205, 228, 229, 236 and 237.All these numbers have binary expressions in the form of xxx x x where x is either 0 or 1. That is,the pattern consisting of a repetition of 2 blacks and 1 white is conserved under the rules keeping theconfigurations 2, 3 and 5 unaltered: (001) → (001), (010) → (010) and (100) → (100).Similarly, there are 32 rules which conserve a periodic pattern containing one black and two whitecells. The binary expressions of these rules have the form of x x xxx . That is, the pattern formed byrepetitions of 2 whites and 1 black is conserved under the rules keeping the configurations 4, 6 and 7unaltered: (011) → (011), (101) → (101) and (110) → (110). One can notice that there are four rules(76, 77, 204 and 205) which belong to both groups and allow the existence of both types of three-periodicpatterns.Our next task is to identify rules which allow recovery of three-periodic pattern perturbed by arandom and rare noise. Our simulations show that there are only two such rules: rule 133, which allowsrecovery of three-periodic patterns formed by two blacks and one white, and its complement rule 94,which allows recovery of the second type of three-periodic pattern. Note that the mirror image of rule133 is itself and its mirror complement is rule 94. Let’s consider the recovery of perturbed three-periodicpattern under rule 133 (whose binary representation is 10000101). Consider the fragment of periodicpattern . . . 001001001. . . and assume that the noise strikes it at either 2 nd or 3 rd or 4 th position. Theperturbed sequence becomes either . . . 011001001. . . or . . . 000001001. . . or . . . 001101001. . . respectively.Simple analysis shows that it takes 3, 2 or 3 time steps respectively to recover for each of the casesconsidered. Thus, the recovery under the rule 133 is local. All results obtained for rule 133 are naturallyextended to rule 94.We found no rules which allow the recovery of three-periodic patterns by means of propagating waves.Furthermore, there are no rules which allow the formation of threeperiodic patterns from random initialconditions.It is easy to show that there are three types of strictly four-periodic patterns which can form in ECA.Indeed, the four-periodic patterns in a generic form can be represented as ( a a a a ) n where elementsai can be either 0 or 1 and the chain has 4 n cells. We have 16 possible cases for four-periodicity, namely:0000, 0001, 0010, 0011, 0100, 0101, 0110, 0111, 1000, 1001, 1010, 1011, 1100, 1101, 1110 and 1111. Thefirst and last cases (0000 and 1111) are trivial and can be considered as one-periodic. Cases 0101 and1010 make the classical two-periodic pattern and we have studied them previously. All other cases, afterre-indexing, fall into one of the three cases: (0111) n , (0011) n or (0001) n . There are 16 rules which cankeep the preset 4-periodic pattern of each type. The four-periodic pattern of the type (0111) n is kept10y rules in the form 110 x xxx (where x is either 0 or 1); patterns of the type (0011) n - in the form x x x x , and patterns of the type (0001) n - in the form xxx x The study presented here was motivated by the problem of biological pattern formation governed bylocal (cell-to-cell) signals. As our primary interest was about segmentation in fly embryos, we havefocused on the formation of periodic patterns. Furthermore, since the segmentation takes place alongthe embryonic anterio-posterior axis and is essentially a one-dimensional process, we have modelled itusing a one-dimensional chain of logical elements represented by ECA. In the framework of this model,we have identified all rules allowing the existence or formation of periodic patterns. At the same time,our study was mainly focused on properties of ECA and, we believe, that it contributes to an extensiveresearch [15] performed by a large community of scientists in this field.In the presented study, we have explored whether formation of morphogenetic patterns in developingtissue can be modelled using ECA. The two-state model used in this study represents a chain of locallyinteracting cells, where cells in state “1” express some particular gene while in state “0” - don’t. Inter-actions of logical elements in ECA can be considered as representing contact (membrane-to-membrane)interactions between cells in biological tissues. The modelled interactions can be seen as regulating thedifferentiation of cells. For example, the two-periodic stationary pattern forming in the model representsa chain of cells where each second cell expresses a certain gene. We have identified sets of interaction ruleswhich allow existence, recovery from noise and formation from random initial conditions of various peri-odic patterns in CA. Particularly, we have found that 64 (out of 256) rules allow existence of two-periodicpatterns, 33 of them allow recovery of two-periodic patterns from rare and random noise and only sixrules (15, 85, 30, 86, 135 and 149) allow formation of two-periodic patterns from any initial conditions.Furthermore, our study of three- and more- periodic patterns has shown that although there are manyrules allowing their existence, only two rules allow the recovery of three-periodic patterns from randomand rare noise and none allow their recovery from noise or formation from random initial conditions.Here we have modelled the formation of periodic patterns in one-dimensional medium of constantsize with fixed boundary conditions. Although biological development takes place in three-dimensionalobjects the formation of simplest (i.e. periodic) patterns can be considered as one-dimensional process.Patterning in a growing tissue implies that the size of the modelled medium is changing over time andmodelling the formation of periodic patterns in growing medium is a task for future studies. Also,developing communities of bacteria often form circular chains to model which one should apply periodicboundary conditions. Most of the results obtained in the model with fixed boundary conditions remainunchanged for the model with periodic boundary conditions. Particularly, all results summarised inTable 1 remain correct. Rules listed in Table 2 allow the recovery of two-periodic patterns provided thatthe medium is composed of even number of cells. This excludes rules 30, 135, 86 and 149, for whichsingle perturbations transform into chaotic patterns. There are extra interesting results obtained for the11eriodic boundary conditions. For example, under rule 142, any perturbation circulates along the loopedchain infinitely many times (while in the case of fixed boundary conditions, cells in the right-hand sideborder do not recover). This result is valid for mirror image and mirror complement of rule 142, namely,for rule 212.It is known that the pattern associated with the expression of segmentation genes in the fly embryo isfour-periodic. As this pattern forms in a few steps and under the influence of maternal, gap and pair-rulegenes, one could consider this pattern as preset four-periodic which is allowed by a number of rules.However, it is also known that the periodic pattern formed by segmentation genes undergoes certaincorrection in position [16] and this wouldn’t be allowed by any of transition rules in two-state model.
A B CCell stateGene 1Gene 2
Figure 5:
Periodic patterns in a four-state model. Four- (A), five- (B) and six- (C) periodic patterns withillustration of how each state corresponds to different expressions of two genes. Cell state is indicated in blackif both genes are expressed; in dark grey - Gene 1 is expressed and Gene 2 is not; in light grey - Gene 1 is notexpressed and Gene 2 is expressed; white - both genes are not expressed.
Four-periodicity can easily be obtained if the model allows at least three distinct states for cells.While the three-state model doesn’t have direct biological implementation, the four-state model seemsto have better perspective as it can be viewed as modelling cells whose differentiation is associated withexpression of pair of genes. The fourperiodic pattern can correspond to the alternation of expressions oftwo genes as illustrated in Fig 5A. The four-state model is considerably more sophisticated comparedwith the two-state model and can be used to reproduce the formation four- and even more- periodicpatterns (Fig 5). While the number of possible transition rules in the two-state model is relatively small(256), this number is extremely large for the four-state model (4 , that is, every transition rule is acombination of 64 = 4 elementary rules, while each elementary rule has four forms). This numbermakes the analysis of 4-state models challenging. Our preliminary simulations show that the formationof periodic stationary patterns in four-state model is extremely sensitive to the initial conditions, i.e.these patterns form only when very special initial conditions are met. This may explain the multi- (four-)level of segmentation in the fly embryo. The model can account for interactions between segment polaritygenes with specific initial conditions set by the above three levels of patterning.Mathematical study of the segmentation in fly embryo performed in [6] and based on gene networkmodelling has indicated that all nuclei, during normal development, can be only in four distinct states.This points to the direction of further extension of the presented work: our extensive analysis of ECAcan be considered as a step towards similar analysis of four-state CA. Also, we note that more sophisti-cated CA (two-dimensional and including randomization of transition rules) have recently been used tosuccessfully model the dynamics of skin patterning in growing lizard [7]. Acknowledgements.
This work was funded by BBSRC grant BB/K002430/1. Authors are gratefulto Dr Igor Potapov for fruitful discussions. 12 eferences [1] Lewis Wolpert, Elliot Meyerowitz, Peter Lawrence, Rosa Beddington, Thomas Jessell, and JimSmith.
Principles of development . Oxford University Press, USA, 2 edition, 2002.[2] James D Murray.
Mathematical Biology II: Spatial Models and Biomedical Applications . 3rd edition,2003.[3] O. Vasieva, M. Rasolonjanahary, and B. Vasiev. Mathematical modelling in developmental biology.
Reproduction , 145(6):R175–84, 2013.[4] Sarah Hake and Fred Wilt.
Principles of Developmental Biology . W. W. Norton and Company, NewYork, NY, 2003.[5] D St Johnston and C Nsslein-Volhard. The origin of pattern and polarity in the drosophila embryo.
Cell , 68(2):201219, January 1992.[6] R´eka Albert and Hans G. Othmer. The topology of the regulatory interactions predicts the ex-pression pattern of the segment polarity genes in drosophila melanogaster.
Journal of TheoreticalBiology , 223(1):1–18, jul 2003.[7] Liana Manukyan, Sophie A Montandon, Anamarija Fofonjka, Stanislav Smirnov, and Michel CMilinkovitch. A living mesoscopic cellular automaton made of skin scales.
Nature , 544(7649):173179,April 2017.[8] J von Neumann. The general and logical theory of automata. 1948.[9] John G. Kemeny. Theory of self-reproducing automata. john von neumann. edited by arthur w.burks. university of illinois press, urbana, 1966. 408 pp., illus. $ 10.
Science , 157(3785):180–180,1967.[10] Gabriel Wainer, Qi Liu, Olivier Dalle, and Bernard P. Zeigler. Applying cellular automata and devsmethodologies to digital games: A survey.
Simulation & Gaming , 41(6):796–823, 2010.[11] Xin-She Yang and Y. Young. Cellular automata, pdes, and pattern formation. 2010.[12] Stephen Wolfram.
A new kind of science . Wolfram Media Inc., 2002.[13] Stephen Wolfram. Statistical mechanics of cellular automata.
Rev. Mod. Phys. , 55:601–644, Jul1983.[14] S. Wolfram.
Theory and Applications of Cellular Automata: Including Selected Papers, 1983-1986 .Advanced series on complex systems. World Scientific, 1986.[15] Janko Gravner and David Griffeath. Robust periodic solutions and evolution from seeds in one-dimensional edge cellular automata.
Theoretical Computer Science , 466:64 – 86, 2012.[16] Johannes Jaeger, Svetlana Surkova, Maxim Blagov, Hilde Janssens, David Kosman, Konstantin NKozlov, Manu, Ekaterina Myasnikova, Carlos E Vanario-Alonso, Maria Samsonova, David H Sharp,and John Reinitz. Dynamic control of positional information in the early drosophila embryo.
Nature ,430(6997):368371, July 2004. 13
Rule 30
In this appendix, we analyse the evolution of perturbed periodic pattern under the rule 30 and show thatthe pattern recovers from perturbation by driving the perturbed area to the right at a speed of at least2 cells (one strip) per 9 time-steps. As a starting point, we note that a single perturbation can result toformation of a large perturbed area and therefore, we will consider a sequence in the form (01) kxxx . . . (or1(01) kxxx . . . ), that is a sequence starting with k stripes followed by an arbitrary com-bination of 0s and1s. We check how many time steps it takes for any such combination to transform to (01) n +1 xxx . . . , thatis, to form an extra strip. We note that the transformation (01) nxxx . . . → (01) n + 1 x . . . is equivalentto the transformation 1 xxx . . . → x . . . and in the following will refer to the latter case (to reducenotations). Below, we split all sequences having the form 1 xxx . . . into (ten) groups so that the sequenceslisted in each following group take one more step for the required transformation. • st group : sequence is already in the form 101 x . . . , and thus the transformation took place inzero time steps. • nd group : sequences resulting in the formation of 101 xx . . . in one time step have the form1100 x . . . (i.e. 1100 x . . . → xx . . . in one time step). • rd group : sequences resulting in the formation of 1100 xx . . . in one time step (and thereforeresulting in the formation of 101 xxx . . . in two time steps) have the form 10000 x . . . . • th group : sequences resulting in the formation of 10000 xx . . . in one time step (and there-fore resulting in the formation of 101 xxx . . . in three time steps) have the forms 11111 xx . . . and111101 x . . . . • th group : sequences resulting in the formation of 11111 xxx . . . and 111101 xx . . . in one timestep (and therefore resulting in the formation of 101 xxx . . . in four time steps) have the forms10010 xxx . . . and 1001100 x . . . . • th group : sequences resulting in the formation of 10010 xxx . . . and 1001100 x . . . in one timestep (and therefore resulting in the formation of 101 xxx . . . in five time steps) have the forms11011 xxx . . . , 110101 xx . . . , x . . . , xx . . . and 1110011 xx . . . . • th group : sequences resulting in the formation of 11011 xxx . . . and 110101 xx . . . in one time step(and therefore resulting to the formation of 101 xxx . . . in six time steps) have the form 10001 xxx . . . ,while no any sequence results to formation of 11010000 x . . . , xx . . . and 1110011 xx . . . . • th group : sequences resulting in the formation of 10001 xxx . . . in one time step (and thereforeresulting in the formation of 101 xxx . . . in seven time steps) have the forms 11101 xxx . . . and111100 xx . . . . • th group : sequences resulting in the formation of 11101 xxx . . . and 111100 xx . . . in one timestep (and therefore resulting in the formation of 101 xxx . . . in eight time steps) have the forms1001101 xxx . . . and 100111 xxx . . . . • th group : sequences resulting in the formation of 1001101 xxx . . . and 100111 xxx . . . in one timestep (and therefore resulting in the formation of 101 xxx . . . in nine time steps) have the forms11010001 xxx . . . , xxx . . . and 111001 xxx . . . .14e note that all sequences having the form 1 xxx . . . have appeared in one of the above ten groupsand therefore the transformation never takes more than 9 steps. A summary of the time required for thetransformation 1 xxx . . . → x . . . for all possible sequences is given in the Table A.1.10000 xxx . . . . . . . . . x . . . xx . . . xxx . . . xxx . . . xxx . . . xx . . . xx . . . x . . . xx . . . x . . . x . . . xxxxx . . . xx . . . xx . . . xxx . . . xxx . . . Table A.1:
Summary of the number of time steps required for sequences in the form 1 xxx . . . to transform tothe sequence 101 x . . . for rule 30. Transforming sequences are indicated in the 1 st and 3 rd columns, while thenumber of time steps required for their transformation in the 2 nd and 4 th columns. Rule 15
In this appendix, we analyse the formation of two-periodic patterns under the rule 15 in a chain imposedto the random initial conditions. We note that under this rule, periodic patterns form first on theleft border and then expand to the right. Therefore, we start with the consideration of various initialcombinations on the left side of the chain and analyse how the periodic structure forms and expands tothe right. First, we consider initial conditions with the digit 1 on the left border of the chain. • Case A: xxx . . .
It is easy to see that none of these three digits will change over time. The firstdigit 1 is fixed (fixed boundary). The second digit, 0, stays unaltered following the configurationrule (101) → (101). The third digit, 1, also doesn’t change no matter what digit is on the forthposition following either (010) → (0
0) or (011) → (0 • Case B: xxx . . .
The second digit, 0, remains unaltered following the configuration: (100) → (1 → (0
0) or (001) → (0 xxx . . . ( case A ) in one time step: 100 xxx . . . → xxx . . . . • Case C: xxx . . .
The second digit, 1, turns into 0 following (110) → (1 → (1
0) or (101) → (1 xxx . . . which is case B . Thus, it takes two time steps to transfer from case C to(stable) case A: xxx . . . → xxx . . . → xxx . . . . • Case D: xxx . . .
The second digit, 1, turns into 0 following (111) → (1 → (1
0) or (111) → (1 xxx . . . (whichis case B ) in one time step and therefore obtain 101 xxx . . . (which is case A ) in two time steps:111 xxx . . . → xxx . . . → xxx . . . . • Cases E-H:
These cases deal with the initial conditions in the form 0 xxx . . . , i.e. with 0 on theleft border. It is easy to show that the first three digits in the sequence 010 xxx . . . do not changeover time, which is similar to case A considered above. The sequence 011 xxx . . . transfers into010 xxx . . . in one time step (similar to case B ), while the sequences 000 xxx . . . and 001 xxx . . . transfer into 0 xxx . . . in two time steps (similarly to cases C and D ).It is also easy to see that the pre-existing strips on the left side of the chain are conserved, that is,(10) n xxx . . . → (10) n xxx . . . and (01) n xxx . . . → (01) n xxx . . . . This observation, in the combina-tion with the above cases A-H , indicates that the periodic pattern forming on the left expands tothe right ((10) n xxx . . . → (10) n +1 xxx . . . or (01) n xxx . . . → (01) n +1 xxx . . . ) and every next stripappears in at least two time steps. Thus, the slowest way of formation of periodic structure impliesits expansion to the right at a speed of 1 cell per time-step. The formation of periodic pattern mayremain incomplete if its expansion to the right leads to occurrence of xxx
011 or xxx
100 where xxx is a periodic fragment (i.e. xxx = (10) n ) and therefore stable. This follows from the configurationrules (011) → (011) and (100) →→