*Perspective* **Size-Regulated Symmetry Breaking in Reaction-Diffusion Models of Developmental Transitions**

#### **Jake Cornwall Scoones 1,†, Deb Sankar Banerjee 2,† and Shiladitya Banerjee 2,\***


Received: 21 April 2020; Accepted: 6 July 2020; Published: 9 July 2020

**Abstract:** The development of multicellular organisms proceeds through a series of morphogenetic and cell-state transitions, transforming homogeneous zygotes into complex adults by a process of self-organisation. Many of these transitions are achieved by spontaneous symmetry breaking mechanisms, allowing cells and tissues to acquire pattern and polarity by virtue of local interactions without an upstream supply of information. The combined work of theory and experiment has elucidated how these systems break symmetry during developmental transitions. Given that such transitions are multiple and their temporal ordering is crucial, an equally important question is how these developmental transitions are coordinated in time. Using a minimal mass-conserved substrate-depletion model for symmetry breaking as our case study, we elucidate mechanisms by which cells and tissues can couple reaction–diffusion-driven symmetry breaking to the timing of developmental transitions, arguing that the dependence of patterning mode on system size may be a generic principle by which developing organisms measure time. By analysing different regimes of our model, simulated on growing domains, we elaborate three distinct behaviours, allowing for clock-, timer- or switch-like dynamics. Relating these behaviours to experimentally documented case studies of developmental timing, we provide a minimal conceptual framework to interrogate how developing organisms coordinate developmental transitions.

**Keywords:** symmetry breaking; pattern formation; reaction-diffusion; developmental transitions

#### **1. Introduction**

In developmental systems, it is important for the mechanistic constituents to "know" about the size of the living system as a whole [1,2]. This is most apparent in developmental transitions, which in many cases only proceed when cells or tissues have reached a critical size. Such control strategies allow living systems to couple developmental time to their size and geometry. What is the physical basis of these phenomena?

One can envisage two broad classes of size-control mechanisms: The first proposes size-dependent transitions are under external regulation: in developing tissues, this could be manifested as a cell-intrinsic clock, whereby a transition is achieved after a pre-defined time interval [3], or a gradient-based mechanism, whereby a distance critically far from the source of signalling molecules triggers a transition [4]. Alternatively, size-regulation could be an emergent property of collective decision-making, akin to quorum sensing [5–7]: communication between mechanistic constituents allows the system to sense its size. A priori, these control mechanisms have several conceptual benefits: synchrony in a transition is more robust, given decisions are made collectively, and such decision-making does not rely on a subset of constituents (e.g., source cells in gradient generation), instead being decentralised [1,8]. Can we find examples of emergent size-regulation from collective decision-making in living systems?

Many developmental transitions couple size- and temporal-control to a change in polarity regime: at a critical size, the system may spontaneously break symmetry to polarise, depolarise, bipolarise, or even radically change its pattern. We hypothesise that size regulation of such developmental transitions are an emergent property of the many mechanisms of polarisation. This view provides a framework for understanding developmental time [9], placing a critical emphasis on system size.

Here, we outline a minimal reaction–diffusion model for size-dependent polarisation in developing systems, arguing that the underlying regulatory motifs can be understood via a substrate-depletion feedback motif coupled with the growth of the system. We then overview cases of size-dependent symmetry breaking across scales, focusing on biochemical systems. We consider size-dependent decision-making within individual cells through to analogous processes in developing multicellular systems, proposing that our minimal model can help unify these divergent processes within a common theoretical framework. We conclude by speculating on the role of these mechanisms in coupling size-dependent transitions to developmental time.

#### **2. Reaction–Diffusion as a Framework to Understand Size-Regulated Symmetry-Breaking**

Pattern forming reaction–diffusion (RD) systems [10] are widely used to characterise the complex networks of molecular and cellular interactions that underlie biological symmetry breaking. In these systems, pattern formation can arise spontaneously driven by feedback motifs between diffusible molecules (intracellular polarity proteins or extracellular morphogens). Since Turing's insight in 1952, many different RD motifs have been proposed as the physical basis for the emergence of developmental patterns across scales [11]—from polarity establishment at the scale of a single cell [12] to pattern formation on the scale of a whole organism [11,13].

Perhaps the most famous phrasing of an RD system is the activator–inhibitor circuit, originally developed by Gierer and Meinhardt [14]. In activator–inhibitor systems, an activator molecule promotes its own production as well as the production of its fast-diffusing inhibitor that suppresses autocatalytic production of the activator. This motif has remarkable explanatory power across contexts, being used to describe the spontaneous establishment of hair follicle spacing [15], left-right asymmetry establishment in vertebrates [16], skeletal patterns in growing limbs [17–19], as well as pole-to-pole oscillation of Min proteins during bacterial cell division [20], and self-organisation of Rho GTPases in the animal cell cortex [21,22]. Substrate depletion models can also yield the spontaneous emergence of periodic patterns [23,24]. In these models, the activator consumes its own substrate to promote its autocatalytic production, leading to out-of-phase patterning of the activator and the substrate molecules (Figure 1a). For example, a substrate-depletion model was been used to explain lung branching, explaining out of phase patterns of gene expression between Shh (the activator) and FGF (the substrate) [25].

Substrate-depletion models are particularly relevant in the study of intracellular pattern formation and polarity establishment as feedback can be phrased in a mass-conserved manner. In such models, patterns emerge by the redistribution of polarity proteins [26,27]: proteins that form a polarity patch engage in self-recruitment, acting as activators, but this positive feedback is limited by a finite pool of (typically cytoplasmic) subunits. Variants on mass-conserved substrate depletion models have been used to understand PAR polarity establishment in the *C. elegans* embryo [28–30] and Cdc42 polarisation in *S. cerevisiae* [31–33]. Further, such models exhibit dynamic regimes, for example helping to explain oscillations in the *E. coli* Min-protein system [34–36].

While the mechanistic constituents and precise feedback architectures of RD mechanisms differ, many rely on a central motif of local activation and long-range inhibition [37]. This concept has acted as an important heuristic in framing models of biological pattern formation, but importantly unifies diverse RD systems within a common mathematical framework. Specifically, recent theoretical

models have demonstrated that most RD models of pattern formation can be approximated by the same mathematical formulation: the Swift–Hohenberg equation [38]. Strikingly, other mechanisms of periodic pattern formation that rely on cell movement [39] or mechanical instabilities [40] also rely on local-activation and long-range inhibition [38,41]. Therefore, many dynamical features of these models are applicable across systems and length scales.

In this perspective, we restrict our focus to RD models for biochemical pattern formation, elucidating the biological significance of a common dynamical feature shared across many motifs: the role of a critical system size for symmetry breaking [23,42,43]. To demonstrate this, we develop a mathematical framework for a mass-conserving RD system with a feedback motif similar to activator–substrate models for cell polarisation (Figure 1a).

**Figure 1.** Size-regulated symmetry breaking in activator–substrate model. (**a**) Pattern formation in a model of positive feedback coupled to a finite constituent pool. (**b**) Patterns form above a critical system size (*L*∗), corresponding to the largest mode where the homogeneous state becomes unstable and the system breaks symmetry. The order parameter O for symmetry breaking is defined as O(*L*) = *L* <sup>0</sup> |*dS*/*dx*|/ max *<sup>L</sup> L* <sup>0</sup> |*dS*/*dx*| , where O is zero for a homogeneous state and O = 1 for a symmetry broken patterned state. All parameters other than system size was kept constant in this analysis and initial perturbations were so chosen that *N*/*L* is constant, where *N* = *P*(*x*, *t*) + *S*(*x*, *t*) *dx* is the total pool size. Parameters: *DP* = 1, *DS* = 0.05, *k*on/*k*off = 20, *S*<sup>2</sup> <sup>0</sup> = 10 and *N*/*L* = 1.5. (**c**) Phase diagram in the plane of autocatalytic activity *κ* and the Hill saturation parameter *S*<sup>2</sup> <sup>0</sup>, showing three different phases: homogeneous state (black), symmetry broken saturated state (green), and symmetry-broken unsaturated state (red). Colormap (green to red) denotes the average value of the reaction rate *Fav* = *F*(*x*)*dx*, computed in the high density region, with *Fav* = 0 in the saturated state and *Fav* = 0 in the unsaturated state. *Fav* = 0.01 (blue points) defines the crossover value from the saturated state to the unsaturated state. Parameter values are the same as (**b**) except for *DS* = 0.01 and *L* = 3. Inset: Profiles of *S*(*x*). (**d**–**e**) Order parameter O for the unsaturated (**d**) and the saturated (**e**) regimes, showing that the symmetry of the homogeneous state is broken beyond a critical system size. We used periodic boundary conditions for all numerical simulations, unless otherwise specified. For initial conditions, we assumed a homogeneous *P*(*x*) and a sinusoidal *S*(*x*) profile of large wavelength.

#### **3. A Minimal Model for Size-Regulated Symmetry-Breaking**

To analyse the role of system size on the timing of symmetry breaking in a biochemical system, we consider a minimal model for a mass-conserved substrate-depletion system. Specifically, we model the spatiotemporal dynamics of a regulatory structure *S* in a living system of size *L* and coupled to a finite pool of building blocks. Let *P*(**x**, *t*) denote the concentration of building blocks in the subunit pool at location **x** at time *t*, and *S*(**x**, *t*) is the concentration of building blocks incorporated in the regulatory structure. *S* increases in amount by depleting the subunit pool *P*, and *S* can undergo dissociation into *P* (Figure 1a). The coupled dynamics of *S* and *P* are given by

$$
\partial\_t \mathcal{S} = D\_s \nabla^2 \mathcal{S} + k\_{\text{on}} P f(\mathcal{S}) - k\_{\text{off}} \mathcal{S} \,, \tag{1}
$$

$$
\partial\_t P = D\_p \nabla^2 P - k\_{\text{on}} P f(S) + k\_{\text{off}} S \, , \tag{2}
$$

where *Dp* and *Ds* are the diffusion constants of the subunits and the structure *S* (*Ds Dp*), *k*on parameterises the association rate of subunits to *S*, and *k*off is the constant rate of disassembly of *S*. The function *f* defines the size-dependence of the autocatalytic production rate of *S*. We assume the functional form *f*(*S*) = *Sn*/(*S<sup>n</sup>* <sup>0</sup> + *<sup>S</sup>n*), where the constant *<sup>n</sup>* > 1 controls the strength of cooperative assembly of *S*. The total amount of *P* and *S* remains conserved at all times, i.e., *N* = *P*(*x*, *t*) + *S*(*x*, *t*) *dx* = constant, and is assumed to scale linearly with the system size *L*.

As the structure *S* grows by locally depleting the pool *P*, localized patterns of *S* will exist in low density regions of *P* as long as *DP DS* (Figure 1a). The symmetry of the homogeneous state breaks and patterns appear above a critical domain size, making this transition size-dependent (Figure 1b). The critical size can be obtained from linear stability analysis as

$$L^\* = \left(\frac{2D\_S D\_P}{D\_S \partial\_P F + D\_P \partial\_S F}\right)^{\frac{1}{2}}\tag{3}$$

where *F* = *k*on*P f*(*S*) − *k*off*S*, and the derivatives are evaluated at the homogeneous steady state. This property of size dependent symmetry-breaking can serve as a decision-making rule to enact developmental state transitions when the system size reaches a critical value *L*∗. As the system size gets larger more discrete structures will emerge (Figure 2a). Such sequential pattern formation may regulate size-dependent developmental transitions, as we discuss later.

The *S*-*P* model introduced above (for *n* = 2) is similar to previously studied mass-conserved RD models such as wave-pinning [26] and the Turing-like autocatalytic model [31]. These activator– substrate models exhibit two distinct dynamic regimes [44]: the wave-pinning regime, characterised by wide mesa-like patterns and saturated subunit association kinetics and the Turing regime that yields narrow concentration peaks by virtue of competition between structures. The Turing regime operates below saturation (*S S*0), where a winner-take-all competition between structures asymptotically results in a single concentration peak [31,45,46]. By contrast, in the regime above saturation (*S S*0), the structures can co-exist for very long timescales, with the timescale of coarsening determined by parameters such as the diffusion coefficients or the reaction fluxes [47]. The *S*-*P* model exhibits both the saturated and the unsaturated regimes that can be obtained by tuning the strength of autocatalytic activity (*κ* = *k*on/*k*off) and the Hill saturation parameter *S*<sup>2</sup> <sup>0</sup> (Figure 1c). Both these dynamical regimes exhibit size-dependent symmetry breaking (Figure 1d,e), as well as sequential pattern formation for increasing system sizes (Figure 2). In the latter case, we assume that the subunit density is constant for increasing *L*, *L* <sup>0</sup> (*S* + *P*) *dx* ∝ *L*, as macromolecular contents often scale with cell size [48]. For appropriate choice of model parameters in the unsaturated regime, an increase in total subunit pool size coupled with local depletion of the subunit pool gives rise to coexisting peaks of the same height over biologically relevant timescales, much shorter than the long timescale of structure coarsening.

**Figure 2.** Sequential pattern formation in growing domains. (**a**) Kymographs of *S*(*x*, *t*) for increasing system size. Figures show spontaneous pattern formation via symmetry breaking of the homogeneous state above a critical system size *L*∗. As we simulate larger systems (Equations (1) and (2)), multiple patterns emerge in a size-dependent manner. Simulations were done as described in Figure 1b. (**b**) Number of patterns as a function of system size, for both the saturated and unsaturated regimes of the *S*-*P* model. Parameter values for (**a**,**b**) are the same as Figure 1b,c, respectively. Small amplitude random (uniform distribution) initial conditions for *S* and *P* were used for (**b**).

#### **4. Critical Size for Polarisation Can Be Utilised to Enact Cell State Transitions**

Our minimal model demonstrates how a positive-feedback motif coupled to features of system size, such as a limiting cytoplasmic pool, can yield size-regulated symmetry breaking of regulatory structures. In this section we explore biological realisations of our model around the bifurcation from unpolarised to polarised, arguing that cells may utilise these size-dependent properties to coordinate state transitions.

#### *4.1. Cell Size Dependent Transition from Asymmetric to Symmetric Division in the Early C. elegans Embryo*

The polarisation of the early *C. elegans* embryo has become a paradigm in biological symmetry breaking [49]. Anterior-posterior (AP) polarity in *C. elegans* is established before the first cell division [50]. Polarity establishment is achieved by the segregation of two groups of partitioningdefective (PAR) proteins to the anterior versus posterior [28–30]. Initially anterior PARs (aPARs) cover the entire membrane of the egg, but upon fertilisation at the posterior, serving as the symmetry breaking cue [51], aPARs segregate anteriorly, and posterior PARs (pPARs) posteriorly. Segregated PARs coordinate polarised division, whereby the division plane is set by the boundary of the two PAR domains [52,53].

PARs bind the membrane from a common, finite cytoplasmic pool and diffuse freely. Symmetry breaking is achieved by phosphorylation-dependent mutual inhibition between aPARs and pPARs [54,55], patterning the cell membrane into two polarisation domains. This double-negative feedback structure reinforces biases in the localisation of PARs and thus plays a similar role as the positive feedback motif in our minimal model. Indeed explicit substrate-depletion models yield phenomenologically identical results [56].

The boundary between PAR domains is regulated by the relative diffusivities of aPARs and pPARs, as well as the relative off rates. Boundary length is set by *LD* <sup>=</sup> <sup>√</sup>*D*/*k*off, where *<sup>D</sup>* is the aPAR diffusion constant and *k*off is the aPAR dissociation rate from the membrane. Therefore, provided diffusion and dissociation rates are independent of system size, *LD* will also be independent of cell size. Modelling confirms that this holds true regardless of structural differences in the model [56]. This length-scale thus sets a minimum cell size that can sustain polarised PAR domains (grey region in Figure 3a): below this critical size, diffusion overwhelms the capacity for PAR segregation, resulting in homogeneous localisation of either aPARs or pPARs (pink and blue regions of Figure 3a)

The physical critical size limit may be used by developing *C. elegans* embryos to coordinate a developmental transition. Sequential divisions of generating the first three posterior cells (P1-3) are asymmetric, each generating two daughters of different fates via the segregation of germline determinants. This pattern shifts to a symmetric mode by the third division of P4 (Figure 3a), generating the two founding cells of the germline lineage (Z2/Z3) [57]. Divisions are fast, meaning cell volume declines progressively, falling beneath the theoretical critical cell size for polarisation by P4 [56]. By quantifying division symmetry using 3D reconstructions of PAR distributions, Habatsch et al. [56] found that the polarisation regime shifts from asymmetric to symmetric by P4. The timing of this regime shift can be changed by reducing embryo size through genetic (ima3 RNAi) or physical (laser mediated extrusion) perturbations. Thus, a reduction in cell size coordinates a developmental transition in *C. elegans* embryos from asymmetric to symmetric cell division.

**Figure 3.** Size-regulated symmetry breaking in single cells. (**a**) A phase-diagram for the PAR system, considering polarisation state as a function of the circumferential length of the embryo ("Length"), and the ratio of aPAR to pPAR pool size ("AP ratio"). The diagram demonstrates a bipolar state (grey region) becomes unstable below a critical circumferential embryo length (pink and blue regions). Schematics for each of the three states are overlayed, with aPARs denoted in pink and pPARs denoted in blue. This bifurcation point quantitatively matches the critical size for which dividing P-cells in the early *C. elegans* embryo transition from asymmetric to symmetric division. Figure adapted from the work in [56]. Adjacent is the feedback motif that drives pattern formation. (**b**) In the budding yeast (*S. cerevisiae*), cell cycle commitment to Start is linked to the localisation of Cdc42 effectors at the presumptive bud site [58]. Cdc42 polarity establishment is related to the duration of the *G*<sup>1</sup> phase of the cell cycle, which ends at a critical cell size [59]. Models have shown that a growth process with positive feedback leads to Cdc42 polarisation at a single site [31].

#### *4.2. Size-Dependent Polarity Establishment in Budding Yeast*

Cell cycle commitment to budding in *S. cerevisiae* follows from Cdc42 polarity establishment at the presumptive budding site (Figure 3b). The small Rho GTPase Cdc42-GTP forms a polarity patch to mark the bud location [60,61]. This polarity pattern emerges from an autocatalytic positive feedback via Bem1 on the clustering of slowly diffusing membrane bound Cdc42-GTP [12], while the cytosolic Cdc42-GDP diffuses fast (Figure 3b). Polarity establishment in the system can be captured by mass-conserved substrate-depletion model, with a slow-diffusing activator and a fast-diffusing

substrate [26,31,45]. With appropriate choice of parameters, activator-substrate models would predict the formation of a polarised cluster beyond a critical cell size [44] (Figure 1). Thus, the establishment of Cdc42 polarisation can be linked to a critical cell size, consistent with models of critical cell size threshold at the termination of G1 phase of the cell cycle [59,62]. Molecular rewiring experiments have shown that when the Bem1 is tweaked to diffuse very slowly, multiple Cdc42 polarity patches are formed [63]. As the onset of pattern formation depends on *L*/*LD*, with *LD* the diffusion length, slowing down diffusion is equivalent to increasing the system size, so multiple patterns emerge in accordance with predictions from increasing domain length in our minimal RD model (Figure 2a,b).

#### **5. Sequential Pattern Formation and Polarisation Can Be Coordinated by a Growing Domain**

Across scales of biological organisation, structures of a bipolar or iterative nature are abundant, and their development is often sequential. While clock-and-wavefront models have successfully explained sequential patterning of axial segmentation in vertebrates [64] and recently also some invertebrates [65], sequential pattern formation can also arise from collective decision-making. Our minimal mass-conserved substrate-depletion model illustrates how an RD mechanism can give rise to sequential periodic patterning: smaller domain sizes can sustain only a single pattern via an instability of the homogeneous state, whereas domain growth provides sufficient space to accommodate multiple patterns (Figure 2a,b). Here, we implicitly assume that the RD system relaxes much faster than the timescale of domain growth, and that subunit concentration remains unchanged. When domain growth rate is comparable to the reaction rate, different dynamic patterns emerge as discussed in Section 6. In spite of differences in their mechanistic bases, sequential patterning through domain growth is common among many RD models (activator–inhibitor and substrate-depletion) [23]. In the following, we expand our focus beyond just mass-conserved models, to illustrate how local-activation and long-range inhibition can explain how growth can couple developmental tempo to state transitions.

#### *5.1. Neuronal Sequential Bipolarisation Coordinated by Membrane Growth*

Neuronal polarisation is critical for brain development. Polarisation commences as soon as neurones complete their final division, by a process of neurite formation and selection. In vitro studies have suggested neurones acquire a bipolar phenotype, generating a leading neurite, key in guiding migration, and a trailing neurite which later acquires axonal fate [66,67]. Bipolarisation in vitro is achieved stochastically, whereby the position of the first neurite is seemingly random, with the second being positioned opposite to the first [68]. How are these patterns coordinated?

Menchón et al. [69] proposed an activator–inhibitor Turing model for cell polarisation. They argued that the necessary feedback architecture for a Turing instability is manifest in developing neurones: integral membrane proteins (the polarisation cue) undergo cooperative self-recruitment, i.e., local activation, and also recruit more diffusive endocytosis modulators which facilitate their removal, i.e., long-range inhibition (Figure 4a). Indeed, in the right parameter regime in a finite domain, simulations suggest neurones can spontaneously break symmetry. The polarity regime is critically dependent on membrane size: a subcritical size prohibits symmetry breaking (like in *C. elegans*); an intermediate size allows for a single polarity axis and larger sizes allow for a bipolar phenotype (Figure 4a). Sequential and "mirrored" polarisation can be achieved by membrane growth. A growing domain leads to a time-dependent bifurcation, whereby the cell transitions from a unipolar to bipolar stability regime. The "mirroring" of the second neurite on the first can be rationalised in terms of the feedback circuit: the region of membrane furthest from the first neurite will display the lowest concentration of inhibitor. Neurones thus coordinate the developmental timing of bipolarity through a size-dependent process.

**Figure 4.** Sequential pattern formation in developmental systems. (**a**) In in vitro cultured neurones, polarity arises sequentially. A second polarity axis is formed after cell growth, and its orientation is "mirrored" off the first. A putative symmetry breaking circuit is presented adjacently [69], considering a membrane-protein (MP) activator coupled to modulators of endocytosis (ME), representing the effects of small GTPases. (**b**) Gastruloids polarise and elongate only when initialised with a critical number of cells [70]. For seed numbers beyond this initial bifurcation value, gastruloids can self-organise more axes. T/Brachyury expression is localised to the protrusion in monopolarised gastruloids, and is speculated to also be localised to further protrusions in multipolar variants. A potential feedback circuit is drawn adjacently, which remains to be investigated. (**c**) An activator–substrate model for lung branching [25], based on autocatalytic production of the signalling molecule Shh (activator) at the lung bud tip, via consumption of the substrate molecule Fgf10. (**d**) A dot-stripe mechanism is proposed to pattern the joints of developing digits: a Turing-like dot-forming system specifies the positions of bones and orients through repression a Turing-like stripe-forming system to specify joints. Modelled on a growing domain, sequential joint specification emerges, with joints forming near the developing tip. A coupled Turing scheme is described adjacent, considering a dot-forming substrate-depletion module (A,S) coupled to a stripe-forming activator–inhibitor module (B,I).

#### *5.2. Size-Dependent Sequential Patterning in Mammalian Development—Insights from Gastruloids*

Unlike in *C. elegans*, establishment of anterior–posterior polarity in the epiblast of mammalian embryos occurs well after the first cell division, an axis that lays the ground plan for the commitment of germ layers during gastrulation. In mice, AP symmetry breaking has long thought to be coordinated by the positioning of extra-embryonic cues to the posterior and hence specifying the future primitive streak [71]. This view of sequential polarity hand-off has been thrown into question in recent years by several in vitro systems, suggesting that epiblast has the capacity to break symmetry spontaneously in the absence of extra-embryonic cues [70,72–75]. While the precise genetic constituents of this symmetry breaking are under contention [76,77], several that argue some form of reaction–diffusion system is at play, citing for example the co-expression of morphogens with their extra-cellular antagonists, e.g., Wnt and its antagonist Dkk [78,79].

Consistent with this hypothesis, in vitro systems display size-dependence in symmetry breaking capacity and patterning modality [70,75]. This is seen in gastruloids, small aggregates of embryonic stem cells (ESCs) that can spontaneously break symmetry [70], axially elongating and displaying polarised expression of primitive streak marker T/Brachyury. In refining their protocol for generating gastruloids, van den Brink et al. [70] found that seeding microwells with different numbers of ESCs yielded qualitatively different phenomenology (Figure 4b): critically small (≤200 cells) aggregates could not break symmetry, aggregates of intermediate size (∼300–400 cells) could subsume a unipolar state, aggregates of double that size (∼800 cells) displayed two oppositely positioned poles, and critically large aggregates (>1600 cells) generated many poles. These results are consistent with a Turing-like system controlling polarity, whereby critically small domains cannot sustain an instability, whereas increasingly large domains can maintain progressively more "peaks". This is furthered by the observation that bipolar gastruloids polarise sequentially, with the second pole seemingly emerging after growth, protruding from the opposite edge of the structure.

#### *5.3. Sequential Patterning of Phalanges in Developing Digits Is Coordinated by Coupling Patterning to Growth*

An analogous mechanism may explain the sequential specification of joints in developing digits of tetrapod limbs. Digit patterning and growth are concomitant, with joints being laid down sequentially as progenitor cells are added to the distal tip. Guided by gene expression patterns, as well as mutant phenotypes, joint patterning has been proposed to be governed by a coupled Turing system [80] (Figure 4d). An activator–substrate system specifies the positions of bones (phalanges) by prescribing a series of "dots" of gene expression, which represses a second activator–inhibitor system to specify joints as "stripes" of gene expression at alternate positions. Simulations on a static domain recapitulate both wild type and mutant expression patterns, but patterning occurs simultaneously across the entire digit. However, simulated on a growing domain, adding new cells distally, leads to a shift in dynamics in favour of sequential patterning. Therefore, here too, the coordination of developmental timing and growth may be an emergent property of patterning by collective decision-making.

#### **6. Regulating Pattern Size and Lifetime in Growing Systems**

Pattern-forming systems that align with activator–substrate or activator–inhibitor motifs are able to undergo sequential transitions in pattern concomitant with domain growth (Figures 2 and 4). As domain growth continues, patterns undergo further bifurcations to establish periodicity. While these motifs allow irreversible transitions in pattern, it may be desirable for systems to sense intermediate sizes, and for these transitions to be reversible. One can conceive of a timer-like set-up in a biphasic scheme: in the assembly stage, symmetry is broken at a critical size, and in the proceeding dissolution stage, patterns are lost at some larger size. To investigate the emergence of timer-like behaviour, we couple an activator–substrate system to a growing domain (Figure 5a).

Specifically, we consider isotropic growth of the domain [81] (all parts of the domain grow in a similar fashion) and the subunit pool density *P* grows homogeneously with a constant rate. Due to domain growth, both the system size *L*(*t*) and the total amount of building block pool, *N*(*t*) = (*S* + *P*)*dx*, are time-dependent. The growth of the system introduces local flow and dilution of both *S* and *P* [81]. The coupled dynamics of *S* and *P* are given by (Figure 5a)

$$
\partial\_t S + \frac{\dot{r}}{r} \left( \mathbf{x} \partial\_x S + S \right) = D\_s \partial\_x^2 S + k\_{\rm on} P f(S) - k\_{\rm off} S \ , \tag{4}
$$

*Cells* **2020**, *9*, 1646

$$
\partial\_t P + \frac{\dot{r}}{r} \left( \mathbf{x} \partial\_{\mathbf{x}} P + P \right) = D\_{\mathcal{P}} \partial\_{\mathbf{x}}^2 P - k\_{\text{on}} P f(\mathcal{S}) + k\_{\text{off}} \mathcal{S} + \mathcal{G} \tag{5}
$$

where *G* is the growth rate of the subunit pool and the domain growth function *r*(*t*) is defined as *L*(*t*) = *L*(0)*r*(*t*), where *L*(0) is the initial system size. We write the assembly rate function as *f*(*S*) = *κ*<sup>0</sup> + *Sn*/(*S<sup>n</sup>* + *S<sup>n</sup>* <sup>0</sup> ), where *κ*<sup>0</sup> defines the size-independent rate of assembly of *S*. Motivated by exponentially growing cells and tissues, we specifically consider the case of exponential growth in system size, such that *r*(*t*) = *eα<sup>t</sup>* with *α* the growth rate. As the macromolecular composition of cells scales with the cell size, we assume that the total amount of building blocks, *Pdx*, grows at a rate proportional to system size *L*, resulting in a constant rate of growth *G*. While the formation of patterns occurs beyond a critical system size *L* > *L*∗, the stability of the pattern depends on the interplay between the rates of growth-induced dilution, synthesis of the subunit pool *P*, and autocatalysis of *S*.

#### *6.1. Case 1: Transient Polarity Pattern Due to Growth-Induced Dilution*

When the subunit pool grows at a rate much slower than the rate of system growth *<sup>G</sup>*˜ *<sup>α</sup>* (where *G*˜ = *GL*(0)), the structure is formed transiently and dissolves after a critical time *Tc*. The initial slow growth of the system size allows the formation of pattern beyond a critical size *L*∗. However, as *dL*/*dt* increases rapidly (due to exponential growth) and becomes much faster than the rate of pool synthesis, the subunit density starts decreasing. Below a critical density of subunits, the pattern dissolves and the system reaches a homogeneous state (Figure 5b, left). Transient structure formation has been observed in slime mould [82] and during mammalian development [83], and modelled using stochastic RD systems [84]. Here, we argue that system growth can also induce such transient polarity formation.

The patterned state makes a transition to a homogeneous state at a time *T<sup>c</sup>* when the system size reaches *Lc*. The critical time for the transition to the homogeneous state, *Tc*, is determined by the parameters of the feedback motif (*k*on/*k*off, *κ*0) and the growth rates *α* and *G*˜ (Figure 5c). The lifetime of the pattern, *Tc*, and the system size at transition to the homogeneous state, *Lc*, can be tuned independently of each other by modulating *κ*0, *α* and *G*˜ (Figure 5d). Controlling the lifetime of polarity patterns is essential for regulating developmental transitions, and further experimental studies are essential to uncover such control mechanisms.

#### *6.2. Case 2: Pattern Scaling Due to Proportional Growth of System Size and the Subunit Pool*

When the subunit pool grows at a rate comparable to the rate of growth of system size, *<sup>G</sup>*˜ <sup>∼</sup> *<sup>α</sup>*, the subunits can reach a homeostatic density in time. As a result, the patterns formed during growth do not dissolve. Strong autocatalytic growth prevents delocalisation of the early pattern and prevents the possibility of period doubling as seen in Schnakenberg kinetics and Gierer and Meinhardt model [81]. This leads to a dynamic pattern scaling behaviour where the size of the pattern scales with the size of the system (Figure 5b, middle). This mechanism of scaling is notably different from the morphogen gradient scaling [85]. Here, pattern scaling is a consequence of system growth where the polarity pattern that does not change qualitatively with system size.

#### *6.3. Case 3: Pattern Splitting*

When the rate of autocatalysis is sufficiently high, *κ*<sup>0</sup> = 0, and the subunit pool and system size grow at similar rates, *G*˜ *α*, dynamic pattern splitting can emerge in the context of growth (Figure 5b, right). Here, slower positive feedback weakens the long-range inhibition arising as the consequence of pool depletion, allowing new peaks to emerge, in contrast to case 2. Adaptive benefits of pattern splitting may be multiple, for example, allowing for the emergence of sequential patterning, and in maintaining relative stasis in local patterns upon domain extension.

**Figure 5.** Pattern scaling, splitting and transient pattern formation in growing domains. (**a**) Feedback motif for an activator–substrate system coupled to a growing domain. (**b**) (Left) When subunits are produced at a rate slower than the rate of domain growth, growth-induced dilution leads to transient pattern establishment. (Middle) When the production of subunit pool occurs at a rate comparable to system size growth, the pattern formed grows in proportion to system size, exhibiting a dynamic scaling behaviour. This is different from sequential pattern formation as the polarity is preserved during growth. (Right) In the case of strong autocatalysis of *S*, the pattern spontaneously splits. (**c**) Phase diagram for pattern formation as functions of pool growth rate relative to the system, *G*˜/*α*, and *k*on/*k*off. Colormap denotes the inverse of the pattern lifetime, 1/*Tc*. (**d**) Time evolution of structure size *<sup>S</sup>*tot = *<sup>L</sup>* <sup>0</sup> *Sdx* (blue) and system size *L* (red) for the case of transient pattern formation. The lifetime of the pattern *T<sup>c</sup>* is coupled to the system size at transition to the homogeneous state, *Lc*. They can be tuned independently of each other, for example, by changing growth rate *α* (case B) where only the transition time changes, or by changing autocatalysis rate (case C) where *T<sup>c</sup>* remains the same but transition happens at a different system size. (**e**) Tunability of pattern lifetime can be utilised as a control mechanism for symmetric and asymmetric cell division (in terms of polarity protein content). When the pattern is transient (left) the dissolution of the structure will make the daughter cells symmetric in fate, containing the same amount of polarity proteins. If the pattern persists (right) then the division will lead to asymmetric fate inheritance. Parameter values: *DP* = 1, *DS* = 0.005, *κ*<sup>0</sup> = 0.85, *S*2 <sup>0</sup> = 10, *α* = 0.01, *L*(0) = 1, and total pool density (*P* + *S*)*dx*/*L* = 2, with *G* and *k*on/*k*off variable.

#### **7. Using Growth as a Timer: Transient Symmetry Breaking at Intermediate Size**

Analysis of our minimal model shows that a timer-like control of pattern is an emergent property of symmetry-breaking systems that couple pool size to system volume. We speculate that transient symmetry breaking may serve as a control strategy to mediate shifts between symmetric and asymmetric cell division in stem cell homeostasis. In particular, the biphasic nature of transient symmetry breaking scheme can help rationalise the equal distribution of determinants in the face of unequal nature of cell division. Suppose this system underlies the establishment of cell polarity required for asymmetric division. A critical cell size would allow polarity to be established, preventing precocious cell division. After the cell size timer has elapsed, and the cell enters the dissolution phase, polarity proteins will return to the fast-mixing pool (Figure 5e, left). If cell polarity can have some effect on differential daughter cell fate, such a scheme would dissolve any bias prior to division.

The significance of a mechanism like this can be understood in cases where cell lineages undergo switches between symmetric and asymmetric cell division. Cases 1 and 2 of our model are identical besides the relative rates of pool synthesis, system size growth and autocatalysis. Therefore, tuning the coupling between pool production rate and system size can regulate transitions between reversible polarisation (case 1) and irreversible polarisation (case 2), wherein polarity is maintained through a cell division event, maintaining a bias in determinants and seeding differential daughter cell fate (Figure 5e, right). In situations where such switches between symmetric and asymmetric division are dynamic, regulating the extent of growth-induced dilution to effect switches between transient versus irreversible polarisation may be an optimal control strategy.

Given the switch between reversible and irreversible symmetry breaking in our minimal model is governed by a single parameter change, it is tempting to speculate that such a mechanism may be responsible for the mixed modes of stem cell proliferation, which in many systems are seemingly stochastic. Under such a scenario, the switch between these division modes may be noisy at the level of individual cell decisions, given the requirement of being poised near the transition point. However, given fate decisions and patterns of cell division are known to be influenced by signals emanating from stem cells or their progeny in many systems, such a model would confer plasticity and robustness in stem cell homeostasis at the level of the population. We stress that this mechanism remains a theoretical prediction and are intrigued as to whether such a control strategy is indeed utilised in nature.

#### **8. Overcoming Size Constraints: Scaling Patterns in Growing Systems**

Not all biological systems that display symmetry breaking also show size-dependent pattern formation. Indeed, it may be adaptive for systems to canalise their patterning mode irrespective of size. This is a feature of case 2 of our minimal model (Figure 5b, middle), which features commensurate growth of pool size and system size, leading to scale invariance upon domain growth: if the system breaks symmetry to form a single structure at a smaller size, upon isotropic growth, the system maintains a single structure which grows in proportion to the domain as a whole. This motif utilises the symmetry breaking capacity of reaction–diffusion systems but subverts the feature of intrinsic wavelengths characteristic of traditional Turing circuits. In this section, we delineate two potential modes of scale-invariant symmetry breaking systems—one with history dependence, and one without—and argue that such systems display adaptive features in certain contexts.

#### *8.1. Autocatalysis as a Mechanism to Preserve Patterns in the Face of Growth*

Patterning in developing systems is almost invariably proceeded by growth, which is often proportional to the initial pattern. Traditional morphogen gradient hypotheses [86] have implicitly assumed that the tissue is initially patterned when it is small and subsequently undergoes growth, facilitating proportional extension of the pattern. This two-phase model of patterning, where cell fates are assigned during an initial patterning phase, face the challenge of noise: while growth can lock in the lower positional error entailed by patterning in small fields of cells, this error cannot be reduced

through growth. Accordingly, small errors in boundary position can be amplified upon growth, demanding the read-out of positional information at early stages is exquisitely tuned. If, however, fates are assigned in a self-organised manner, as in symmetry-breaking the systems we overview in our minimal model, this hard limit on noise in boundary positioning can be surpassed. Provided the symmetry breaking system scales with domain size, absolute noise in boundary position if anything reduces with growth; self-organised systems such as these continually refine boundaries throughout growth, rather than amplifying noise in initial specification.

Case 2 of our minimal model allows for pattern scaling via proportional growth of pool and domain size, and autocatalysis, which in effect instils history-dependence in pattern formation, thus helping to preserve proportions. Therefore, the pattern generated at small domain sizes is preserved upon elongation. Given the diffusion length scale shortens with respect to relative domain size upon growth, boundaries sharpen over time. In the context of a developing field of cells, this autocatalysis could represent positive feedback in master transcription factors or indeed epigenetic changes, which allow cells and their progeny to remember past states. Therefore, such a model may help provide alternative mechanisms for scaling of patterns with growth, whereby initial stages establish the crude pattern (e.g., number and position of structures), which is in turn refined over time. The hallmark of actively scaling processes such as these is the reduction of noise in boundary positioning, i.e., violating the data-processing inequality [87].

#### *8.2. Expander-Coupled Systems Can Scale Patterns to Domain Size Irrespective of History*

While symmetry-breaking schemes incorporating autocatalysis show benefits of maintaining patterns with growth, certain systems may require patterning to be scale-invariant without the requirement of time-dependence. This is exemplified in regenerating systems, which are able to regrow organs or entire organisms in the correct proportions, in spite of drastically different starting sizes. Recent theoretical work has advanced understandings of how scale-invariant symmetry breaking could operate. Werner et al. [88] proposed that a third component is required, analogous to expanders in morphogen gradient scaling, which dynamically modulates patterning wavelength as a function of system size by tuning levels of pattern forming molecules. This model demonstrated time-independent scaling across several orders of magnitude differences in domain size. While the model is based on a traditional activator-inhibitor model, the scheme is generalisable to other modes of expander-mediated modulation and other symmetry breaking motifs such as substrate depletion.

#### **9. Discussion**

In this perspective, we presented a minimal model for symmetry breaking to serve as a unifying framework to understand pattern formation in the context of timing and growth. We argue that systems that display positive feedback in activator recruitment, drawn from a limiting pool, can yield spontaneous symmetry breaking. This basic scheme is mathematically akin to other RD mechanisms including activator–inhibitor or substrate-depletion motifs, all relying on a common logic of local activation and long-range inhibition. Thus, the insights gleaned from the phenomenological behaviour of this system is applicable to diverse systems.

Across the cases of the minimal model we consider, we observe a hard size limit on pattern formation: below a critical size, diffusive dispersion overwhelms the capacity to break symmetry. Given developing systems across scales typically display patterning and growth occurring in unison, if this critical size is within biologically meaningful length scales, such behaviour can elicit qualitative changes in patterning: growth above a critical size leads to a bifurcation, whereby the system transitions from unpolarised to polarised. Viewing growth as a control parameter of the system that increases system size at a predictable rate, developmental systems can utilise this bifurcation to enact developmental transitions at the right place and time. As a generic by-product of symmetry breaking systems, we predict that this time-keeping mechanism may be more abundant than anticipated. We note that

this feature is the most generic among RD models of pattern formation, and among the different cases of our minimal model: the diffusion length-scale sets a physical limit on pattern formation.

Beyond this first bifurcation, our mass-conserved RD model predicts different dynamic behaviours depending on the regulatory motifs. These include sequential pattern formation, transient pattern formation, pattern scaling, and pattern splitting in growing systems. In line with the well-established literature on domain size in Turing patterns [23], our model predicts clock-like sequential pattern formation (Figures 2 and 4): as the system grows larger, given patterning wavelength is intrinsic to the system, the domain can accommodate multiple structures. An important dynamical consequence of this is temporal ordering: growth elicits consecutive bifurcations, resulting in sequential patterning, shown to be instrumental in neuronal cell (bi)polarity [69], and joint patterning in digits [80]. Alternatively, growth-induced pool dilution can drive systems back towards an unpolarised state, allowing for transient pattern formation at intermediate size (Figure 5). Thus, an alteration in growth regulation can yield timer-like dynamics, which we hypothesise may be important in orchestrating switches between asymmetric and symmetric stem cell division modes. A qualitatively different behaviour upon continued growth is scale invariance, whereby the proportions of the pattern are maintained upon domain elongation. Scale-invariant systems show switch-like dynamics, becoming time-independent after the first bifurcation. We argue that such behaviour could allow patterned tissues to maintain proportions upon proliferation, where self-organisation continually refines boundary position instead of stretching noise in initial specification.

Our reaction–diffusion framework for understanding developmental time in terms of size-dependent symmetry breaking is generalisable beyond the systems that couple increases in size to developmental transitions via biochemical circuits. First, decreases in system size can also be utilised by developmental systems to temporal transitions. For example, the transition from asymmetric to symmetric division in the P-lineage of *C. elegans* can be understood in terms of sequential reductions in cell volume pushing the system over the critical cell size threshold for polarisation. Second, the organising principle of Turing-like pattern formation—local-activation and long-range inhibition—extends beyond systems based solely on chemical cross-talk [38]: pattern formation can emerge from cell–cell interactions or mechanical instabilities [89–91]. While we restricted our focus in this paper to biochemical systems, future work should attempt to unify these results with mechanically driven size-dependent symmetry breaking. Indeed, we may see strong parallels in how nature utilises chemical or mechanical instabilities to regulate the timing of developmental transitions. We hope that our proposed strategies for time-keeping in natural living systems can also provide inspirations for engineering of synthetic circuits with tunable dynamics.

**Author Contributions:** Conceptualization, J.C.S., D.S.B. and S.B.; methodology, D.S.B.; validation, J.C.S., D.S.B. and S.B.; formal analysis, D.S.B.; investigation, J.C.S., D.S.B. and S.B.; resources, S.B.; data curation, J.C.S., and D.S.B.; writing—original draft preparation, J.C.S., D.S.B. and S.B.; writing—review and editing, J.C.S., D.S.B. and S.B.; visualization, J.C.S., D.S.B. and S.B.; supervision, S.B.; project administration, S.B.; funding acquisition, S.B. All authors have read and agreed to the published version of the manuscript.

**Funding:** SB acknowledges funding from Royal Society University Research Fellowship URF/R1/180187, Human Frontiers Science Program (HFSP) grant number RGY0073/2018, and Carnegie Mellon startup funds.

**Acknowledgments:** The authors thank Andrew Goryachev and the anonymous reviewers for many useful comments and suggestions.

**Conflicts of Interest:** The authors declare no conflicts of interest. The funders had no role in the design of the study; in the collection, analyses or interpretation of data; in the writing of the manuscript; or in the decision to publish the results.

#### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

#### *Article*
