Robustness and timing of cellular differentiation through population-based symmetry breaking

Functional diversification of cell types during mammalian development is characterised by the transition from an initially homogeneous group of multilineage primed cells towards a heterogeneous population of differentiated cell types (Zhang and Hiiragi, 2018; Simon et al., 2018). To ensure robust development, the onset of the differentiation event must be accurately timed, and the number distribution of each cell type must be correctly established.

The experimental observation that the expression of mutually exclusive genetic markers distinguishes the differentiated cell types from each other, and from the multilineage primed state, has promoted the hypothesis that cell types correspond to one of multiple stable gene expression states that arise from intracellular gene regulatory networks. Switching between these distinct states that dynamically represent stable attractors, specifically from the attractor encoding the multilineage primed state to the differentiated states, has been assumed as the dynamical basis of differentiation (Kauffman, 1969; Andrecut et al., 2011; Wang et al., 2011; Enver et al., 2009). The most common motif that accounts for bistability, i.e. the co-existence of two stable gene expression patterns on a single-cell level, is a two-component toggle-switch gene network (Thomas, 1981; Cherry and Adler, 2000; Snoussi, 1998). The addition of self-activating loops to each of the toggle-switch genes gives rise to a third stable state through which the multilineage primed co-expression state has been typically explained (Huang et al., 2007; Bessonnard et al., 2014; Jia et al., 2017). Such single-cell multistable circuits have been used to describe, for example, the Cdx2-Oct4 (also known as Pou5f1) switch in the differentiation of totipotent cells of the early embryo (Niwa et al., 2005), as well as the Gata6-Nanog switch in the differentiation of cells in the inner cell mass (Bessonnard et al., 2014; Chickarmane and Peterson, 2008; Schröter et al., 2015). The differentiation outcomes have thereby been typically analysed from a single-cell perspective, assuming that the intrinsic transitions towards one of the co-existing stable states are either driven stochastically (Gupta et al., 2011), modulated by an external signalling input, or stemming from cell-to-cell heterogeneity (De Mot et al., 2016; De Caluwé et al., 2019). However, the respective cell-type proportions strongly depend not only on the parameters of the system but also the initial conditions, indicating that the underlying mechanism of reliable proportioning cannot be directly inferred from this view.

On the other hand, experimental evidence suggests that paracrine signals have a crucial role for reliable cell-type establishment (Monk, 1997; Nichols et al., 2009; Yamanaka et al., 2010; Youk and Lim, 2014; Hart et al., 2014; Saiz et al., 2020). This has been particularly recognised in lateral inhibition models, such as the Delta-Notch system (Sprinzak et al., 2010; Collier et al., 1996), in which the intercellular communication realised through toggle switches distributed between the cells enables differentiation into distinct cell types with specific proportions and spatial organisation (Collier et al., 1996; Formosa-Jordan and Ibañes, 2014; Teomy et al., 2019 preprint; O'Dea and King, 2012). A similar model relying on cell-cell communication via fibroblast growth factors has been recently proposed for the cell-type specification in the mouse blastocyst (Saiz et al., 2020).

The emergence of specific gene expression profiles due to cell-cell communication is also well exemplified by the quorum-sensing mechanism in bacteria, in which the timing of emergence is regulated through the concentration of secreted signalling molecules (Taga and Bassler, 2003; You et al., 2004; Danino et al., 2010). Thus, these examples indicate that novel attractors can emerge in populations of communicating cells that are not present for isolated cells. This principle of attractor emergence in coupled systems has been extensively investigated in both natural and synthetic genetic networks (McMillen et al., 2002; Taga and Bassler, 2003; Kuznetsov et al., 2004; Garcia-Ojalvo et al., 2004; Ullner et al., 2007). For example, emergent oscillatory solutions (Ullner et al., 2007, 2008; Koseska et al., 2010) have been proposed as dynamical mechanisms of stem cell differentiation with self-renewal (Suzuki et al., 2011). These theoretical and experimental findings suggest that cooperative dynamics could underlie the emergence of differentiated cell types in cell populations. However, the principles that govern self-organised differentiation timing, while at the same time leading to specific cell-type proportions that can reliably adapt upon perturbations, remain unresolved.

We propose a dynamical mechanism that underlies the differentiation of cell populations into cell types with stable proportions driven by the growth of the communicating population itself. Within this model, the expression profiles of both the multilineage primed and the differentiated cell types can be captured without any change in the parameters, and the transition between them can be established at a specific size of the population. The same mechanism also enables robust recovery of the acquired proportions upon perturbation of the system, accounting for autonomous scaling. We show that these properties of a communication-based cellular system arise from a collective symmetry breaking through a subcritical pitchfork bifurcation (PB). Dynamically, this novel collective state of mutually exclusive gene expression states in the population is represented as a single inhomogeneous steady state (IHSS). This renders the individual differentiated cell types dependent on each other, and the respective gene expression profiles different from those that can be attained in isolated cells. We show that the differentiation timing occurs through the growth of the population beyond a critical size, when the parameters of the undifferentiated cell population correspond to organisation before the PB. These findings indicate that the proposed collective inhomogeneous solution is a generic dynamical mechanism that describes how heterogeneous cellular types emerge from a homogeneous population of precursor cells and are maintained via cell-cell communication.

Cell differentiation during early mammalian embryogenesis often occurs from a multilineage primed state (mlp), which is maintained for several cell division cycles before differentiation into distinct cell types occurs (Fig. 1A) (Saiz et al., 2016; Hatakeyama et al., 2004; Soldatov et al., 2019). This implies that gene expression states transit from an initial homogeneous pattern across the precursor cells towards a heterogeneous expression state representing the differentiated cell types.

Employing this model, we generated a numerical lineage tree starting from one cell, using stochastic simulations. The population growth was implemented in a simplified manner: after a given time period that mimics cell cycle length, all cells divide and the number of cells is doubled. The initial gene expression states of the daughter cells are inherited from the final state of the mother cell. The cells are placed on a grid with no flux boundary conditions, and the communication between them is short-ranged (within distance 2), i.e. between adjacent and second-adjacent cells on the grid (depicted for N=8 cells in Fig. 1C, left). Although we mainly focus on this communication type, the effect of other coupling types is also explored further below: global all-to-all, local (only between adjacent cell) and distance-based probabilistic coupling, in which interaction links between cells are established with decreasing probability as the distance between them increases (Fig. 1C). Cell divisions occur along the horizontal and vertical axes of the grid alternately, such that the grid doubles in size after each division event, yielding lattices of 1×1, 1×2, 2×2, 2×4, etc. (Fig. S1). The collective system state constituted of u+, v + or mlp cell types was characterised in every time instance by categorising the cell states, and was used to plot the temporal evolution of the lineage tree (Fig. 1D). Furthermore, cellular proportions in the system were estimated from the collective state at each time point, and their temporal evolution is shown in the panel above the lineage tree.

The simulations demonstrated that even in the presence of gene expression noise (σ=0.1, Materials and Methods), for a population size of up to N=4 cells, the dynamical state of the system was homogeneous such that individual cells had equivalent u/v co-expression patterns resembling an mlp state (Fig. 1D). However, when the population reached a threshold size of N=8 cells, the expression patterns in single cells collectively transitioned to u- or v-dominated expression (u+ or v+ cells), indicating that the initial symmetry of the system had been broken and cells differentiated. Interestingly, defined proportions of u+ and v + cells emerged upon differentiation, which were stably maintained over many rounds of cell division (Fig. 1D, top).

That the mlp state could be maintained for smaller population sizes, despite the presence of gene expression noise, indicates that stochastic events do not trigger the transition from a homogeneous to a heterogeneous cell population. As the model parameters also do not change, these results rather suggest that the timing of the transition event emerges from the growth of the population, upon which distinct proportions of the generated heterogeneous cell types are established.

To uncover the dynamical mechanism that underlies the transition between the mlp and the u+/v+ cell types, we performed a bifurcation analysis. The results in Fig. 1D indicated that differentiated cell types emerge with growing population size. However, a direct identification of the transition type using the number of cells as a bifurcation parameter cannot be performed, as it is not an explicit parameter in the model (Eqn 1). As a first step, we therefore chose α_u as a representative bifurcation parameter to identify and characterise the solutions of the system, whereas the impact of population size will be explored further below.

The bifurcation analysis showed that for N=2 coupled cells, the system exhibits qualitatively different dynamics for different α_u values. For lower α_u values, the system exhibits monostable behaviour (Fig. 2A, top, green circle for α_u=2.3). Here, both cells populate the same state (u₁=u₂) as evident from the projection that falls on the diagonal in the u₁−u₂ state space (Fig. 2A, bottom left). This homogeneous steady state (HSS) captures the mlp state with a characteristic precursor u/v co-expression pattern. On the level of a single-cell system, this is the only stable regime over the full α_u domain (Fig. S2A), meaning that isolated cells will maintain the mlp state.

However, at a critical α_u value, the symmetry of the mlp HSS is broken via a subcritical PB. From the point of view of dynamical systems, the HSS (Fig. 2A, black line) loses its stability at this parameter value and a pair of fixed points is generated, giving rise to IHSSs stabilised via saddle-node bifurcations (SN) (Fig. 2A, purple lines). There is a co-existence between the HSS and the IHSS before the PB, rendering it subcritical. The IHSS is a single attractor, here a six-dimensional point (u₁, v₁, s₁, u₂, v₂, s₂), that describes a heterogeneous state. This state consists of a high u-expression level (u+) in one cell (u₂) and a low u-expression level (v+) in the other cell (u₁), as shown by the two-dimensional (u₁, u₂) projection (Fig. 2A, bottom right). This means that the IHSS is a symmetry-broken collective state (Koseska et al., 2013), as it describes simultaneously both cell types with mutually exclusive gene expression patterns (Fig. S2B,C). As there is no preference which cell will acquire the u+ or v + cell type, both possibilities (u₁<u₂ and u₁>u₂) are present as branches (upper and lower branch in Fig. 2A, top, respectively). Therefore, the same fixed point will be manifested as an upper branch for the u+ cell, but as a lower branch within the equivalent bifurcation diagram for the v + cell. The u₁−u₂ state space projection also demonstrates that the IHSS solutions are reflections of one another over the diagonal (Fig. 2A, bottom right). This demonstrates that PB provides a unique mechanism for a dynamical transition from a homogeneous (mlp, HSS) to a single but heterogeneous (u+/v+, IHSS) state of the population, in which the differentiated cell types always jointly emerge. The described IHSS solution of a coupled system is fundamentally distinct from a bistable system on a single-cell level. There, the u+ and v + cell types are described by two different steady states. Thus, in the absence of cell-cell communication, each cell in a population can independently transition to one of them. In principle, all cells within a population could acquire u+ cell type, without the v + cell type ever occurring.

We next investigated which type of network topology could give rise to stable heterogeneous population states. We found that the IHSS emerges for a number of different network topologies that are characterised by inhibitory coupling, i.e. topologies with an effective negative feedback on u via the signalling molecule s (Fig. S3A to C; Eqn 2). Moreover, IHSS could be also obtained for oscillatory gene expression dynamics in isolated cells (Fig. S3D, Eqn 3). The bifurcation analysis additionally implies that heterogeneous cell types will be reached even when starting from identical initial conditions, representing a true symmetry breaking event.

Interestingly, the value for α_u=2.3 used in Fig. 1D lies before the PB, what we will term herein as a subcritical organisation of the system. The bifurcation diagram for N=2 cells (Fig. 2A) can therefore explain why in the lineage tree in Fig. 1D the cells maintained the mlp state and did not switch to the u+/v+ cell type, but it does not explain how the transition can occur for increasing population size N. We therefore generated a bifurcation-like diagram by varying N as in the lineage tree, while keeping all parameters fixed (with α_u=2.3) and using short-range communication (Fig. 2B). We identified the distinct collective states for each N via their u+/v+ cell-type proportions. The existence of steady states was estimated with an exhaustive stochastic search, with different initial conditions and noise intensities (Materials and Methods). The identified IHSSs with equivalent proportions were grouped, and the average u-values were depicted separately for the u+ and v + cell types within each distribution group. This yields paired upper and lower branches, analogously to the bifurcation diagram in Fig. 2A. Whereas for N=2 coupled cells only the mlp HSS was detected (Fig. 2B green; corresponding to the green circle in Fig. 2A), for N=4 coupled cells stable inhomogeneous solutions with a distinct proportion could be additionally identified (red/blue u+/v+ horizontal stack bar markers). This resembled the subcritical PB observed in Fig. 2A. The identified inhomogeneous solution corresponds to one cell having high- (u+) and three cells having low-expressing u state (v+), or 1u+/3v + distribution. However, for larger population sizes, the homogeneous mlp state lost its stability and only IHSSs reflecting the two cell types with specific proportioning were identified. In other words, the growth of the population to N=8 triggered a transition from the precursor to the differentiated cell state. This demonstrates that when α_u assumes a sufficiently low value to set the system before the bifurcation point, the growth of the population can trigger a dynamical transition, resembling the one that occurs when α_u is increased (Fig. 2A).

To explore how the population size N_SB at which there is an occurrence of symmetry breaking depends on α_u, we performed an equivalent of a two-parameter bifurcation analysis. The diagram shows that the symmetry-breaking transition could be triggered over a distinct α_u parameter region for short-range coupling (Fig. 2C, solid line). Depending on the α_u value before the PB, the symmetry-breaking point will be realised at a distinct size of the population (N_SB). Thus, for subcritical organisation, differentiation timings can emerge in a self-organised manner. Similar results were also obtained for local coupling and the probabilistic distance-based coupling (Fig. 2C). On the other hand, for global all-to-all coupling, the symmetry-breaking transition could only be triggered stochastically with an increase in N. Therefore, these results suggest that under local and short-range coupling, as the population grows in size, the PB shifts its position towards lower α_u values. This is likely caused by the relative change in the effective communication range of the cells from global to a more local one with the size increase. The PB shift thereby enables transition of the system state from an HSS to an IHSS. Altogether, the analysis renders the number of cells as an effective bifurcation parameter that in conjunction with subcritical organisation drives the timing of cellular differentiation.

To investigate how the u+/v+ cell-type proportions emerge and are stabilised as the size of the population increases, we analysed the IHSS manifestation for N>2 cells in terms of bifurcation structure. For N=4 cells, the short-range and global coupling are equivalent, due to the small system size. The bifurcation analysis showed that at α_u=2.3, although for N=2 only the mlp HSS was stable (Fig. 2A), for N=4 coupled cells there was co-existence of the HSS and IHSS (Fig. 3A, equivalent to the result in Fig. 2B). The observed IHSS distribution was 1u+/3v +, with one cell having high- (u+) and three cells having low-expressing u state (v+), which for increasing α_u values was followed by 2u+/2v +, with two cells in each state, and 3u+/1v +, with three cells with high- and one with low-expressing u state (blue, purple and red branches in Fig. 3A, respectively). All of these distributions are associated with stable attractors that emerge from the same PB. In general, for N globally coupled cells, N−1 different distributions of cells between the u+ and v + cell types are stable ([k]u+/[N−k]v +, for ⁠) (Koseska et al., 2010). These distributions are always sequentially ordered towards an increasing proportion of u+ cells for increasing α_u values. IHSS branches with adjacent distributions ([k]u+/[N−k]v + and [k+1]u+/[N−k−1]v +) overlap, whereas the more dissimilar u+/v+ distributions are separated along the bifurcation parameter domain (Fig. 3A). Thus, for a given α_u, either a single IHSS distribution with distinct cell-type proportions, or adjacent branches with similar ratios, will be predominantly visited (Fig. 3A, top).

Intrinsic cell-to-cell variability in terms of circuit parameters does not affect the ordering of IHSS distributions and thereby cell-type proportioning. This was demonstrated by a bifurcation analysis on a coupled system of four cells with non-identical α_v values (Fig. S4A). Although in this scenario there are already slightly different mlp values in each cell, the parameter range over which the IHSS branches are stable is expanded (Koseska et al., 2009), and the overlapping intervals between adjacent distribution branches are reduced (compare Fig. S4A with Fig. 3A). Thus, the role of cell-to-cell variability is fundamentally different in coupled multicellular systems compared with multistable single-cell systems: it does not cause the symmetry-breaking event, but rather enhances its manifestation. Overall, these results indicate that the reliable cell-type proportions that emerge as the system transits to the differentiated state at a critical population size, are an inherent property of the distribution of IHSS solutions.

To investigate whether for larger population sizes the u+/v+ cell type proportions are indeed related to the sequential ordering of IHSS distributions, we compared the proportions for increasing α_u values at N=4, 8 and 32 short-range coupled cells. Multiple realisations of the numerical simulations were performed when starting from initial conditions randomly drawn from a normal distribution with a mean equal to the α_u-specific mlp state and s.d. σ_ics=0.2μ_ics (Fig. 3B). For N=4, the results were identical as those obtained from the bifurcation analysis: at α_u=2.3, 1u+/3v + proportion was detected, transiting to 2u+/2v + and 3u+/1v + as α_u was increased (Fig. 3B, top), corresponding to the sequential ordering of the branches and the probabilities for visiting them (Fig. 3A).

For N=8 short-range coupled cells, only the IHSS solutions were stable at α_u=2.3 (Fig. 3B, middle), showing again that the differentiation occurs at the critical transition from N=4 to N=8 cells (red arrow in Fig. 3B; as in Fig. 2B). The PB was shifted to a lower α_u value (α_u=2.281), thus enabling the differentiation timing to be realised in a self-organised manner. Furthermore, the u+/v+ cell type proportions were stabilised as the population increased (Fig. 3B, bottom, for N=32 cells). This shows that for a given parameter organisation, defined proportions can be reliably established and reproduced through multiple simulation realisations. We also corroborated the reliability of the timing mechanism by estimating the temporal evolution of the fraction of u+ cells across the different stages of the lineage tree (Fig. 3C). This showed that differentiation timing at N=8 cells, following the third cell cycle, was achieved in 94% of the cases, but also, the u+/v+ cell type proportions were reliably and stably maintained thereafter. The reliable cell-type proportioning was also manifested for the three additional coupling scenarios (Fig. 1C). In these cases, the increase in the proportion of u+ cells with an increase in α_u for N=32 cells is analogous to the progression of branches in the generic bifurcation analysis (N=4, Fig. 3A), although the HSS was destabilised at different α_u values for each of them (Fig. S4B-D). The variability between the cell-type proportions for a fixed α_u value is nevertheless with a ≤5% s.d. For local coupling, specifically, a 50%−50% ratio was maintained in a large α_u interval, indicating that the probability of visiting a different IHSS manifestation increases only for higher α_u values (Fig. S4C).

Therefore, this analysis showed that although different initial conditions can lead to different IHSS distributions, in most cases these are distributions with similar cell-type proportions. Thus, (1) reliable cell-type proportioning is a result of the sequential ordering of IHSS solutions in parameter space and (2) cell-cell variability enhances the manifestation of the identified symmetry breaking solution.

Our demonstration that the proposed symmetry-breaking mechanism is population-based suggests that robust cell-type proportions would be generated despite differences in initial conditions, and that they can be dynamically recovered upon perturbation. To probe the robustness property, we investigated the influence of different initial conditions distributions (changing the variance or shifting the mean value) that determine the initial gene expression states in single cells, as well as the effect of gene expression noise intensity. The results were obtained for a population of fixed size, N=32 cells, under the four distinct coupling types (Fig. 1C), and a fixed α_u value before the PB as before (α_u=2.3). Sampling the single-cell initial conditions from a normal distribution with increasing s.d. around the mlp value (Fig. 4A, top) produced distributions with reliably conserved u+/v+ cell-type proportions under short-range coupling (Fig. 4A, ratios with ≤8% s.d.). This demonstrates that variance in initial gene expression at the single-cell level does not majorly affect the differentiation outcome. Rather, cell types are established in characteristic proportions within a coupled system, even when starting from broad distributions of initial conditions. However, between the different coupling types, different stable u+/v+/mlp proportions were generated for this fixed α_u, in agreement with previously estimated values (Fig. 3B; Fig. S4B-D): for short-range (Fig. 4A), 0.5 for local (Fig. S5A) and for probabilistic distance-based coupling (Fig. S5G), whereas the HSS remained stable against moderate perturbations for global coupling (Fig. S5D). Furthermore, stochastic realisations with a stepwise shift in the initial mean value from high v-expression to high u-expression state (Fig. 4B, top) also resulted in reliable cell-type proportions (Fig. 4B; Fig. S5B,E,H). The proportions were also reliable for simulations when multiplicative white noise intensity was increased (Fig. 4C; Fig. S5C,F,I, Materials and Methods). We also observed a manifestation in which, besides populating u+/v+ states within the IHSS solution, few cells also populated the mlp state, resembling a chimera-like state (Kuramoto and Battogtokh, 2002; Abrams and Strogatz, 2004). This was observed mainly for the probabilistic distance-based coupling (Fig. S5G-I), and for the short-range coupling in rare instances (Fig. 4A). As chimera states have been predominantly characterised for systems of coupled oscillators, an additional study would be required to dynamically classify these solutions.

We next explored whether this population-based symmetry-breaking mechanism could also underlie the ability of the early embryo to re-establish exact cell-type proportions following perturbations (Martinez Arias et al., 2013). To investigate this, we performed a numerical experiment in which cells were separated according to their type after the fourth cycle of the lineage tree (N=8 cells, Fig. 1D), forming two single-type subpopulations of different sizes that could further continue to grow and divide (Fig. 4D, left and right). The initial conditions of the two new lineage trees were therefore not positioned in the vicinity of the stable steady states. The subpopulation of two coupled cells with high u-expression reverted to the only stable attractor for this system size: the mlp HSS. However, after two cell cycles (reaching N=8 cells), both differentiated cell types re-emerged (Fig. 4D, left). The other subpopulation of N=6 cells with low u-expression initially transiently revisited the mlp state (attracted and then repelled by the unstable saddle-type HSS) before both cell types stably re-emerged and the population settled in an IHSS during the first cell cycle (Fig. 4D, right).

The difference in timing between the two cases again points to a cell-number dependence in the triggering of the symmetry breaking (Fig. 2B). The cell-type ratios for both subpopulations were stabilised to a steady value similar to that of the full system before separation, and differed among the two subpopulations by ∼6%. This autonomous scaling and regenerating capability of the coupled system is a direct consequence of the properties of the IHSS solution: dynamically, it is not possible to stably populate the u+ without populating the v + cell type. Thus, even when cells are separated such that only the cells of one type remain, the cell division and the cell-cell communication through which IHSS is established in the first place, will enable the system to recover both cell types with reliable ratios.

To confirm that both cell types and their proportions are indeed generated in a communication-dependent manner, we next investigated how inhibiting the cell-cell communication would affect the proportioning. We considered two different simulation scenarios, by implementing an increasing: (1) inhibition of α_s strength to mimic decreased production of the signalling molecules, and (2) inhibition of α_u,s strength, which effectively uncouples the dynamics of the intracellular circuit from the extracellular signalling. For this, α_s or α_u,s were effectively decreased by given multiplicative factors (1−s_inh,out) and (1−s_inh,in), respectively. The simulations showed that in the first scenario, decreasing the production of the signalling molecule results in an increase of the u+ cells proportions, abruptly transiting to a homogeneous population of high u-expressing cells at 25%−30% of s_inh,out (Fig. S6A). Decreasing α_u,s on the other hand, which decreases the overall transcription rate of u, reduces u+ cells proportions to abruptly transit to a homogeneous population of low u-expressing cells at of s_inh,in (Fig. S6B). The single cells are now weakly coupled and, for the given parameters, the system is organised in the monostable low u-expressing state. The lack of robustness in cell-type proportioning observed under these conditions thus underlines the importance of population-level coordination in the system.

In summary, the presented results demonstrate that a subcritical organisation in conjunction with cell division within a communicating population ensures not only the robustness of the proportions of differentiated cell types with respect to initial conditions and noise, but can also enable re-establishment of this distribution upon perturbation.

As the differentiation mechanism described here relies on cell-to-cell communication, we investigated whether the differentiated cells were randomly distributed or organised into spatial patterns. We performed stochastic lineage-tree simulations, with 13 cell cycles as in Fig. 3C, reaching a grid size of 64×64 cells. The u+/v+/mlp proportions were estimated from the stable steady states at the end of each cell cycle. To characterise the spatial distribution of cells, we quantified the extent of spatial clustering as the cluster radius of the u+ cells. Thus, for all u+ cells, the fraction of u+ cells in their surroundings with increasing distance was estimated (Fig. S7). The cluster radius was thereby the distance at which the fraction of u+ cells dropped to halfway between the fraction at zero distance (equal to 1), and the global fraction of u+ cells (dashed line in Fig. S7).

To systematically assess how the coupling range, i.e. the extent of signalling communication, contributes to the formation of patterns with distinct spatial frequencies, we used deterministic and probabilistic coupling with different ranges (1-, 2-, 3-, 5-, and 10-range, Materials and Methods). For the deterministic coupling schemes, distinct regular spatial patterns were formed, the features of which were strongly linked with the respective coupling range: local coupling generated checkerboard-like patterns, short-range coupling generated regular patterns of small u+ cell clusters (see also Fig. S8A, Movie 1), whereas stripe-like patterns were observed for increasing coupling (10-) range (see also Fig. S8B, Movie 2). Example patterns from the last stage of the lineage tree are shown in Fig. 5A. The clustering radius in these cases increased with the increase of the coupling range (Fig. 5B, left). In all the deterministic coupling scenarios, the patterns and their proportions were reliably reproduced for independent realisations of the lineage-tree simulations (Fig. 5C, left).

Interestingly, the spatial frequency of the formed stable patterns was preserved as the size of the population grew. This can be particularly exemplified for the case with 10-range coupling, in which doubling in the number of stripes followed the horizontal cell division events (Movie 2). This effectively rendered the width of the stripes, i.e. the u+ cluster radius and the distance between them, stable across the cycles (Fig. S8B). However, initially, at around the third cycle, a lower stable fraction of u+ cells was established and maintained by simple state propagation through cell division, which resulted in growth of the u+ cluster radius, until the critical system size was reached (eighth cycle) and a stable stripe-like pattern was formed (Fig. S8B). The pattern-directed u+ cluster size and cell-type proportions were further stably maintained. The simulations also showed that the cell population size at which the transition to the differentiated state was triggered was also preserved (Fig. S8A-D, lower left). Even more, the formed patterns and their characteristics were also preserved for a fixed population size as in the final grid, in which cells were randomly initialised (Fig. S9A-C).

On the other hand, when probabilistic distance-based coupling of different ranges was used, the lineage-tree simulations showed arrangement formation with less regular patterns (Movie 3 for 10-range coupling). Nevertheless, larger cluster size was again characteristic for larger communication ranges (Fig. 5B, middle; Fig. 5D). The proportioning was also reliably achieved over the different realisations (Fig. 5C, middle), in which the u+ fraction was maintained mainly by propagating the cell type from mother to daughter cells during division. Therefore, the u+ cluster radius is directed by the population growth, but is ultimately constrained by the communication range (Fig. S8C). That population growth highly regulated the formed patterns is also reflected through the formation of random spatial arrangements when populations of fixed size were randomly initialised (Fig. S9B-D).

The limiting case of these observations is the global all-to-all coupling scenario, which in most of the cases resulted in the formation of a single u+ cell-type cluster when population growth was considered (Fig. 5D, right; Fig. S8D, Movie 4). Such spatial arrangement is a direct consequence of the daughter cells inheriting their cell type from the mother cell along the lineage tree. In this case, the cell-type proportions were maintained globally, and the states of adjacent cells did not need to readjust to maintain the respective proportions locally. The relative size of the u+ cluster with respect to the population size was thereby only constrained by the cell type proportions. In contrast, for a fixed population size and starting from random initial conditions, a random arrangement of cell types across the grid was observed (Fig. S9D). Thus, in systems with global all-to-all or probabilistic coupling, the growing of the population is crucial for the observed spatial organisation.

Therefore, these results suggest that the IHSS not only represents a dynamical mechanism for generating stable proportions of differentiated cell types, but also enables their reproducible arrangement in regular patterns, the frequency of which is in turn constrained by the communication range.

We have argued here that intercellular communication, an integral process in developing mammalian embryos (Saiz et al., 2020; Lorthongpanich et al., 2012), gives rise to mutually exclusive differentiated cell types (Koseska and Bastiaens, 2017) as the population grows in size. We demonstrated that both the homogeneous undifferentiated and heterogeneous differentiated cell types, as well as the transition between them, can be described without changing the model parameters. This is valid when isolated cells are characterised with monostable steady state or oscillatory gene expression dynamics. The proposed simple but scalable mechanism enables robust cell-type proportions to emerge autonomously, but also describes their reliable recovery upon perturbation. This population-based heterogeneous solution is thereby distinct from the concept that single-cell multistability together with cell-to-cell variability are necessary to describe how differentiated cell types emerge. Such a single-cell scenario does not provide a mechanism for robust proportioning between differentiated cell types, but rather relies on tight regulation of initial states and signalling.

These two conceptual views have been recently subjected to an experimental test (Raina et al., 2020 preprint) using an in vitro model for robust proportioning of epiblast and primitive endoderm-like cell types in mouse embryonic stem cells. Starting from a broad range of initial gene expression profiles, the wild-type populations achieved robust epiblast and primitive endoderm cell-type proportions, in contrast to a communication-deficient mutant. The experiments also showed that the proportions could be reliably re-established in the wild type upon disproportionate removal of one cell type, in line with our predictions. Therefore, these observations corroborate the proposed hypothesis that balancing between robustness of differentiated cell-type proportions, while maintaining the plasticity of the system, such that it can recover to the exact proportioning upon perturbation, requires a population-based symmetry-breaking mechanism as realised by the IHSS.

Important insights regarding symmetry-breaking mechanisms unquestionably came from Turing's seminal work (Turing, 1952), and have been widely explored to describe the emergence of spatial organisation during development (Raspopovic et al., 2014; Economou et al., 2012). The IHSS similarly provides a dynamical transition from homogeneous mlp to a heterogeneous state of differentiated cell types within a population. Even more, the IHSS enables a reproducible spatial arrangement of the two differentiated cell types for a broad range of coupling scenarios, irrespective of the initial conditions. We have demonstrated that the formed pattern type is tightly related to the coupling range and the specified stable IHSS distribution. The spatial patterns identified in our simulations are broadly consistent with experimentally observed patterns in cell differentiation paradigms that could potentially be governed by an IHSS. Cell differentiation in the inner cell mass (ICM), for example, is triggered by a short-range communication signal (Raina et al., 2020 preprint). In mouse embryos and ICM organoids, the two cell types differentiating from the ICM form small clusters (Mathew et al., 2019; Fischer et al., 2020), similar to what we obtain in simulations with short-range coupling.

What is unique about the bifurcation transition presented here is not only the reliability of the differentiated cell type proportions and spatial organisation, but also the timing of the differentiation event that emerges as a result. The mechanism of differentiation timing is a unique property for organisation of the parameters of the system before the subcritical PB. A similar mechanism of population-based symmetry breaking, but via supercritical PB, has been previously suggested for the Delta-Notch lateral inhibition model, when the strength of the local interaction between the two cells is varied (Ferrell, 2012). However, in this case, the differentiation timing cannot be realised and the manifestation of only a specific cell-type proportioning with likely a spatial checkerboard-like pattern can be explained. A conceptually equivalent model to Delta-Notch has been recently described for cell differentiation in the mouse blastocyst, by relying on non-autonomous switching of gene expression circuits at a specific size of the population to recapitulate the emergence of differentiated cell types (Saiz et al., 2020). Thereby, the issue of self-organised differentiation timing characteristic for early embryo development has not been addressed.

The IHSS mechanism proposed here can be also taken one step further. An extension of the proposed mechanism can be envisioned that describes pluripotency and multipotency of stem cells. Conceptually, this would correspond to a finite cascade of subsequent PBs occurring simultaneously on both branches of the existing IHSS solutions (Zakharova et al., 2013; van Kekem and Sterk, 2019). IHSS therefore represents a cooperative dynamical mechanism through which a growing homogeneous cell population can break the symmetry, a prerequisite for novel information regarding different cellular types to emerge. Organisation in the vicinity of this dynamical transition allows the comprehensive capture of not only the differentiation timing, but also how robustness and accuracy during development are generated.

To model the coupled system (Eqn 1), single cells were distributed spatially on a regular two-dimensional lattice, the dimensions of which are specified throughout the figures. The cells were coupled through paracrine signalling communication, by secreting and sensing the signalling molecule s within a certain distance. Four different communication ranges, R, were mainly considered, forming distinct network couplings (Fig. 1C): globally connected network (all-to-all communication, R=∞); locally connected network (cells communicate only with adjacent cells on the lattice, R=1, i.e. within a unit lattice distance); short-range connected network (cells communicate with cells within distance R=2, i.e. with adjacent cells and second-adjacent cells on the lattice); and probabilistic distance-based coupling [cells establish communication with other cells with probability ⁠, where d is the Euclidean cell-cell distance and R=1 in the default case]. No-flux or insulation boundary conditions on the lattices were used, hence cells near the boundaries had fewer communicating neighbours, as exemplified by the schematic for short-range coupling on a 2×4 grid in Fig. 1C. In all cases, is the external amount of signal perceived by cell i, depicted as the averaged secreted signal from its communicating neighbourhood N(i) (including itself) at every time instance, thus assuming an immediate mixing of s.

Each gene regulation term takes the renormalised Hill-type function form, xⁿ/(1+xⁿ) for activation or 1/(1+xⁿ) for repression, and is therefore sensitive to values of the order of 1 of the transcription factor input x. The corresponding maximal transcription rates are α_u, α_v, α_s and α_u,s, whereas β, γ, δ and η are the Hill coefficients. The reaction rates are globally scaled by λ. Values of α_u=2.3, α_v=3.5, α_s=2, α_u,s=1, β=γ=δ=η=2 and λ=50 were used throughout the study, unless noted otherwise. For the case of non-identical cells (Fig. S4A), the α_v parameter was uniformly varied between the cells in the range from −1% to 1% of its default value.

For the stochastic simulations (Fig. 1D, Fig. 2B, Fig. 3C, Fig. 4B-D, Fig. 5; Fig. S5B,C,E,F,H,I, Figs S8, S9, Movies 1-4), a stochastic differential equation model using Eqn 1 was constructed by adding a multiplicative noise σXΔW_t, where ΔW_t is the Wiener process term, i.e. a normally distributed random variable with zero mean and variance Δt, whereas X is a state variable (X∈{u₁, v₁, s₁, …, u_N, v_N, s_N}). The model was solved with Δt=0.01 using the Milstein method (Mil'shtein, 1974), by adding a second-order approximation term ⁠. In the cases when the noise intensity σ was not varied, it was set to 0.1.

To discriminate between the multilineage-primed- (mlp), u-positive (u+) or v-positive (v+) cell types for a given realisation, each marginal cell state vector (u_i, v_i, s_i) within the converged state of the system (IHSS or HSS) was individually categorised as one of three cell types and the three-term u+/v+/mlp ratio (proportions) in the population was subsequently calculated. The reference mlp state vector (u_mlp, v_mlp, s_mlp) was predetermined for a given parameter set, i.e. for a specific value of α_u, from the respective steady-state value of the 1-cell monostable system realisation (as in Fig. S2A), as the bifurcation analysis demonstrated that the mlp HSS for single-cell and cell-to-cell coupled systems is equivalent. Every marginal cell state (u_i, v_i, s_i) for of a deterministic realisation was categorised as mlp cell type when its value fell within 1% around the predetermined mlp state vector, whereas cell states whose u_i/v_i ratio was larger than the u_mlp/v_mlp ratio of the mlp state were assigned u+, and otherwise v+ cell types. Transient states during the stochastic realisations in Fig. 1D, Fig. 4D and Movies 1-4 were categorised as mlp type if they fell within 5% around the deterministic mlp state. End states of all stochastic realisations were allowed to converge to their deterministic attractor state in noise-free fashion before categorising (with 1% s.d. around the mlp state), as in Fig. 3C, Fig. 4B,C, Fig. 5; Figs S5B,C,E,F,H,I and Figs S8, S9.

Initial conditions for all cells and variables were independently randomly sampled from a normal distribution ⁠, typically around the corresponding α_u-specific mlp state vector (u_mlp, v_mlp, s_mlp) as the mean (μ_ics), and with σ_ics=0.1μ_ics or σ_ics=0.2μ_ics, as denoted in the respective figures. α_u-specific mlp states correspond to the steady-state values in the single-cell system (Fig. S2A). For Fig. 4A and Fig. S5A,D,G, σ_ics was increasing from 0 to 50% of the respective mlp value (μ_mlp). For Fig. 4B and Fig. S5B,E,H, μ_ics values were shifted stepwise along the line segment from a typical v + cell state (μ_v+) towards a typical u+ cell state (μ_u+), using μ_ics=kμ_u++(1−k)μ_v+ for ⁠, and with σ_ics=0.1μ_ics.

Finally, for the results from ten repetitions, in which the identified IHSSs with equivalent proportions were grouped, each proportion was plotted as a stacked sub-bar within a bar plot, the width of which corresponded to the fraction of occurrence of that proportion in the simulations. Proportions were plotted in ascending order of their u+ cell-type fractions. The results from 100 repetitions (Fig. 5; Fig. S9) were plotted equivalently, without vertical lines separating the sub-bars of the proportions.

In Fig. 2B, the different branches of the IHSS (and the proportions of cells in them) were estimated by analogy to Fig. 2A, but here using the number of cells as a bifurcation parameter. For this, exhaustive scanning was performed to locate the different fixed point attractors in the state space for each N. The scanning process involved 20 repetitive executions with different noise intensities (varying from 0 to 0.3). Each repetition consisted of 30 alternating cycles of stochastic execution (for exploration), followed by deterministic execution (for convergence to an attractor), when the reached state was tested for stability and subsequently recorded. For every distinctly detected attractor, the u+/v+/mlp proportion of cells was estimated, and attractors with equivalent proportions were grouped. The average u values were calculated separately for both u+ and v + cell types for each proportion, and plotted as branches (u+ and v + cells for IHSS, or mlp cells for HSS; colour-coding as in Fig. 2A), in analogy with the bifurcation diagrams in Fig. 2A and Fig. 3A. Chimera-like states were omitted from the diagram for simplicity.

The generation of lineage trees was performed using stochastic simulations in which the system doubles in size at regular time intervals, starting from a single-cell system (Fig. 1D), unless otherwise specified (Fig. 4D). At every cell division, the final state of the mother cell is passed on to the initial conditions of the daughter cells. Cell divisions occur along the horizontal and vertical axes on the grid alternately, sequentially yielding lattices of 1×1, 1×2, 2×2, 2×4, 4×4, 4×8, 8×8, etc., as demonstrated in Fig. S1. Cellular states were categorised in every time instance to plot the single-cell temporal evolutions in the lineage trees (Fig. 1D, Fig. 4D). Furthermore, cellular proportions in the system were estimated from those values, and their temporal evolution was shown in the panels above the lineage trees. The fraction of u+ cells was plotted to yield Fig. 3C, Fig. 5C, Fig. S8 (lower left plots) and Fig. S9C.

In the cell-type separation case (Fig. 4D), the states of the cells at the end of the fourth cycle (from Fig. 1D) were categorised and the differentiated cells were then separated: u+ cells were given as seeds to a new lineage tree (1×2 grid), whereas v + cells were seeds for a separate one (2×3 grid). Following this, multiple cell divisions were again performed and the cell proportions were estimated.

To characterise the spatial organisation of the system, lineage tree simulations with an extended duration of 13 cell cycles were performed, reaching a final size of a 64×64 grid. Both for deterministic and probabilistic distance-based coupling, communication with ranges of 1, 2, 3, 5 and 10 were considered. Additionally, global coupling was also used.

Spatial organisations were analysed at the end of each cell cycle, after the collective state reaches a steady state in a noise-free fashion. Fractions of u+/v+/mlp cells were quantified for the final spatial configuration (64×64 grid), and the results from 100 different repetitions were grouped as stacked bars. Moreover, the fractions of u+ cells were quantified for the end state of each cell cycle to track the development of the lineage tree (Fig. 3C; Fig. S8). The percentage of the numerical simulation realisations that showed symmetry breaking following a certain cell division cycle are also denoted on the plots (Fig. S8, lower left plots). Histograms depicting the proportion of numerical simulation realisations upon which a specific u+ cell-type fraction was reached in the final state of the lineage tree are also shown (Fig. S8, lower left plots).

To characterise the spatial clustering, the average u+ cluster radius was estimated from the spatial configurations by calculating for all u+ cells the fraction of u+ cells in the surround within increasing distance (Fig. S7). The cluster radius is therefore the distance by which the fraction of u+ cells drops to halfway between the fraction at zero distance (equal to 1), and the global fraction of u+ cells (dashed lines in Fig. S7), which is analogous to calculating half-life from exponential decay. The results from the realisations were plotted as single points in ascending order for each coupling range type in Fig. 5B and Fig. S9B. As with the fraction of u+ cells, the u+ cluster radius was also tracked over the development of the lineage tree, with a histogram being plotted for the final state as well. Both measures were also estimated when considering a population with a fixed size (64×64 grid), and starting from random initial conditions (μ_ics=μ_mlp, σ_ics=0.1μ_ics) (Fig. S9).

The lineage tree simulations shown in Movies 1-4 were performed with a reduced number of time points (2000) per cell cycle to yield shorter movies. Spatial configurations with stochastic categorisation of the cell types were saved for every tenth time point, and the saved frames were combined in a movie of 60 frames per second.

For Fig. S6A, deterministic simulations were performed for a fixed population size on a 4×8 grid, short-range coupling (other parameters as in Fig. 1D) and starting from random initial conditions (μ_ics=μ_mlp, σ_ics=0.1μ_ics). To model the decreased cell-cell communication strength, an inhibitory multiplicative term (1−s_inh,out), where s_inh,out∈{0, 0.05, 0.1, 0.15, 0.2, 0.25, 0.3}, was introduced to the first term in the equation for s in Eqn 1. This enabled the effective reduction of the production strength of s by a given percentage. The corresponding results from ten independent realisations for each condition are presented in Fig. S6A.

Similarly, for Fig. S6B, an inhibitory multiplicative term (1−s_inh,in), where s_inh,in∈{0, 0.0153, 0.0307, 0.046, 0.0613, 0.0767, 0.092}, was introduced to the second term in the equation for u in Eqn 1. Other settings were the same as for Fig. S6A.

The numerical bifurcation analysis was performed using the XPP/AUTO software (http://www.math.pitt.edu/∼bard/xpp/xpp.html). All simulations were performed using custom-made code in MATLAB (MATLAB and Statistics Toolbox Release R2020b, MathWorks).

We thank Philippe Bastiaens for numerous discussions, together with Peter Bieling and Dhruv Raina for critically reading the manuscript.

The peer review history is available online at https://dev.biologists.org/lookup/doi/10.1242/dev.197608.reviewer-comments.pdf