Cellular Automata (CA) as a minimal model of pattern formation via local cell-cell signaling. (a) Example of a one-dimensional CA of length L = 15, two states and fixed boundary conditions. Given a pattern of cells in the blue (0) or white (1) state at time t, an update function f (the “update rule”) determines the pattern at time t + 1. This rule maps each of the 8 possible configurations of a cell with its two nearest neighbors to an output, the new state of the cell. Each rule is uniquely identified by the string of output values for all possible input configurations in a prescribed order. (b) By analogy to pattern-forming developmental systems, we define a genotype to phenotype to fitness relation for the CA models. All cells in the system use the same update rule, which specifies the dynamical behaviour of the system. The string of output values encoding the rule is thus interpreted as a genotype, whereas the pattern formation process produced by the rule is the phenotype. While the fitness could generally depend on various aspects of the pattern formation process, we assume here that fitness depends only on the final pattern: The fitness is defined as the fraction of the final pattern that matches the target pattern. In the shown example, the French Flag pattern is the target pattern, but other axial patterns are also considered.

Evolutionary search for high-fitness rules and multiple alignment. (a) Schematics of the evolutionary algorithm (EA). Starting with a population of rules, the algorithm uses mutation and selection steps to find rules with higher fitness. (b) Left: Phenotype of a high-fitness rule found with the EA. Center: Multiple alignment of the genotype of this rule with those of its high fitness neighbors (Φ > 0.9) up to a Hamming distance of 6, totalling 379 different rules (the first rule in the alignment is the EA output). The consensus rule is the result of the multiple alignment. Right: Phenotype of the consensus rule.

Dissecting the pattern formation mechanisms of the consensus rule. (a) Sequence logo for the consensus rule and its high-fitness neighbors (Φ > 0.9) up to a Hamming distance of 5 (in total 387 rules). The horizontal axis represents the 27 entries within the update rule (the input pattern is indicated for each entry to facilitate interpretation). The information content (decrease in Shannon entropy compared to a uniform distribution) in bits for each entry can be read off on the vertical axis as the total height of the stack of symbols. The relative height of the symbols represents their relative frequencies, with the most frequent output at the top. An entry of the update rule is considered as highly variable, if its entropy is more than 1 bit (information content less than log2(3) − 1 = 0.58). (b) Table of all single mutations of the consensus rule. Vertical axis shows to which state the output is mutated. The fitness of the mutant rule is color-coded (white indicates fitness above 0.9). Black squares mark the consensus output. (c) The condition that the target French flag pattern must be maintained once it is reached (steady-state) fixes seven entries of the update rule. (d) Partition of consensus rule entries into two patterning modules. The sorting module contains 12 entries and effectively implements a swap between white and blue states when white is to the left. The bulldozer module contains 15 entries and creates a red bulldozer state that moves to the right, while erasing any white and blue states to its right and alternately seeding blue and white states to its left. These two characteristics of the red bulldozer state are schematically represented by the grey boxes with the pointy arrow to white and blue states on the left, and the blunt-ended arrow to white and blue states on the right. (e-g) T-distributed Stochastic Neighbor Embedding (tSNE) plots of 3-state CA rules with fitness larger than 0.85, based on the genotypes S of the rules. Crosses represent the position of the consensus rule (red) and its symmetric partner (pink). The plots are colored based on rule fitness (e), overlap with conserved entries of the sorting module (f) or the bulldozer module (g). The left side of the plots in (f) and (g) shows high agreement with the sorting and bulldozer modules, while the right side shows high agreement with their respective symmetric partner rules (fig. S3c and d). Therefore, most high scoring rules combine the sorting and bulldozer module to achieve the French flag pattern.

Pure sorting strategy. (a) Pure sorting with only three states faces a dilemma at local configurations of three cells in an inverted French flag pattern. According to the sorting logic, both lateral states, blue and red, should move to the center, while the white state should move both to the right and left. (b) The dilemma is resolved with an extra state (colored black, numbered 3 in genotypes S) representing the superposition of blue and red. A resolution module is formed by the shown set of updates. (c) Example of how the resolution module operates. A random initial pattern with local inverted French flag configuration is sorted, maintaining the initial number of cells in each state. To undo the duplication of the white state in the first time step, the resolution module lets the black state chase the extra white state on its right and, when they meet, exchange both states for a blue and a red state, resolving the ambiguity of the black state and erasing the duplicated white state. The arrows link this example to the general scheme underlying the resolution module shown in (d). An ambiguous position is always converted into a black cell and a duplicated white state. If the state of the cell to the right of the duplicated white is blue, that blue state is swapped with the one inside the black state (case 1). Otherwise, the red state takes the place of the duplicated white, resolving both the ambiguity of the black state and the duplication of the white state (case 2).

Three different axial patterning strategies in 4-state rule space. In addition to the mixed patterning strategy of the 3-state consensus rule, the introduction of a fourth, transient state also enables the formation of a French flag pattern via a pure sorting strategy (left), or via a full erase-and-reconstruct patterning strategy (middle). These new patterning strategies are implemented by combining the three pairwise sorting modules for the French flag color states with the resolution module or the hidden bulldozer module, respectively. The resulting phenotypic patterning behavior (bottom) is shown for two different system sizes, starting from random initial conditions in 4-state space. For comparison, a 4-state version of the 3-state consensus rule implementing the mixed patterning strategy is also shown (right). See fig. S1 for the full definitions of the three shown rules.

Comparison of the patterning accuracy of the three patterning strategies of Fig. 5. Accuracy is measured via the mean fitness value (averaged over 1000 randomly chosen initial conditions). The maximal fitness of one is achieved only when a perfect French flag pattern is formed for all initial conditions. The mean fitness increases with the system size, L, for all patterning strategies, but it approaches the maximal fitness with different power laws (inset). Linear fits on a log-log scale (inset) yield the exponents −.509 ± .001 for the Mixed rule, −1.00 ± .01 for the Landscaping rule, and −.469 ± .001 for the Bubble rule. The mean fitness oscillates with a period of three, since perfect French flag patterns can be formed only when L is divisible by three.

Robustness of the three patterning strategies of Fig. 5. (a) Maintenance of the French Flag pattern under noisy updates according to the Mixed, Landscaping, and Bubble rule. (b) Estimated average position of the two boundaries within the tripartite pattern, as a function of the update error probability p, for each of the three rules. Each point was calculated as an average over 1000 iterations starting from the French flag pattern, and the boundary position was determined using a sliding window of three cells and a majority rule to associate each cell to the most prevalent state in its window (the boundaries are then set to the first position where cell states change from blue to white, and from white to red). (c) Examples of simultaneous growth and patterning under the Mixed, Landscaping, and Bubble rule, respectively, for two different growth rates (cell division probability per time step). (d) Fitness as a function of final system size for the three rules and additionally a modified version of the Landscaping rule (see main text). The vertical dashed lines mark the system sizes beyond which the average time of pattern formation (fig. S1) is expected to exceed the time required to reach a growth speed of one cell per time step. (e) Fitness as a function of the probability of asynchrony pA for the same four rules. (f) Examples of the patterning dynamics at a probability of asynchrony pA = 0.1 in a system of size L = 34.

Undifferentiated initial condition. Here, the black state is used as an undifferentiated state, in which all cells within the system are initialized (except for the fixed boundaries). (a) The Landscaping rule solves the French flag problem also for this initial condition. (b, c) An evolutionary search yields a different patterning strategy of two fronts starting from the boundaries and moving inward at speeds of 1 and 1/2 cells per time step, respectively, with the faster front emerging either form the left (b) or the right (c).

Tunability of the patterning process. (a-c) An evolutionary search for rules that form a French Flag pattern with a stripe width ratio of 1:2:1 from random initial conditions identified a rule (genotype 201 220 122 202 110 112 201 110 000), which (a) achieves this target ratio increasingly well for larger system sizes (dots show the average proportion of each state, dashed lines the aimed ratios, and solid lines the average value from L = 901 to L = 1000). (b) The patterning behavior of this rule is similar to that of the consensus rule, however (c) the bulldozer state (blue) produces white and red cells in a 2:1 ratio in its wake, as it travels from right to left (shown here using an initial condition with only two isolated blue cells in a white background). (d-e) Mutations of the Bubble rule can tune the resulting stripe widths in fine-grained steps, for instance to (d) increase the width of the blue stripe, in (e) 24 mutations that gradually increase the production of blue states from black states. (f-i) Generalization of the mixed and landscaping strategies to form axial patterns with more than three stripes, as illustrated here using either 5 or 10 states.

Results and analysis of evolutionary algorithm. For all simulations in this figure, the algorithm was started with 100 random rules and ran for 1000 iterations. Rule fitness was estimated by an average over 50 random initial conditions with CA size of L = 100. Each iteration, rules were mutated and, if less fit than the mutant, replaced by the mutant rule. Then, the 30 least fit rules were substituted by random ones. (cf. alg. S1) (a) Average and best fitness of the 100 rules in the list at each iteration as a function of evolutionary algorithm iterations. (b) Comparison between the ‘full’ EA using both mutations and replacement by random rules to modify the rule list, as well as modifications of the EA either only mutating the rules or only randomizing the least fit ones. (c) Distribution of maximum rule fitness found by the algorithm over 1619 runs (1620 for only mutation of rules, 1621 for only insertion of random rules) with the parameters specified before. (d) Phenotype of the fittest rule at the end of the run in (a), also fittest rule of the Full EA distribution in (c), for length L = 30, 100, and 300. (e) Dependence of the fitness of the rule in (a, b, d) on system length L for an average over 50 random initial conditions as in (a, b, c). Length L = 100 for which the fitness is usually calculated is marked with the red, dashed line.

Kymographs of consensus rule (a) and its symmetric partner (b) with genotype S = 000011012121011012111022002.

Consensus rule for different parameters. The consensus rule was calculated for 2 different 2 different neighborhood sizes and 3 different fitness thresholds. The differences to the consensus rule in the marked in red, with all consensus rules differed from the one selected by at most 5 mutations.

Properties of consensus rule and its modifications. (a) Number of red states in final pattern as a function of the number of reds in the initial pattern for consensus rule, consensus with mapping 9 changed to white and to blue. The number of final reds for the consensus rule is always slightly larger or equal than the number of initial reds, while for both modifications of the consensus the numbers are the same. (b) Same data as in (a) but averaged for each system length. The dashed black lines in the figure are the lines y = x and . (c) Example of how the bulldozer state in the consensus rule sometimes seeds more whites than blues and how modifying the output of mapping 15 to blue makes the seeding ratio 1 z 1. (d) Scaling of the fraction of blue, white, and red states with system length. Mean fitness is also plotted in black. (e) Dependence of the variance of the final difference in the number of whites and blues on the position of the leftmost red state. Each point is the variance over 1000 random initial conditions. For the consensus rule, this variance is capped because mapping 9 introduces a red state close to the left boundary even if there isn’t one in the initial state. Mutating this mapping makes the variance scale linearly with red state position, consistent with the widths being independent binomial random variables. (f) Scaling of the mean fitness with system size L for the consensus rule, its modification, and the fittest rule (rule 200 222 122 202 110 112 211 110 100). While the consensus fitness converges to 17/18, the fitness of the modified and the fittest rules converge to 1 as L → ∞. The analytical solution plotted is the equation . (g,h) Mean difference between the final and initial number of 0s (g) and 1s (h) as a function of the initial number of 0s (g) and 1s (h), respectively. 10000 samples of initial conditions with only blue and white states with a proportion given by the value on the horizontal axis and random sequence. Error bar corresponds to the standard deviation of the sample of each point

tSNE Plot of all rules with fitness larger 0.85 as described in fig. 3 and Fitness. (a) Number of rule entries that are the same as (b) Time to steady state, (c) Relevant entries (entries with information content larger log2(3) − 1 bits (cp. fig. 3 b) for a module that sorts red and white cells or (d) a blue bulldozer state moving to the left

Complete rule tables for the engineered rules in the k = 4 space.

Simulation examples of the prototypical rules for different system sizes L ∈ {9, 33, 99, 999}

Steady State Times for Consensus, Bubble, Landscaping, and Mixed rule. (a) Sampled distributions for L = 100 with the rough estimates of the average from the main text. (b) Scaling of average steady state time with length (dots), together with the rough estimates (dashed lines) and the fits (solid lines).

Average fitness of prototypical rules in growing systems. Average fitness as a function of the system size for the prototypical rules in 3 and 4 state space for the Consensus rule, Bubble rule, Landscaping rule and Mixed rule, as well as a modification of the consensus rule in entry 9 to 1 (from 2) and modification of the Landscaping rule in entry 36 to 1 (from 3). In each case, the fitness as a function of the initial and final length are shown for a system growing with a probability of 0.2% per time step per cell. Sample sizes are 10000 for L ≤ 10, 1000 for 10 < L ≤ 100 and 100 for L > 100. Further, the fitness for the corresponding non growing (static) system is plotted. Sample sizes in this case are 10000 for L ≤ 10, 1000 for L > 10. The vertical lines correspond to the length at which on average the exponential growth of the system exceeds the length at which the exponential growth of the system is faster then the pattern formation process defined by 0.002aLcrit = 1 with a being the average number of steps needed to complete pattern formation (cp. S1).

Simulation of a growing system for the Bubble (a, b) and Mixed (c, d) rule for initial system sizes L = 80 (a, c) and L = 160 (b, d). The system is growing with a division probability of 0.2% per time step per cell.

Simulation of a growing system for the Landscaping rule (a, b) and the Landscaping rule modified at position 36 (to 1) (c, d) for initial system sizes L = 80 (a, c) and L = 160 (b, d). The system is growing with a division probability of 0.2% per time step per cell. The modified entry suppresses the formation of new bulldozer states due to cell divisions

Simulation of a growing system for the consensus rule (a, b) and the consensus rule modified at position 9 (to 1) (c, d) for initial system sizes L = 80 (a, c) and L = 160 (b, d). The system is growing with a division probability of 0.2% per time step per cell. The modified entry suppresses the formation of new bulldozer states due to cell divisions

Rules for different target ratios (b, w, r) = (3, 1, 1): 200220122210110102110110100 and (b, w, r) = (2, 1, 1): 211220022210110102111110100. Left column: average proportion of each state (dots) together with the aimed ratios (dashed lines) as well as the asymptotic values (average value from L = 901 to L = 1000, solid lines). Middle column: exemplary kymograph from random initial conditions for L = 33. Right column: from two times adjacent blue and red cells (3, 1, 1), and randomly initial conditions with only blue and white cells (2, 1, 1), respectively.

Conservation of States for the Bubble Rule: Difference between number number of states in the blue (a), white (b), and red (c) color in final and initial state when starting from initial conditions without black states for L = 99 and 100 randomly drawn samples for every possible total number of blue, white and red initial conditions.

Modifications of the Bubble rule lead to different stripe width. Example (a, c) and average relative stripe width in final pattern as a function of a chosen mutation trajectory (b, d) for modification towards white (a, b) and red (c, d) state, respectively. In each case, the average stripe width can be adjusted using modifications

Distribution of different states after the 2L steps for L = 99 of the Consensus (a), Mixed (b) and Bubble (c-e) rule. Empirical distributions where sampled from 100000 samples, binomial distributions plotted with indicated probability p (see Analytical models). Although in some cases the binomial appears to be extremely close to the empirical distribution, the Kolmogorov-Smirnov-Tests fail in all cases (p << 0.001).