Emergent Patterns from Probabilistic Generalisations of Lateral Activation and Inhibition

The combination of laterally activating and inhibiting feedbacks is well known to spontaneously generate spatial organisation. 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 generalisations of lateral activation and inhibition. By doing so, we identify novel 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 realised 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 sign 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 organisation [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 oper-ating 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 new pattern-generating modules-new 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]; whereas the new feedback modules that we introduce generate a number of surprising static and dynamic patterning modes that are, to the best of our knowledge, new. 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 generalised in a different way.
The model runs on a 2D 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 (this is 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 (this is short-range activation).
A noise level T L is associated with the long-range inhibition: as T L increases, the longrange inhibition is wiped out. Similarly, a noise level T S is associated with the short-range activation. A parameter, β, 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 β = 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 (1) below).
Perturbing β away from zero, then rather colour configurations accumulate an excess of one particular colour-white is more abundant for β > 0, while black is more abundant for β < 0-and so β is called a symmetry breaking parameter, see Figure 1a.
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 8 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 configuration-dependent normalising factor, or equivalently the sum of numerators in (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 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 (see Figure 1b and 1c and Movies S1 and S2 and [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 (1) we have the following necessary and sufficient mean-field approximation for stationary attractors, For lateral inhibition only and no symmetry break (high noise on the short-range  (2)).

Methods
Each model was simulated on a 2D square grid of (l + 2r L  Simulations of the core model show that the local solutions predicted by the meanfield approximation (2) are invariably realised: for colour symmetry and lateral inhibition only (T S = ∞, T L 1, β = 0), attractors are approximately straight stripes, whereas curved labyrinths are generated by lateral activation and inhibition together (T S = T L 1, β = 0), see Figure 1b, Movies S1 and S2, and [6]. Now breaking the colour symmetry by increasing β > 0, then stripes bifurcate to white hexagonal nets while the predicted linear increase of w Lx with β 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 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 Supplementary Material. The wavelength of stripes, and of labyrinths when r S r L , is 4r L /3 see Supplementary Material and [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 wave number 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). Therefore, we searched for interesting new attractors by retaining the three rules (i)-(iii) while systematically generalising the model in different simple ways. Importantly, all generalisations in this study introduce new linearly independent terms into the mean-field approximation for stationary interfaces: we hypothesized that only such generalisations can produce attractors with interesting new patterns, see Supplementary Material. We consider each linearly independent term in the mean-field 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 β 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 simu- Sx associated with the short range, and symmetry preserving nonlinear Sx which compete with the short range activation. All other notation is left unchanged. The reaction kinetics are: where Z is the colour configuration-dependent normalising factor. We focus on pertur-

Gyrating labyrinths from a competing nonlinear activation on the long range
A second natural extension of the model, analogous to the extension of Section 3.1, includes in the exponent symmetry breaking nonlinear terms ±β L (w − b) 2 Lx associated with the long range, and symmetry preserving nonlinear activatory terms ±σ L (w − b) 3 Lx which compete with the long range inhibition. All other notation is left unchanged. The reaction kinetics are: where Z is the colour configuration-dependent normalising factor. Again, we focus on

Multi-colour hexagonal lattices and labyrinths from multi-colour lateral activation and inhibition
The 2 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 i Jx is the density of colour C i within range r J . We simulated 3 and 4 colours with activation at the interface and 1.
long-range inhibition only; 2. short-range activation only; 3. both long-range inhibition and short-range activation; and 4. symmetry preserving nonlinear terms that compete with the short-range activation (|σ S | > 1) or the long-range inhibition (|σ L | > 1). The dynamics are portrayed in Figure 4.
As depicted in Figure 4a, lateral inhibition only (T S = ∞, T L 1) generates stationary hexagonal lattices. Lateral activation only (T S 1, T L = ∞) 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 2 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 |σ S | > 1, multi-colour labyrinths bifurcate to patterning modes that are analogous to labyrinthine highways for 3 and 4 colours, see Figure 4b.
However, increasing the long-range symmetry preserving nonlinear competition |σ L | > 1 causes only 4-colour labyrinths to bifurcate to gyrating 4-colour labyrinths whereas 3colour labyrinths bifurcate to attractors that appear not to gyrate.

Travelling stripes and spirals and reorganising labyrinths from cyclic symmetry breaking
When the model is extended to more than 2 colours, a new mode of symmetry breaking is possible which, unlike in Sections 3.1 and 3.2, does not enforce an accumulation of one particular colour. As in Section 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: for 3 colours or 4 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 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 γ S (or γ L ) drives the cycle in the opposite direction on the short (or long) range. A similar mode of symmetry breaking has been called cyclic dominance and studied in a spatial and probabilistic context within [16].
For lateral inhibition only (T S = ∞, T L 1), cyclic symmetry breaking with |γ L | ≈ 1 causes stationary multi-colour lattices to bifurcate to travelling stripes (3 colours) or travelling lattices (4 colours), see Figure 5 and Movie S7. For lateral activation only (T S 1, T L = ∞), cyclic symmetry breaking with |γ S | ≈ 1 causes multi-stable uniform block attractors to bifurcate to cyclic spiralling attractors, see Figure 5 and Movie S8.
Strikingly, the focal points of the spirals appear to be stuck rigidly in one place and rarely wander within the domain. Simulations indicate that the number of spiral foci at the end of the simulation N varies linearly with l 2 /πr 2 S , see Figure 5 inset panel, top plot; N appears to be unaffected by the simulation's initial condition, see Figure 5 inset panel, middle and bottom plots.
For short-range activation and long-range inhibition together (T S ≈ T L 1), cyclic symmetry breaking on the short or long range (|γ S | ≈ 1 or |γ L | ≈ 1) causes multi-colour labyrinths to continually reorganise, see Movies S9 and S10.     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 new 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 recognised 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 long-range 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. [18], and the dynamic protrusions of filopodia, e.g. [19,20].
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 (β = 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 .

Acknowledgements
This work began under the supervision of Tom Duke at the London Centre for Nanotechnology and CoMPLEX, University College London. Tom Duke is deeply missed.
We thank David Wright for running a number of simulations and for stimulating discussions, Michael Cohen for introducing us to discrete models of lateral inhibition, and Pau Formosa-Jordan for a critique of the manuscript.

Funding
This work began while L.W. was supported by an EPSRC fellowship. , dynamics converge to stationary stripes; for short-range activation and long-range inhibition together (T L ≈ T S 1), stripes bifurcate to stationary labyrinths; weakening the long-range inhibition (T L T S ), eventually labyrinths bifurcate to bistable attractors. (c) For colour symmetry breaking (β = 0), stripes transit to hexagonal nets while labyrinths transit to irregular arrangements of pores; these transitions are represented by plots of f b , the summary statistic for 2-colour symmetry breaking transitions (upper plot is for lateral inhibition only; lower plot is for short-range activation and long-range inhibition together). The correspondence between patterns and f b is colour coded by red/orange/green. Simulation parameters and numbers of instances are listed in Table 1.  Table 1. (a) A nonlinear inhibitory lateral feedback that competes with the activatory feedback on the short range generates labyrinthine highways (left) and Kagome lattices (right). (b) Colour symmetry. Increasing |σ S | > 1 causes labyrinths to bifurcate to labyrinthine highways as represented by the sudden increase in f S (the black curve superposed on the blue curve is for l, r L , and r S rescaled by a factor of 1.5 while other parameters are held constant). The bifurcation occurs only when r S /r L < 1/2 (lower plot). (c) Colour symmetry break. In all plots σ S = 2. Increasing β causes labyrinthine highways to bifurcate smoothly to Kagome lattices for r S /r L = 1/2 (upper blue plot) or interwoven short and long-scale hexagonal nets for r S /r L = 1/3 (upper black plot). Increasing β S causes labyrinthine highways to bifurcate to train tracks (lower plots; blue (black) curve is for r S /r L = 1/2 (1/3)). In all simulations, T S = T L = 1/16; other parameters are in Table 1. (a) Colour symmetry. Labyrinths bifurcate to continually gyrating labyrinths as |σ L | 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 σ L = 2. As either β or β L increases, gyrating labyrinths bifurcate to uniform block colour attractors.