ABSTRACT

During mammalian development and homeostasis, cells often transition from a multilineage primed state to one of several differentiated cell types that are marked by the expression of mutually exclusive genetic markers. These observations have been classically explained by single-cell multistability as the dynamical basis of differentiation, where robust cell-type proportioning relies on pre-existing cell-to-cell differences. We propose a conceptually different dynamical mechanism in which cell types emerge and are maintained collectively by cell-cell communication as a novel inhomogeneous state of the coupled system. Differentiation can be triggered by cell number increase as the population grows in size, through organisation of the initial homogeneous population before the symmetry-breaking bifurcation point. Robust proportioning and reliable recovery of the differentiated cell types following a perturbation is an inherent feature of the inhomogeneous state that is collectively maintained. This dynamical mechanism is valid for systems with steady-state or oscillatory single-cell dynamics. Therefore, our results suggest that timing and subsequent differentiation in robust cell-type proportions can emerge from the cooperative behaviour of growing cell populations during development.

INTRODUCTION

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.

RESULTS

Timing and stable proportions of differentiated cell types are generated in growing populations

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.

Fig. 1.

Emergence of stable proportions of differentiated cell types in growing cell populations. (A) Schematic of the transition from multilineage primed cell types towards differentiated cell types during early embryogenesis. (B) Network topology of the cell-cell communication system (Eqn 1). (C) Different communication scenarios on a grid: short-range, global, local and distance-based probabilistic coupling. (D) Lineage tree depicting the transition from a multilineage primed state (mlp) to the u+/v+ states in a growing population with short-range communication, determined from a stochastic simulation. Green/red/blue: mlp/u+/v+ cells, respectively. Top: respective cell-type proportions. Parameters: αu=2.3, αv=3.5, αs=2, αu,s=1, β=γ=δ=η=2 and λ=50.

Fig. 1.

Emergence of stable proportions of differentiated cell types in growing cell populations. (A) Schematic of the transition from multilineage primed cell types towards differentiated cell types during early embryogenesis. (B) Network topology of the cell-cell communication system (Eqn 1). (C) Different communication scenarios on a grid: short-range, global, local and distance-based probabilistic coupling. (D) Lineage tree depicting the transition from a multilineage primed state (mlp) to the u+/v+ states in a growing population with short-range communication, determined from a stochastic simulation. Green/red/blue: mlp/u+/v+ cells, respectively. Top: respective cell-type proportions. Parameters: αu=2.3, αv=3.5, αs=2, αu,s=1, β=γ=δ=η=2 and λ=50.

To investigate under which conditions such transition occurs, we considered a generic cell-cell communication system, in which the gene regulatory circuits in single cells are coupled to each other by paracrine signalling molecules s (Fig. 1B). The intracellular circuit consists of two genes u and v that inhibit the transcription of each other, and u additionally positively regulates the secretion of s. In turn, the extracellular concentration of the signalling molecules, sext, regulates the circuit dynamics (inhibiting u production; Fig. 1B). As the communicating signals are secreted by the cells themselves, their concentration is a variable in the system. This effectively establishes a single joint dynamical system that is described by the following dimensionless equations:
formula
(1)
where is the single cell index and N is the number of cells. The vector (ui, vi, si) represents the state of each cell within the 3N-dimensional state of the coupled system (u1, v1, s1, …, uN, vN, sN). The external amount of signal perceived by cell i from its communicating neighbourhood N(i) (including itself) is defined as , i.e. the average amount of secreted signalling molecule at the given time. The parameters αu, αv, αu,s and αs represent the respective promoter strengths, whereas β, γ, η and δ are the Hill coefficients. Aiming to conceptualise the problem rather than to provide a quantitative description, the system was scaled to yield a minimal number of parameters: rate constants proportional to the promoter strengths were on the same order of magnitude, and the Hill coefficients were set to 2 (Materials and Methods).

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.

Collective state with heterogeneous cell types emerges with precise timing due to subcritical organisation

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 (u1=u2) as evident from the projection that falls on the diagonal in the u1u2 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.

Fig. 2.

Subcritical organisation before the PB enables timing of cellular differentiation. (A) Top: bifurcation diagram for N=2 coupled identical cells (inset), using u2 as a representative variable. Solid lines indicate HSS (black) and IHSS (purple). Dashed lines indicate unstable steady states. Dotted lines indicate organisations in parameter space, with corresponding stable mlp HSS (solid green circle, αu=2.3), u+/v+ IHSS (solid red/blue circles, αu=2.52) and unstable HSS (green cross, αu=2.52). Bottom: corresponding u1u2 state space organisation. Left: stable mlp HSS (u1=u2) for αu=2.3. Right: stable IHSS (u1>u2 and u2>u1) for αu=2.52. (B) Bifurcation-like diagram depicting the emergence of IHSS distributions with increasing number of cells N. Green boxes indicate mlp HSS. Red and blue stack bar markers represent distinct u+/v+ IHSS proportions. Lines connect distributions with the same proportions. (C) The symmetry-breaking population size threshold (NSB) at which the HSS becomes unstable in relation to αu. For probabilistic distance-based coupling, the median (n=200) is shown. Green-to-red/blue arrow denotes the symmetry-breaking mlp-to-u+/v+ transition at N=8 cells (orange star) for αu=2.3. In all panels, remaining parameters are the same as used in Fig. 1D.

Fig. 2.

Subcritical organisation before the PB enables timing of cellular differentiation. (A) Top: bifurcation diagram for N=2 coupled identical cells (inset), using u2 as a representative variable. Solid lines indicate HSS (black) and IHSS (purple). Dashed lines indicate unstable steady states. Dotted lines indicate organisations in parameter space, with corresponding stable mlp HSS (solid green circle, αu=2.3), u+/v+ IHSS (solid red/blue circles, αu=2.52) and unstable HSS (green cross, αu=2.52). Bottom: corresponding u1u2 state space organisation. Left: stable mlp HSS (u1=u2) for αu=2.3. Right: stable IHSS (u1>u2 and u2>u1) for αu=2.52. (B) Bifurcation-like diagram depicting the emergence of IHSS distributions with increasing number of cells N. Green boxes indicate mlp HSS. Red and blue stack bar markers represent distinct u+/v+ IHSS proportions. Lines connect distributions with the same proportions. (C) The symmetry-breaking population size threshold (NSB) at which the HSS becomes unstable in relation to αu. For probabilistic distance-based coupling, the median (n=200) is shown. Green-to-red/blue arrow denotes the symmetry-breaking mlp-to-u+/v+ transition at N=8 cells (orange star) for αu=2.3. In all panels, remaining parameters are the same as used in Fig. 1D.

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 (u1, v1, s1, u2, v2, s2), that describes a heterogeneous state. This state consists of a high u-expression level (u+) in one cell (u2) and a low u-expression level (v+) in the other cell (u1), as shown by the two-dimensional (u1, u2) 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 (u1<u2 and u1>u2) 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 u1u2 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 NSB 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 (NSB). 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.

Reliable proportioning of differentiated cell types is a dynamical consequence of the sequential ordering of IHSS solutions

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+/[Nk]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+/[Nk]v + and [k+1]u+/[Nk−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).

Fig. 3.

Reliable cell-type proportions are established through stable IHSS distributions. (A) Bifurcation analysis for N=4 globally/short-range-coupled cells (inset). Solid lines indicate stable HSS (black) and IHSS (blue, purple and red) branches. Three stable IHSS distributions with increasing u+/v + cell-type ratios appear sequentially (key). Dashed lines indicate unstable steady states. Dotted lines indicate organisations in parameter space. The grey shaded area indicates the parameter range of HSS/IHSS co-existence for subcritical organisation. Top: probabilities for visiting IHSS distributions (y-axis) at the respective αu values, estimated from 200 independent realisations. The initial conditions were randomly drawn from a normal distribution around the corresponding αu-specific mlp state (μics, σics=0.2μics). (B) u+/v+/mlp cell-type proportions for increasing αu values, when a short-range coupled population of N=4 (top), N=8 (middle) and N=32 cells (bottom) on a 2×2, 2×4 and 4×8 grid is considered, respectively. The width of each sub-bar within a bar reflects the fraction of occurrence of the respective u+/v+/mlp proportion in ten independent realisations. Initial conditions were the same as in A (top). Red arrows indicate the progression of the distribution of stable cell-type proportions with increase of system size (N=4,8,32 cells) for αu=2.3. (C) Temporal evolution of the fraction of u+ cells at each step of a lineage tree under short-range coupling and αu=2.3. Results from 100 independent numerical simulations are shown. Other parameters are the same as in Fig. 1D.

Fig. 3.

Reliable cell-type proportions are established through stable IHSS distributions. (A) Bifurcation analysis for N=4 globally/short-range-coupled cells (inset). Solid lines indicate stable HSS (black) and IHSS (blue, purple and red) branches. Three stable IHSS distributions with increasing u+/v + cell-type ratios appear sequentially (key). Dashed lines indicate unstable steady states. Dotted lines indicate organisations in parameter space. The grey shaded area indicates the parameter range of HSS/IHSS co-existence for subcritical organisation. Top: probabilities for visiting IHSS distributions (y-axis) at the respective αu values, estimated from 200 independent realisations. The initial conditions were randomly drawn from a normal distribution around the corresponding αu-specific mlp state (μics, σics=0.2μics). (B) u+/v+/mlp cell-type proportions for increasing αu values, when a short-range coupled population of N=4 (top), N=8 (middle) and N=32 cells (bottom) on a 2×2, 2×4 and 4×8 grid is considered, respectively. The width of each sub-bar within a bar reflects the fraction of occurrence of the respective u+/v+/mlp proportion in ten independent realisations. Initial conditions were the same as in A (top). Red arrows indicate the progression of the distribution of stable cell-type proportions with increase of system size (N=4,8,32 cells) for αu=2.3. (C) Temporal evolution of the fraction of u+ cells at each step of a lineage tree under short-range coupling and αu=2.3. Results from 100 independent numerical simulations are shown. Other parameters are the same as in Fig. 1D.

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.

The IHSS distributions confer robust cell-type proportioning and mediate its recovery from perturbations

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.

Fig. 4.

Robustness in cell-type proportioning and plasticity in its re-establishment upon perturbations. (A-C) Cell-type proportions upon increase in the s.d. around the mlp initial conditions (as a fraction from the mean) (A); a shift in the mean of the distribution (μics) value from high v-expression (μv+) to high u-expression cell state (μu+), with μics=u++(1−k)μv+ for and σics=0.1μics (B); and an increase in the noise intensity (Materials and Methods) (C). Top: respective perturbations; 1-sigma ellipses of the distributions are depicted. Bottom: proportions from ten independent numerical realisations per condition are shown, estimated for N=32 short-range coupled cells on a 4×8 grid and organisation before the PB, αu=2.3. (D) Lineage trees generated from homogeneous sub-populations of differentiated cells, separated at the fourth step of the lineage tree in Fig. 1D (N=8 cells). Left: lineage tree seeded from N=2 cells that previously had adopted a high u-expression state (u+). Right: lineage tree from the N=6 cells (2×3 grid) that before separation adopted the high v-expression state (v+). Upper panels show the respective cell-type proportions. Parameters are the same as in Fig. 1D.

Fig. 4.

Robustness in cell-type proportioning and plasticity in its re-establishment upon perturbations. (A-C) Cell-type proportions upon increase in the s.d. around the mlp initial conditions (as a fraction from the mean) (A); a shift in the mean of the distribution (μics) value from high v-expression (μv+) to high u-expression cell state (μu+), with μics=u++(1−k)μv+ for and σics=0.1μics (B); and an increase in the noise intensity (Materials and Methods) (C). Top: respective perturbations; 1-sigma ellipses of the distributions are depicted. Bottom: proportions from ten independent numerical realisations per condition are shown, estimated for N=32 short-range coupled cells on a 4×8 grid and organisation before the PB, αu=2.3. (D) Lineage trees generated from homogeneous sub-populations of differentiated cells, separated at the fourth step of the lineage tree in Fig. 1D (N=8 cells). Left: lineage tree seeded from N=2 cells that previously had adopted a high u-expression state (u+). Right: lineage tree from the N=6 cells (2×3 grid) that before separation adopted the high v-expression state (v+). Upper panels show the respective cell-type proportions. Parameters are the same as in Fig. 1D.

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−sinh,out) and (1−sinh,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 sinh,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 sinh,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.

IHSS leads to reproducible spatial patterns in growing cell populations

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).

Fig. 5.

Robust cell-type proportions are correlated with reproducible spatial patterns. (A) Exemplary spatial patterns for different coupling ranges: local, short-range, 3-range and 10-range coupling. (B) u+ cluster radius estimated at the final stable state of the lineage tree (64×64 size) for deterministic (left), probabilistic (middle) and global all-to-all coupling (right). Results from 100 independent stochastic realisations per coupling type are plotted in ascending order. (C) Corresponding u+/v+/mlp cell-type proportions. Bars are analogous to those in Fig. 3B. (D) Exemplary spatial patterns for range 1, 3 and 10 probabilistic coupling and global all-to-all coupling. Parameters are the same as in Fig. 1D.

Fig. 5.

Robust cell-type proportions are correlated with reproducible spatial patterns. (A) Exemplary spatial patterns for different coupling ranges: local, short-range, 3-range and 10-range coupling. (B) u+ cluster radius estimated at the final stable state of the lineage tree (64×64 size) for deterministic (left), probabilistic (middle) and global all-to-all coupling (right). Results from 100 independent stochastic realisations per coupling type are plotted in ascending order. (C) Corresponding u+/v+/mlp cell-type proportions. Bars are analogous to those in Fig. 3B. (D) Exemplary spatial patterns for range 1, 3 and 10 probabilistic coupling and global all-to-all coupling. Parameters are the same as in Fig. 1D.

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.

DISCUSSION

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.

MATERIALS AND METHODS

Modelling a generic cell-cell communication system

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, xn/(1+xn) for activation or 1/(1+xn) 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ΔWt, where ΔWt 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∈{u1, v1, s1, …, uN, vN, sN}). 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 (ui, vi, si) 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 (umlp, vmlp, smlp) 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 (ui, vi, si) 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 ui/vi ratio was larger than the umlp/vmlp 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 (umlp, vmlp, smlp) 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=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.

Estimating IHSS distributions as a function of the number of cells

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.

Lineage tree generation

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.

Analysis of spatial patterns

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.

Decreasing cell-cell communication strength

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−sinh,out), where sinh,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−sinh,in), where sinh,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.

Additional cell-cell communication systems exhibiting population-based symmetry breaking

To demonstrate the generality of the population-based symmetry-breaking mechanism via a PB, we also tested several minimal network topologies with effective negative feedback coupling on u via s (Fig. S3A-C). Their dynamics were described using the following system of equations:
formula
(2)
For Fig. S3A, αu,s=0, αs=1, αv,s=1, αs,v=0 and αv=2.75. For Fig. S3B, αu,s=0.5, αs=0, αv,s=0, αs,v=3 and αv=3. For Fig. S3C, αu,s=0, αs=0, αv,s=1, αs,v=2 and αv=3.

Paradigmatic model mimicking the vertebrate neurogenesis process

It has been previously demonstrated that the presence of time delays in models of lateral inhibition can result in significant oscillatory transients before patterned steady states are reached. The impact of local feedback loops in a model of lateral inhibition based on the Notch signalling pathway, elucidating the roles of intracellular and intercellular delays in controlling the overall system behaviour, has also been proposed (Momiji and Monk, 2009). Here, we aimed to understand whether population-based PB can provide the dynamical background behind the observed symmetry-breaking phenomenon. As our aim was to demonstrate the validity of this concept, we omitted the molecular details of the Notch pathway example and modelled a generic case in which the gene expression dynamics in each cell were characterised by oscillatory behaviour, whereas intercellular communication between the N cells was realised in a global manner, for simplicity (N=2, Fig. S3D). The dynamics of the system were therefore described as:
formula
(3)
Here, p and q are two genes that mutually inhibit the expression of each other, r controls the production of the signalling molecule, the extracellular concentration of which is denoted as rext. This system has been demonstrated to exhibit synchronised oscillations in a population of communicating cells (Kuznetsov et al., 2004; Koseska et al., 2007). The parameters are as follows: αq=5, αp,r=1, αr=4, β=2, γ=2, η=2, ε=0.01, d=0.008 and de=1.

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).

Acknowledgements

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

Footnotes

Author contributions

Conceptualization: A.K.; Methodology: A.S., A.K.; Software: A.S.; Validation: A.S., C.S., A.K.; Formal analysis: A.S.; Investigation: A.S., A.K.; Writing - original draft: A.K.; Writing - review & editing: A.S., C.S., A.K.; Visualization: A.S.; Supervision: A.K.

Funding

The authors acknowledge funding from the Max-Planck-Gesellschaft.

Data availability

All data and code used in this study are available at www.github.com/astanoev/SymmetryBreaking.

Peer review history

References

Abrams
,
D. M.
and
Strogatz
,
S. H.
(
2004
).
Chimera states for coupled oscillators
.
Phys. Rev. Lett.
93
,
174102
.
Andrecut
,
M.
,
Halley
,
J. D.
,
Winkler
,
D. A.
and
Huang
,
S.
(
2011
).
A general model for binary cell fate decision gene circuits with degeneracy: indeterminacy and switch behavior in the absence of cooperativity
.
PLoS ONE
6
,
e19358
.
Bessonnard
,
S.
,
De Mot
,
L.
,
Gonze
,
D.
,
Barriol
,
M.
,
Dennis
,
C.
,
Goldbeter
,
A.
,
Dupont
,
G.
and
Chazaud
,
C.
(
2014
).
Gata6, Nanog and Erk signaling control cell fate in the inner cell mass through a tristable regulatory network
.
Development
141
,
3637
-
3648
.
Cherry
,
J. L.
and
Adler
,
F. R.
(
2000
).
How to make a biological switch
.
J. Theor. Biol.
203
,
117
-
133
.
Chickarmane
,
V.
and
Peterson
,
C.
(
2008
).
A computational model for understanding stem cell, trophectoderm and endoderm lineage determination
.
PLoS ONE
3
,
e3478
.
Collier
,
J. R.
,
Monk
,
N. A. M.
,
Maini
,
P. K.
and
Lewis
,
J. H.
(
1996
).
Pattern formation by lateral inhibition with feedback: a mathematical model of delta-notch intercellular signalling
.
J. Theor. Biol.
183
,
429
-
446
.
Danino
,
T.
,
Mondragón-Palomino
,
O.
,
Tsmiring
,
L.
and
Hasty
,
J.
(
2010
).
A synchronized quorum of genetic clocks
.
Nature
463
,
326
-
330
.
De Caluwé
,
J.
,
Tosenberger
,
A.
,
Gonze
,
D.
and
Dupont
,
G.
(
2019
).
Signalling-modulated gene regulatory networks in early mammalian development
.
J. Theor. Biol.
463
,
56
-
66
.
De Mot
,
L.
,
Gonze
,
D.
,
Bessonnard
,
S.
,
Chazaud
,
C.
,
Goldbeter
,
A.
and
Dupont
,
G.
(
2016
).
Cell fate specification based on tristability in the inner cell mass of mouse blastocysts
.
Biophys. J.
110
,
710
-
722
.
Economou
,
A. D.
,
Ohazama
,
A.
,
Porntaveetus
,
T.
,
Sharpe
,
P. T.
,
Kondo
,
S.
,
Basson
,
M. A.
,
Gritli-Linde
,
A.
,
Cobourne
,
M. T.
and
Green
,
J. B. A.
(
2012
).
Periodic stripe formation by a Turing mechanism operating at growth zones in the mammalian palate
.
Nat. Genet.
44
,
348
.
Enver
,
T.
,
Pera
,
M.
,
Peterson
,
C.
and
Andrews
,
P. W.
(
2009
).
Stem cell states, fates, and the rules of attraction
.
Cell Stem Cell
4
,
387
-
397
.
Ferrell
,
J. E.
(
2012
).
Bistability, bifurcations, and Waddington's epigenetic landscape
.
Curr. Biol.
22
,
R458
-
R466
.
Fischer
,
S. C.
,
Corujo-Simon
,
E.
,
Lilao-Garzon
,
J.
,
Stelzer
,
E. H.
and
Muñoz-Descalzo
,
S.
(
2020
).
The transition from local to global patterns governs the differentiation of mouse blastocysts
.
PLoS ONE
15
,
e0233030
.
Formosa-Jordan
,
P.
and
Ibañes
,
M.
(
2014
).
Competition in notch signaling with cis enriches cell fate decisions
.
PLoS ONE
9
,
e95744
.
Garcia-Ojalvo
,
J.
,
Elowitz
,
M. B.
and
Strogatz
,
S. H.
(
2004
).
Modeling a synthetic multicellular clock: repressilators coupled by quorum sensing
.
Proc. Natl Acad. Sci. USA
101
,
10955
-
10960
.
Gupta
,
P. B.
,
Fillmore
,
C. M.
,
Jiang
,
G.
,
Shapira
,
S. D.
,
Tao
,
K.
,
Kuperwasser
,
C.
and
Lander
,
E. S.
(
2011
).
Stochastic state transitions give rise to phenotypic equilibrium in populations of cancer cells
.
Cell
146
,
633
-
644
.
Hart
,
Y.
,
Reich-Zeliger
,
S.
,
Antebi
,
Y. E.
,
Zaretsky
,
I.
,
Mayo
,
A. E.
,
Alon
,
U.
and
Friedman
,
N.
(
2014
).
Paradoxical signaling by a secreted molecule leads to homeostasis of cell levels
.
Cell
158
,
1022
-
1032
.
Hatakeyama
,
J.
,
Bessho
,
Y.
,
Katoh
,
K.
,
Ookawara
,
S.
,
Fujioka
,
M.
,
Guillemot
,
F.
and
Kageyama
,
R.
(
2004
).
Hes genes regulate size, shape and histogenesis of the nervous system by control of the timing of neural stem cell differentiation
.
Development
131
,
5539
-
5550
.
Huang
,
S.
,
Guo
,
Y.-P.
,
May
,
G.
and
Enver
,
T.
(
2007
).
Bifurcation dynamics in lineage-commitment in bipotent progenitor cells
.
Dev. Biol.
305
,
695
-
713
.
Jia
,
D.
,
Jolly
,
M. K.
,
Harrison
,
W.
,
Boareto
,
M.
,
Ben-Jacob
,
E.
and
Levine
,
H.
(
2017
).
Operating principles of tristable circuits regulating cellular differentiation
.
Phys. Biol.
14
,
035007
.
Kauffman
,
S. A.
(
1969
).
Metabolic stability and epigenesis in randomly constructed genetic nets
.
J. Theor. Biol.
22
,
437
-
467
.
Koseska
,
A.
and
Bastiaens
,
P. I. H.
(
2017
).
Cell signaling as a cognitive process
.
EMBO J.
36
,
568
-
582
.
Koseska
,
A.
,
Volkov
,
E.
,
Zaikin
,
A.
and
Kurths
,
J.
(
2007
).
Inherent multistability in arrays of autoinducer coupled genetic oscillators
.
Phys. Rev. E
75
,
031916
.
Koseska
,
A.
,
Volkov
,
E.
and
Kurths
,
J.
(
2009
).
Detuning-dependent dominance of oscillation death in globally coupled synthetic genetic oscillators
.
EPL
85
,
28002
.
Koseska
,
A.
,
Ullner
,
E.
,
Volkov
,
E.
,
Kurths
,
J.
and
García-Ojalvo
,
J.
(
2010
).
Cooperative differentiation through clustering in multicellular populations
.
J. Theor. Biol.
263
,
189
-
202
.
Koseska
,
A.
,
Volkov
,
E.
and
Kurths
,
J.
(
2013
).
Transition from amplitude to oscillation death via Turing bifurcation
.
Phys. Rev. Lett.
111
,
024103
.
Kuramoto
,
Y.
and
Battogtokh
,
D.
(
2002
).
Coexistence of coherence and incoherence in nonlocally coupled phase oscillators
.
Nonlinear Phenom. Complex Syst.
5
,
380
-
385
.
Kuznetsov
,
A.
,
Kærn
,
M.
and
Kopell
,
N.
(
2004
).
Synchrony in a population of hysteresis-based genetic oscillators
.
SIAM J. Appl. Math.
65
,
392
-
425
.
Lorthongpanich
,
C.
,
Doris
,
T. P. Y.
,
Limviphuvadh
,
V.
,
Knowles
,
B. B.
and
Solter
,
D.
(
2012
).
Developmental fate and lineage commitment of singled mouse blastomeres
.
Development
139
,
3722
-
3731
.
Martinez Arias
,
A.
,
Nichols
,
J.
and
Schroeter
,
C.
(
2013
).
A molecular basis for developmental plasticity in early mammalian embryos
.
Development
140
,
3499
.
Mathew
,
B.
,
Muñoz-Descalzo
,
S.
,
Corujo-Simon
,
E.
,
Schröter
,
C.
,
Stelzer
,
E.
and
Fischer
,
S. C.
(
2019
).
Mouse icm organoids reveal three-dimensional cell fate clustering
.
Biophys. J.
116
,
127
-
141
.
McMillen
,
D.
,
Kopell
,
N.
,
Hasty
,
J.
and
Collins
,
J. J.
(
2002
).
Synchronizing genetic relaxation oscillators by intercell signaling
.
Proc. Natl Acad. Sci. USA
99
,
679
-
684
.
Mil'shtein
,
G. N.
(
1974
).
Approximate integration of stochastic differential equations
.
Teor. Veroyatn. Primen.
19
,
583
-
588
.
Momiji
,
H.
and
Monk
,
N. A. M.
(
2009
).
Oscillatory Notch-pathway activity in a delay model of neuronal differentiation
.
Phys. Rev. E
80
,
021930
.
Monk
,
N. A. M.
(
1997
).
Cell communities and robustness in development
.
Bull. Mat. Biol.
59
,
1183
-
1189
.
Nichols
,
J.
,
Silva
,
J.
,
Roode
,
M.
and
Smith
,
A.
(
2009
).
Suppression of Erk signalling promotes ground state pluripotency in the mouse embryo
.
Development
136
,
3215
-
3222
.
Niwa
,
H.
,
Toyooka
,
Y.
,
Shimosato
,
D.
,
Strumpf
,
D.
,
Takahashi
,
K.
,
Yagi
,
R.
and
Rossant
,
J.
(
2005
).
Interaction between Oct3/4 and Cdx2 determines trophectoderm differentiation
.
Cell
123
,
917
-
929
.
O'Dea
,
R. D.
and
King
,
J. R.
(
2012
).
Continuum limits of pattern formation in hexagonal-cell monolayers
.
J. Math. Biol.
64
,
579
-
610
.
Raina
,
D.
,
Stanoev
,
A.
,
Bahadori
,
A.
,
Protzek
,
M.
,
Koseska
,
A.
and
Schröter
,
C.
(
2020
).
Cell-cell communication through FGF4 generates and maintains robust proportions of differentiated cell fates in embryonic stem cells
.
bioRxiv
, p.
2020.02.14.949701
.
Raspopovic
,
J.
,
Marcon
,
L.
,
Russo
,
L.
and
Sharpe
,
J.
(
2014
).
Digit patterning is controlled by a Bmp-Sox9-Wnt Turing network modulated by morphogen gradients
.
Science
345
,
566
-
570
.
Saiz
,
N.
,
Williams
,
K. M.
,
Seshan
,
V. E.
and
Hadjantonakis
,
A.-K.
(
2016
).
Asynchronous fate decisions by single cells collectively ensure consistent lineage composition in the mouse blastocyst
.
Nat. Commun.
7
,
13463
.
Saiz
,
N.
,
Mora-Bitria
,
L.
,
Rahman
,
S.
,
George
,
H.
,
Herder
,
J.
,
Garcia-Ojalvo
,
J.
and
Hadjantonakis
,
A.-K.
(
2020
).
Growth-factor-mediated coupling between lineage size and cell fate choice underlies robustness of mammalian development
.
eLife
9
,
e56079
.
Schröter
,
C.
,
Rué
,
P.
,
Mackenzie
,
J. P.
and
Martinez Arias
,
A.
(
2015
).
FGF/MAPK signaling sets the switching threshold of a bistable circuit controlling cell fate decisions in embryonic stem cells
.
Development
142
,
4205
-
4216
.
Simon
,
C. S.
,
Hadjantonakis
,
A.-K.
and
Schröter
,
C.
(
2018
).
Making lineage decisions with biological noise: Lessons from the early mouse embryo
.
Wiley Interdiscipl. Rev. Dev. Biol.
7
,
e319
.
Snoussi
,
E. H.
(
1998
).
Necessary conditions for multistationarity and stable periodicity
.
J. Biol. Syst.
6
,
3
-
9
.
Soldatov
,
R.
,
Kaucka
,
M.
,
Kastriti
,
M. E.
,
Petersen
,
J.
,
Chontorotzea
,
T.
,
Englmaier
,
L.
,
Akkuratova
,
N.
,
Yang
,
Y.
,
Häring
,
M.
,
Dyachuk
,
V.
et al. 
(
2019
).
Spatiotemporal structure of cell fate decisions in murine neural crest
.
Science
364
,
eaas9536
.
Sprinzak
,
D.
,
Lakhanpal
,
A.
,
LeBon
,
L.
,
Santat
,
L. A.
,
Fontes
,
M. E.
,
Anderson
,
G. A.
,
Garcia-Ojalvo
,
J.
and
Elowitz
,
M. B.
(
2010
).
Cis-interactions between notch and delta generate mutually exclusive signalling states
.
Nature
465
,
86
-
90
.
Suzuki
,
N.
,
Furusawa
,
C.
and
Kaneko
,
K.
(
2011
).
Oscillatory protein expression dynamics endows stem cells with robust differentiation potential
.
PLoS ONE
6
,
e27232
.
Taga
,
M. E.
and
Bassler
,
B. L.
(
2003
).
Chemical communication among bacteria
.
Proc. Natl Acad. Sci. USA
100
,
14549
-
14554
.
Teomy
,
E.
,
Kessler
,
D. A.
and
Levine
,
H.
(
2019
).
Ordered hexagonal patterns via notch-delta signaling
.
arXiv preprint arXiv
,
1902.04917
.
Thomas
,
R.
(
1981
).
On the relation between the logical structure of systems and their ability to generate multiple steady states or sustained oscillations
. In
Numerical Methods in The Study Of Critical Phenomena
,
Vol. 9 (ed. J. Della Dora, J. Demongeot and B. Lacolle)
, pp.
180
-
193
.
Berlin: Springer
.
Turing
,
A. M.
(
1952
).
The chemical basis of morphogenesis
.
Philos. Trans. R. Soc. Lond. Ser. B
237
,
37
-
72
.
Ullner
,
E.
,
Zaikin
,
A.
,
Volkov
,
E. I.
and
García-Ojalvo
,
J.
(
2007
).
Multistability and clustering in a population of synthetic genetic oscillators via phase-repulsive cell-to-cell communication
.
Phys. Rev. Lett.
99
,
148103
.
Ullner
,
E.
,
Koseska
,
A.
,
Kurths
,
J.
,
Volkov
,
E.
,
Kantz
,
H.
and
García-Ojalvo
,
J.
(
2008
).
Multistability of synthetic genetic networks with repressive cell-to-cell communication
.
Phys. Rev. E
78
,
031904
.
van Kekem
,
D. L.
and
Sterk
,
A. E.
(
2019
).
Symmetries in the Lorenz-96 model
.
Int. J. Bifurcation Chaos
29
,
1950008
.
Wang
,
J.
,
Zhang
,
K.
,
Xu
,
L.
and
Wang
,
E.
(
2011
).
Quantifying the Waddington landscape and biological paths for development and differentiation
.
Proc. Natl Acad. Sci. USA
108
,
8257
-
8262
.
Yamanaka
,
Y.
,
Lanner
,
F.
and
Rossant
,
J.
(
2010
).
FGF signal-dependent segregation of primitive endoderm and epiblast in the mouse blastocyst
.
Development
137
,
715
-
724
.
You
,
L.
,
Cox
,
R. S.
, III
,
Weiss
,
R.
and
Arnold
,
F. H.
(
2004
).
Programmed population control by cell-cell communication and regulated killing
.
Nature
428
,
868
-
871
.
Youk
,
H.
and
Lim
,
W. A.
(
2014
).
Secreting and sensing the same molecule allows cells to achieve versatile social behaviors
.
Science
343
,
6171
.
Zakharova
,
A.
,
Schneider
,
I.
,
Kyrychko
,
Y. N.
,
Blyuss
,
K. B.
,
Koseska
,
A.
,
Fiedler
,
B.
and
Schöll
,
E.
(
2013
).
Time delay control of symmetry-breaking primary and secondary oscillation death
.
EPL (Europhysics Letters)
104
,
50004
.
Zhang
,
H. T.
and
Hiiragi
,
T.
(
2018
).
Symmetry breaking in the mammalian embryo
.
Annu. Rev. Cell Dev. Biol.
34
,
405
-
426
.

Competing interests

The authors declare no competing or financial interests.

Supplementary information