Emergent patterns from probabilistic generalizations of lateral activation and inhibition

The combination of laterally activating and inhibiting feedbacks is well known to spontaneously generate spatial organization. It was introduced by Gierer and Meinhardt as an extension of Turing's great insight that two reacting and diffusing chemicals can spontaneously drive spatial morphogenesis per se. In this study, we develop an accessible nonlinear and discrete probabilistic model to study simple generalizations of lateral activation and inhibition. By doing so, we identify a range of modes of morphogenesis beyond the familiar Turing-type modes; notably, beyond stripes, hexagonal nets, pores and labyrinths, we identify labyrinthine highways, Kagome lattices, gyrating labyrinths and multi-colour travelling waves and spirals. The results are discussed within the context of Turing's original motivating interest: the mechanisms which underpin the morphogenesis of living organisms.


Introduction
As a complex multicellular organism grows and develops, each one of its cells follows the same set of genomically encoded instructions, yet different cells beget drastically different fates so bringing about the organism's complex structure. Turing was among the first to present a powerful idea pertaining to this phenomenon, when he realized that a spatially homogeneous soup of just two chemically reacting species can spontaneously morph into a structured pattern owing to nothing more than the diffusion of these species, so long as the reaction kinetics are of the appropriate activatory or inhibitory nature and the diffusivities are sufficiently different [1]. Later Gierer and Meinhardt extended this notion to show how processes other than reaction-diffusion can potentially drive pattern formation: they demonstrated that any process which feeds back on itself over two lateral ranges-one short range that quickens or activates the process, the other long range that competitively slows or inhibits it-can spontaneously generate structural organization [2,3]. This combination of feedbacks has come to be known as short-range activation and long-range inhibition and is widely accepted as the key criteria for Turing-type patterning [3][4][5].
Groups of cells can effect lateral activation/inhibition by, for example, secreting ligands that diffuse on average a few cell lengths before degradation or binding to receptors; this binding triggers an intra-cellular signalling cascade which in turn increases/decreases the ligand's expression. Given the vast number of intra-and inter-cellular signals operating within and between cells, it seems probable that other types of lateral feedback beside lateral activation and inhibition also operate during development; furthermore, the structural diversity among living organisms is immense, whereas lateral activation and inhibition alone generate a limited range of patterns. What other types of lateral feedback may be operating to bring about the development of a multicellular organism? What patterns can be generated by recombinations of these feedbacks and what feedbacks are necessary and sufficient to generate a particular pattern?
Here, we study a model of interacting lattice sites that laterally activate and inhibit one another [6]. The lattice sites may be considered as a layer of immobile discrete cells, each executing a dynamic differentiation program to transit between Boolean states according to a common set of deterministic instructions but subject to some level of noise. The model has an important distinguishing feature: its discrete and probabilistic formulation makes it straightforward to introduce and to simulate additional pattern-generating modules-additional classes of lateral feedback operating alongside lateral activation and inhibition-in a systematic way: the discreteness allows us to write down simple equations governing the bifurcation diagrams, while the noise provides an in-built test of robustness against fluctuations and renders the long-term patterning dynamics independent of the model's initial state. For lateral activation and inhibition only, the model's final states of morphogenesis are either striped, hexagonally netted, labyrinthine, spotted or uniform block colour, patterns that are ubiquitous in nature and are well known to be generated by Turing-type models [6]; the additional feedback modules that we introduce generate a surprisingly extended range of static and dynamic patterning modes. The ubiquity in nature of the patterning modes generated by the original model-stripes, hexagons, labyrinths, spots-(see e.g. [7,8]; reaction-diffusion [1,[9][10][11]; directional solidification [12]; granular/fluid flows [13], hydrodynamic instabilities, animal furs, seashells and more [7,14,15]) suggest that concrete applications for these novel patterning modes may soon be discovered.

Model
The core model of discrete and probabilistic lateral activation and inhibition introduced in [6] is now described. In each section of the results, this core model is generalized in a different way.
The model runs on a two-dimensional square grid of lattice sites that are either black (B) or white (W ). From an initial configuration, lattice sites flip their colour, from black to white or white to black, one at a time. The rate at which a particular site flips its colour is determined by the following three rules: (i) a site can flip to a colour only if a neighbouring site has that colour (this is activation at the interface); (ii) the likelihood of flipping to a colour decreases with the density of sites of that colour within a particular long range r L (i.e. long-range inhibition), where r L )1 is measured in units of a lattice site diameter; (iii) the likelihood of flipping to a colour increases with the density of sites of that colour within a particular short range r S (i.e. short-range activation).
A noise level T L is associated with the long-range inhibition: as T L increases, the long-range inhibition is wiped out. Similarly, a noise level T S is associated with the shortrange activation. A parameter, b, determines the strength of the propensity for cells to flip to a particular colour, either black or white, independently of the lateral activation and inhibition. When b ¼ 0, this propensity is zero, then colour configurations converge to attractors that are 50% white and 50% black on average over multiple instances of the simulation; in this case, we say there is colour symmetry and the model's specification is unchanged when white is interchanged with black (see the reaction kinetics (2.1) below). Perturbing b away from zero, then rather colour configurations accumulate an excess of one particular colour-white is more abundant for b . 0, while black is more abundant for b , 0-and so b is called a symmetry breaking parameter (figure 1a). To interpret the rules (i)-(iii) in the context of cellular dynamics, each lattice site would represent a cell that can be in one of a number of discrete states. Its ability to change from a state A to a different state B would be dependent on contact with another cell of that particular state. The rate at which the transition A ! B happens would be regulated by two ligands produced by B cells, one activator of the transition that is degraded within an average length r L (long-range inhibition of type A), and one inhibitor of the transition that is degraded within a length r S (short-range activation of type A).
Precisely, from a given colour configuration, the probability that the next colour flip is at site x to colour C ¼ W or B is non-zero if and only if any of the eight neighbours of x has colour C (the set of neighbouring colours is denoted N x ); this probability has one of two possible forms depending on the colour C x of site x where ðw À bÞ Jx is the fraction of white minus black sites within range r J around x, and Z is the colour configurationdependent normalizing factor, or equivalently the sum of numerators in (2.1) over every non-zero probability colour flip. This completely specifies the model's reaction kinetics, except for initial and boundary conditions which are described in the next section. For the corresponding continuous time definition and a partial derivation of the mean-field equations (see the electronic supplementary material). When the model is run, invariably it converges to a macroscopically stationary attractor that is either patterned, or homogeneously noisy, or uniform block colour (figure 1b,c; electronic supplementary material, Movies S1 and S2, [6]). Patterned attractors are time-invariant or stationary if and only if sites along their interfaces flip from black to white and back again from white to black at equal rates, so from (2.1) we have the following necessary and sufficient mean-field approximation for stationary attractors, ðw À bÞ Lx =T L À ðw À bÞ Sx =T S % b for all sites x on interfaces: ð2:2Þ For lateral inhibition only and no symmetry break (high noise on the short range T S ¼ 1, low noise of the long range T L (1 and b ¼ 0), from (2.2) stationary attractors have w Lx % 1=2 along interfaces (because w Lx ¼ 1 À b Lx ), whereas for both lateral activation and inhibition with no symmetry break (T S ¼ T L (1 and b ¼ 0) then we have w Lx % w Sx along interfaces. Straight interfaces satisfy the condition w Lx % 1=2, whereas w Lx % w Sx can be satisfied by curved interfaces. Moreover, the mean-field approximation predicts a linear increase of w Lx with b for lateral inhibition only (T S ¼ 1 in (2.2)).

Methods
Each model was simulated on a two-dimensional square grid of (l þ 2r L ) Â (l þ 2r L ) sites. Initially, all sites of the simulation are equally and independently likely to be one colour among a list of permitted colours (black or white in the model above; the list is extended in the generalizations below). Thereafter, sites within the central square domain of l Â l lattice sites flip their colour according to the reaction kinetics specified in each section. Outside of this l Â l domain, colours remain fixed for all time-this constitutes the boundary condition.
Simulations of the core model show that the local solutions predicted by the mean-field approximation (2.2) are invariably realized: for colour symmetry and lateral inhibition only (T S ¼ 1, T L ( 1, b ¼ 0), attractors are approximately straight stripes, whereas curved labyrinths are generated by lateral activation and inhibition together figure 1b; electronic supplementary material, Movies S1 and S2, [6]). Now breaking the colour symmetry by increasing b . 0 then stripes bifurcate to white hexagonal nets while the predicted linear increase of w Lx with b holds true for the overall fraction of white sites in the l Â l domain, as shown in the plots in figure 1c which are fully described below; the relatively disordered labyrinths bifurcate to relatively disordered pores. The initial and boundary conditions appear to have no effect on the local structure of the final pattern so long as simulations are run until colour configurations have converged to an attractor [6]. Near to the boundary, stripes and labyrinths tend to align perpendicularly to the boundary and in this way the domain's geometry influences patterns' orientations (see the electronic supplementary material). The wavelength of stripes, and of labyrinths when r S ( r L is 4r L /3 (see the electronic supplementary material) [6]; increasing the domain size l Â l appears to have no effect on this wavelength. In all simulations, the domain size was set such that the wavenumber l=4r L =3 . 8.  We have seen, then, that in this model the combination of the three rules (i) activation at the interface, (ii) short-range activation, and (iii) long-range inhibition generates a range of Turing-type patterns-stripes, labyrinths, hexagonal nets and spots-and, moreover, that in several cases the correspondence between stationary patterning modes and parameter values, i.e. the bifurcation diagrams, can be roughly anticipated simply by inspecting the mean-field approximation (2.2). Therefore, we searched for interesting new attractors by retaining the three rules (i) -(iii) while systematically generalizing the model in different simple ways. Importantly, all generalizations in this study introduce new linearly independent terms into the mean-field approximation for stationary interfaces: we hypothesized that only such generalizations can produce attractors with interesting new patterns (see the electronic supplementary material). We consider each linearly independent term in the meanfield approximation to be a lateral feedback module. Each lateral feedback module has a corresponding free parameter.
In order to quantitatively represent transitions between patterning modes, a representative summary statistic was computed from each instance of the simulation and its variation with the corresponding parameter was plotted. White/black colour symmetry breaking transitions, where b is varied from zero, are quantitatively represented by the variation in f b which represents the fractional excess in white over black sites in the l Â l domain, as shown in the plots in figure 1c. In order to generate these plots, the simulation was run five times for each value of b with all other parameters held constant. The cross-hairs are f b computed for each instance of the simulation; the line-graph connects the averages of f b for different values of b. A similar format is followed for every plot in the article (see table 2). Table 1 lists the corresponding parameter values and the number of repeated instances of the simulation. Other summary statistics are introduced in Results; the corresponding and complete descriptions can be found in the electronic supplementary material. These summary statistics do not depend on the domain size or the end time of the simulation.

Results
In all sections, for a clear portrayal of the model's dynamics it is essential to view the corresponding movies. Parameters for each movie are listed in the electronic supplementary material.

Labyrinthine highways, train tracks and Kagome lattices from a competing nonlinear inhibition on the short range
A natural extension of the model, which introduces two new lateral feedback modules while retaining rules (i) -(iii), includes in the exponent symmetry breaking nonlinear terms +b S ðw À bÞ 2 Sx associated with the short range, and symmetry preserving nonlinear inhibitory terms +s 2 S ðw À bÞ 3 Sx which compete with the short-range activation. Increasing b S enhances the white colour flipping rate wherever there is an imbalance of colour on the short range; increasing js S j increases the strength of an inhibitory response that depends on the colour imbalance on the short range. Within a cellular tissue, b S = 0 could correspond to a scenario where a local imbalance of cellular states drives cells to a particular state, while s S = 0 would imply that a local imbalance promotes a negative feedback loop. All other notation is left unchanged. The reaction kinetics are 3 Þ Sx ÀT À1 L ðbÀwÞ Lx Þ Z : Table 1. Parameters for panels and plots in figures. The 'end time' is the number of transitions per simulation. n is the number of instances simulated for generating statistics. 't', 'm', 'b', 'p' stand for 'top row', 'middle row', 'bottom row', 'panel', respectively, in the corresponding figure or plot. 'V' stands for 'varying' in the corresponding figure or plot. The bifurcation point js S j ¼ 1 can be qualitatively explained as follows. So long as js S j , 1, the net short-range feedback is always activatory for all sites since jw À bj Sx is bounded above by 1-there is no possibility that the shortrange feedback switches to inhibitory. When js S j . 1, the net short-range feedback switches sign to become inhibitory for any configuration such that jw À bj Sx . 1=js S j, then qualitatively new dynamics are possible.
When colour symmetry is broken by perturbing b = 0, for r S /r L ¼ 1/2 labyrinthine highways transit to Kagome lattices-a pattern of intermeshed regular hexagons and triangles where the diameter of the hexagon is twice the side of the triangle-that are best known to feature in the atomic arrangement of the minerals Herbertsmithite and jarosite, see e.g. [16,17]. Whereas for r S /r L ¼ 1/3, labyrinthine highways transit to interwoven short and long-scale hexagonal nets (for r S /r L ¼ 1/3) as depicted in figure 2c; electronic supplementary material, Movies S4 and S5. When the shortrange symmetry breaking term is perturbed b S =0, labyrinthine highways transit to train tracks. As both b and b S are varied, the transitions between patterns are apparently gradual and smooth, as demonstrated by the plots of figure 2c which show gradual and smooth increases of f b (the blue curve is for r S /r L ¼ 1/2; the black curve is for r S /r L ¼ 1/3). breaking nonlinear terms +b L ðw À bÞ 2 Lx associated with the long range, and symmetry preserving nonlinear activatory terms +s 2 L ðw À bÞ 3 Lx which compete with the long-range inhibition. Analogous to the previous section, increasing b L enhances the black colour flipping rate wherever there is an imbalance of colour on the long range; increasing js L j increases the strength of an activatory response that depends on the colour imbalance on the long range. All other notation is left unchanged. The reaction kinetics are

Gyrating labyrinths from a competing nonlinear activation on the long range
S ðw À bÞ Sx ÀT À1 L ððwÀbÞþb L ðw À bÞ 2 À s 2 L ðw À bÞ 3 Þ Lx Þ Z and B[N x : C x ¼W !B withprob: where Z is the colour configuration-dependent normalizing factor. Again, we focus on perturbations to labyrinths: in all simulations T S ¼ T L ¼ 1=16 ( 1. Electronic supplementary material, Movie S6 and figure 3 portray the dynamics. For colour symmetry (b, b L ¼ 0), as depicted in figure 3a and electronic supplementary material, Movie S6, stationary labyrinths bifurcate to continually gyrating labyrinths as js L j increases beyond approximately 1.5 (r S =r L ¼ 1=2) or 2.0 (r S =r L ¼ 1=3). This is quantitatively represented in the plots of figure 3a by the summary statistic f L which represents the time-averaged speed of movement of the interface (see electronic supplementary material for details). f L appears to be unaltered when all length-scales are simultaneously rescaled while other parameters are held constant (compare the blue and green curves in the plot of figure 3a), indicating that the structure of the bifurcation diagram is unaffected by discretization effects of the lattice. The dependency of the bifurcation point on r S /r L may be due to more varied interface geometries that are got by increasing r S /r L , since some variations may be more likely to creep into unstable configurations that must then gyrate to new configurations that are locally stable. For s L ¼ 2, the value of r S /r L must exceed a threshold of approximately 1/3 in order for this bifurcation to occur (not shown). An argument precisely analogous to that in the previous section explains why the bifurcation point for js L j must be greater than 1. In the electronic supplementary material, we demonstrate that gyrating labyrinths persist upon perturbations to the boundary condition.
When colour symmetry is broken by perturbing b = 0, gyrating labyrinths first freeze to be stationary labyrinths. They then bifurcate discontinuously to uniform block colour attractors once the magnitude of b exceeds a threshold as depicted in the upper plot of figure 3b (s L ¼ 2 and r S /r L ¼ 1/2). The bifurcation is similar when perturbing b L = 0 but appears to be continuous.

Multi-colour hexagonal lattices and labyrinths from multi-colour lateral activation and inhibition
The two-colour models in the previous sections are extended to an arbitrary number of colours denoted by C i , i ¼ 1, 2, . . . , n. All symmetry breaking terms are omitted in this section. The reaction kinetics can then be described by a single expression: gives short-range activation and long-range inhibition, and c iJ x is the density of colour C i within range r J . We simulated three and four colours with activation at the interface and (i) long-range inhibition only; (ii) short-range activation only; (iii) both long-range inhibition and short-range activation; and (iv) symmetry preserving nonlinear terms that compete with the short-range activation (js S j . 1) or the long-range inhibition (js L j . 1). The dynamics are portrayed in figure 4.  (a) Colour symmetry. Labyrinths bifurcate to continually gyrating labyrinths as js L j increases beyond a threshold greater than 1 as represented by the increase in f L (black plot is for r S /r L ¼ 1/3; blue and green plots are for r S /r L ¼ 1/2; in the green plot, l, r L and r S have been rescaled by a factor of 1.5 compared with the blue plot). (b) Colour symmetry break. In all plots, s L ¼ 2. As either b or b L increases, gyrating labyrinths bifurcate to uniform block colour attractors. (Online version in colour.) rsif.royalsocietypublishing.org J. R. Soc. Interface 13: 20151077 As depicted in figure 4a, lateral inhibition only (T S ¼ 1, T L ( 1) generates stationary hexagonal lattices. Lateral activation only (T S ( 1, T L ¼ 1) generates multi-stable block colour attractors where finally each site has the same colour and all of the permitted colours are equally likely. Short-range activation and long-range inhibition together (T S % T L ( 1) generate multi-colour labyrinths. These attractors are analogous to two colour stripes, bistable uniform blocks and labyrinths depicted in figure 1.
When multi-colour labyrinths are perturbed by increasing the short-range symmetry preserving nonlinear competition js S j . 1, multi-colour labyrinths bifurcate to patterning modes that are analogous to labyrinthine highways for three and four colours (figure 4b). However, increasing the long-range symmetry preserving nonlinear competition js L j . 1 causes only four-colour labyrinths to bifurcate to gyrating four-colour labyrinths, whereas three-colour labyrinths bifurcate to attractors that appear not to gyrate.

Travelling stripes and spirals and reorganizing labyrinths from cyclic symmetry breaking
When the model is extended to more than two colours, a new mode of symmetry breaking is possible which, unlike in §3.1 and 3.2, does not enforce an accumulation of one particular colour. As in §3.3, because there are no symmetry breaking parameters that are associated with particular colours, the reaction kinetics can be described by a single expression: where again enforces short-range activation and long-range inhibition, and where ðMc Jx Þ j is the jth element of the vector for three or four colours, respectively, and M can be defined similarly for n colours. The circulancy of M drives the dynamics in a symmetry breaking cyclic colour ordering C 1 ! C 2 ! Á Á Á ! C n ! C 1 for g S , g L . 0. Colours can no longer be arbitrarily exchanged for one another without changing the model's specification, yet there is no a priori propensity for one particular colour to accumulate. Reversing the sign of g S (or g L ) drives the cycle in the opposite direction on the short (or long) range. This circulancy condition could correspond to a scenario where cells transit between states as partly governed by a non-transitive rock-paper-scissors-like dynamic; for example, high local density of state C 1 enhances the likelihood of state C 2 possibly via the diffusivity of some C 1 -produced species, similarly high local density of state C 2 enhances the likelihood of state C 3 and high local density of state C 3 enhances the likelihood of state C 1 . A similar mode of symmetry breaking has been called cyclic dominance and studied in a spatial and probabilistic context within [18].  For short-range activation and long-range inhibition together (T S % T L ( 1), cyclic symmetry breaking on the short or long range (jg S j % 1 or jg L j % 1) causes multi-colour labyrinths to continually reorganize (see electronic supplementary material, Movies S9 and S10). These multi-colour reorganizing labyrinths appear sensitive to perturbations to the boundary condition, unlike the cyclic spiralling attractors (see the electronic supplementary material).

Discussion
The models we have presented are simple generalizations of the well-known pattern-producing dynamics of lateral activation and inhibition. Yet, despite their simplicity, to our surprise they have produced a broad range of dynamics and attractors including labyrinthine highways, Kagome lattices, gyrating labyrinths and corresponding multi-species analogies. In some cases, to the best of our knowledge, the models constitute novel phenomenological mechanisms for generating such patterning modes. We anticipate that these attractors are robust to changes in details of the model such as the square geometry of the lattice or the precise formulation of the reaction kinetics and rather that they depend on the class of interactions implemented by the model. Whether the geometry of the square lattice impacts upon results could be better established by adapting our code, which is available on request and uses plugins from the software Processing [19], to run on hexagonal or Voronoi grids or by completing the derivation of the mean-field equations (see the electronic supplementary material) and simulating the resulting system of PDEs. We further anticipate that for each class of interactions all qualitatively distinct attractors were identified: a simple mean-field equation governing the model's output allowed us to search for new patterns systematically. The ease of implementation of simulations, the ease of interpretation of bifurcation diagrams, and the capacity for systematically exploring parameter space are strengths of this approach compared with continuum deterministic reactiondiffusion systems, which can be coupled to produce patterning modes similar to labyrinthine highways, gyrating labyrinths and cycle swirls as elegantly shown in [20,21] (table 2). Since we believe that the model outputs do not depend on details of the implementation, we explore and speculate about the likely applicability of the models to biology. Firstly, our study indicates that such complex patterns-beyond stripes, labyrinths and hexagonal nets-such as labyrinthine highways or Kagome lattices may be generated by a collection of cells with the following properties: (i) a master regulator that once expressed initiates the expression of two diffusing ligands and a membrane-junction signal; (ii) receiving the membrane-junction signal from a neighbouring cell is an absolute requirement for the expression of the master regulator; (iii) the first ligand diffuses in excess of a few cell lengths on average before being degraded; or it binds to cell-surface receptors to trigger a signal that decreases the likelihood of expression of the master regulator; (iv) the second diffusing ligand, which tends to be degraded within a range that is short compared with the range of the first diffusing ligand, upon binding to a receptor either increases or decreases the likelihood of expression of the master regulator depending on whether the number of bound receptors is below or above a threshold. However, it must be recognized that our model is only an abstract representation of (i) -(iv); in particular, it represents only cases where the time-scales of diffusion, degradation and membrane-junction signalling are much faster than the time-scale of the master regulator's response to such signals. There are other possible longrange signalling mechanisms besides diffusion, operating over lengths which span multiple cells, that can account for lateral feedbacks: in animals, these include the active migration of cells, e.g. [22], and the dynamic protrusions of filopodia, e.g. [23,24].
Secondly, our toy model motivates us to speculate about a possible dynamical feature of some developmental programs. In the model, when parameters enforce what we have called colour symmetry (b ¼ 0)-a prerequisite for e.g. stripes and labyrinths-that region of parameter space is comparatively rich in the diversity of final patterns that are generated; consequently, the final pattern responds sensitively and flexibly to sustained changes in parameter values. The broad range of natural patterns suggests that this might be evolutionarily advantageous. While genetic drift would tend to take the system away from colour symmetry, we expect evolution to act as a tuning force which maintains the sensitivity and flexibility of developmental programs. 1 Data accessibility. The datasets and mathematical derivations supporting this article have been uploaded as part of the electronic supplementary material and all movies can be found at https://www. repository.cam.ac.uk/handle/1810/253075.