Statistically derived geometrical landscapes capture principles of decision-making dynamics during cell fate transitions

1Cell fate decisions in developing tissues involve gene regulatory networks comprising multiple genes, many molecular components, and elaborate signaling dynamics. Despite the complexity, the outcome of cellular decisions is relatively simple: cells transition between a limited set of discrete cell fates, each defined by a distinct gene expression profile. These transitions occur in a characteristic sequence regulated by extrinsic signaling. While quantitative models that describe signaling pathways and gene regulatory networks in great detail have been used to investigate cell differentiation and decision making, these suffer from a plethora of parameters and their behavior is difficult to predict without case-by-case simulation. Hence, quantitative methods based on experimental data that represent and permit analysis of developmental processes at the scale of cell fate decisions would provide insight into the underlying principles and allow quantitative and testable predictions.

2A popular and intuitive metaphor for the process of developmental decision making is the Waddington landscape, in which the differentiation trajectory of a cell is conceived as a ball rolling down a landscape of branching valleys, representing specific cell fates. This can be mathematically formalized, with only minor changes, using dynamical systems theory. In this formulation, the relevant dynamical systems are gradient-like and the system’s trajectories, which represent the developmental path of a cell, move downhill in this landscape. Thus, Waddington’s valleys correspond to the attractors of the system and sit at the minima of the landscape. Moreover, variation in the parameters of the dynamical system, caused by changes in the signals the cell receives, alter the landscape and give rise to bifurcations that destroy or create attractors. The transition to a new fate is signified by a cell entering a new basin of attraction, caused either by a signal-induced bifurcation or by a stochastic fluctuation resulting in a cell jumping from one attractor’s pull to another. In both cases, the route from the old to the new cell state is defined by a saddle point in the landscape. This approach enables a rigorous link between the dynamical systems underlying gene regulatory networks and Waddington’s landscapes.

3There are several advantages to a geometrical viewpoint. First, since motion is always downhill, it gives a hierarchical structure to the dynamics and an intuitive understanding of eventual fates. Second, a body of dynamical systems theory indicates that the qualitative structure of the landscape can be described by a relatively small corpus of universal normal forms. In particular, Rene Thom’s catastrophe theory provides a powerful classification scheme of the relevant bifurcations of such systems that is facilitated by the existence of the gradient-like structure. This theory suggests that although a system’s dimension might be large, the bifurcations can be described by low-dimensional systems. For example, in developmental systems, it is common that differentiating cells transition from a progenitor state to one of two progeny fates. Such a decision can be an all-or-nothing one in which all cells make the same choice, or it can be one in which some cells make one choice and some the other one, enabling the allocation of a cell population to both fates. We introduce two 3-attractor landscapes, the binary choice (all-or-nothing) and the binary flip (allocation), that are the simplest archetypes underlying these two decision types.

4To transform a developmental process from a metaphorical landscape description into a geometric model that allows quantitative and qualitative experimental predictions, quantitative data are needed together with methods to connect experimental measurements to a parameterized dynamical landscape. A framework based on catastrophe theory and approximate Bayesian computation (ABC) has been developed and applied to experiments measuring the final outcome of a developmental process. However, to test the full power of this approach and assess its general usefulness, quantitative measurements and perturbations during the differentiation process are necessary. To this end, we took advantage of the differentiation of neural and mesodermal cells from pluripotent progenitors using mouse embryonic stem cells (ESCs) exposed to different combinations and durations of signaling factors (Figure 1A;;; ). ABC methods for parameter estimation based on matching summary statistics of temporally ordered data to those produced from simulations of candidate models naturally lend itself to our data and we use the proportion of cells in cell states at different time points as the summary statistic.

5Initiating the differentiation of ESCs, by withdrawing media supplements that sustain pluripotency, results in cells adopting a state that mimics the post-implantation epiblast. In the absence of further signals, these cells differentiate into neural progenitors with a molecular identity of the anterior nervous system. However, exposure to WNT signaling at the epiblast stage blocks ESCs from adopting an anterior neural fate and instead cells acquire a caudal epiblast identity. This cell type, which in the embryo fuels axis elongation and the formation of trunk tissue, is responsible for generating the progenitors of the spinal cord and paraxial mesoderm. Similarly, in vitro, caudal epiblast-like cells differentiate into spinal cord and mesoderm progenitors, with longer durations of WNT signaling resulting in a higher proportion of mesodermal cells at the expense of neural cells. Each of the cell states in this differentiation process is recognizable from well-defined gene expression profiles that can be assayed using representative marker proteins. Hence, the differentiation of ES cells into neural and mesodermal derivatives offers a well-characterized system in which to develop and test a dynamical landscape model of developing tissues and cellular decision making.

6We used flow cytometry to measure the expression levels of representative marker proteins in differentiating ESCs exposed to different signaling dynamics. These data provided an informative low-dimensional representation of cell states from which we developed a principled statistical approach to identify the attractors and the geometric form of the dynamical landscape. Using the ABC approach based on summary statistics, we parameterized and refined the landscape and its changes with extensive training and validation datasets and tested it against specific experiments prompted by predictions from the model. The result was a quantitative model of the differentiation of neural and mesodermal tissue from pluripotent progenitors. The data suggested two distinct decision-making mechanisms—a “binary choice” as a one-or-other decision that commits cells to either anterior neural or caudal epiblast and a “binary flip” that simultaneously allocates cells to posterior neural and mesodermal fates. These two mechanisms correspond to two fundamentally different landscape geometries with distinct qualitative and quantitative features. We discuss the biological relevance of these different mechanisms and suggest that they represent two dynamical archetypes that play a general role in differentiation dynamics. Taken together, the approach is broadly applicable for the quantitative analysis of cell fate dynamics and determining the design logic of developmental decisions.

7To develop a quantitative landscape model of cell fate decision making, we took advantage of an in vitro system in which pluripotent mouse ESCs are directed to differentiate to distinct neural and mesodermal fates in response to defined signaling dynamics (Figure 1A). ESCs removed from pluripotency conditions and grown in basal media containing FGF2 (FGF) and an inhibitor of endogenous WNT secretion (LGK974) adopt a post-implantation epiblast-like (EPI) identity by day (D)2 of culture (Figure 1Ba) and, subsequently, following removal of exogenous FGF, adopt an OTX2+/ SOX1+ molecular identity characteristic of anterior neural progenitors (AN) (Figure 1Bb). In contrast, if WNT signaling is activated at D2 by the GSK3β antagonist CHIRON99021 (CHIR) concurrently with FGF signaling, EPI cells no longer differentiate to AN but instead, acquire a CDX2+ caudal epiblast (CE) identity at D3 (Figure 1Bc). These differentiate to SOX1+ posterior neural progenitors (PN) in response to withdrawal of FGF and WNT signal activation (Figure 1Be), while sustained WNT signaling drives the progressive differentiation of CE to TBX6+ early paraxial mesoderm (EM) and FOXC2+ late paraxial mesoderm (LM) identities (Figures 1Bf and 1Ca).

8To develop a quantitative description of the differentiation process, we developed a flow-cytometry assay to measure, simultaneously, the expression of multiple marker proteins in individual cells. This enabled us to generate time-course data at single-cell resolution in a manner that was sufficiently scalable to explore the differentiation outcomes resulting from combinatorial modulation of FGF and WNT signaling dynamics. By inhibiting endogenous WNT secretion in all culture conditions with LGK974, we ensured that WNT signaling dynamics were entirely dependent upon the timing and concentration of CHIR addition. Moreover, as WNT signaling is known to induce FGF ligand production in CE cells, we tightly controlled FGF signaling dynamics through a combination of exogenous FGF withdrawal and inhibition of downstream FGF signaling with the small-molecule PD0325901 (PD). In total, we obtained 149 datasets in 4 experimental series, comprising 23 different signaling conditions assayed at 7 timepoints each (supplementary Table S1).

9Based on previous transcriptome analysis, we first defined a minimal set of markers that are sufficient to classify cells into distinct progenitor types. We employed the anterior marker OTX2 to distinguish between AN and PN progenitors, both of which express the neural progenitor markers SOX2 and SOX1. As expected, OTX2 expression was permanently extinguished in SOX2+ SOX1+ PN progenitors following combined activation of FGF and WNT signaling between D2 and D3, whereas its expression remained high in SOX1+ AN progenitors (Figure 1Cb). Moreover, the level of SOX1 was higher in PN cells, whereas SOX2 was higher in AN (Figure 1Cb). We reasoned that SOX2 and SOX1 levels might be sufficient for classification of AN and PN progenitors in the absence of OTX2 measurement. This was subsequently confirmed experimentally (see below). In addition, we confirmed that exposure to WNT signaling induced CE markers T/BRA and CDX2 in many cells at D3/D3.5 (Figure S1C). Continued exposure to WNT signaling resulted in the loss of both T/BRA and SOX2 expression from differentiating CE progenitors at D4 and a high proportion of cells expressed TBX6, characteristic of EM identity (Figure S1C), which itself began to decrease at D5 as cells acquired FOXC2+ LM identity (Figure 1Ca). As SOX2 expression is extinguished as cells commit to paraxial mesoderm identity (; Figure S1B) we reasoned that either the presence of TBX6 and absence of SOX2 expression, or the absence of both TBX6 and SOX2 was sufficient to classify mesoderm progenitors (M).

10In order to quantify the differentiation process, we hypothesized that each cell state represented an attractor of a dynamical landscape and we set out to develop a principled procedure to allocate cells to an identity. To this end, we used an algorithm based on fitting a Gaussian mixture model (GMM) to the multi-dimensional flow-cytometry data (fitgmdist in MATLAB). The GMM consists of a set of weighted multivariate normal distributions. We regarded each of these distributions as defining a cluster. A cell belonged to a cluster if its probability for the corresponding distribution was greater than a predefined threshold. Close to an attractor, the state variables (levels of protein expression in this case) are expected to have an approximately multivariate normal distribution since the dynamics near attractors are expected to be approximately linear.

11To mitigate batch effects, we included a reference set consisting of the three signaling regimes comprising no CHIR, CHIR 2–3, and CHIR 2–5 (Figure 1A) in each experimental series (initial, test, and prediction; see Table S1), and used the corresponding pooled data from these three conditions to fit and define the GMM probability distributions at each timepoint (Figure 2A, Methods S1 Section 2). To choose the number of distributions in the GMM, m, we increased it sequentially from 1 until the scalar distributions for each cluster were clearly unimodal in each of the dimensions (Figure 2B). Cells from all other samples within the experimental series were then classified by assigning them to the clusters defined by the GMM model using a threshold probability of p = 0.65. Otherwise, cells were considered to be in transition between clusters (unclassified transitioning, UT). The stochastic dynamic nature of the system implies that there will always be some cells outside the attractors that cannot be unequivocally classified by this method even when the system is at equilibrium and especially when it is close to bifurcation. The process was insensitive to the choice of threshold probability p (Methods S1 Section 2).

12Once cells had been allocated to clusters, we analyzed marker proteins and assigned cell identity accordingly (Figures 2B and 2D). The clustering was unsupervised and unambiguous. There was a good correspondence between the clusters found by the algorithm and the cell-type-specific marker expression known to be present at each day of the three reference signaling regimes (Figure 1A). The number of clusters needed to obtain unimodal distributions was larger for D3-D4.5, where there was more than one cluster corresponding to the same cell identity (Figure 2C). For example, at D3 and D3.5 there are two clusters corresponding to CE cells: one defined by BRA+ and CDX2+ cells (so-called neuromesodermal progenitors, ) and one with only CDX2+ cells (Figures 2B and 2C; Methods S1 Section 2).

13The clustering algorithm identified the expected cell states in each of the datasets (see Table S2). In line with our preliminary analysis, we were able to assign FOXC2+ LM identity on the basis of a SOX2-/TBX6- marker profile (Figures S1B and S2C). Moreover, the clustering algorithm distinguished between OTX2+ AN and OTX2- PN progenitors based on their distinct levels of SOX1 and SOX2, as we had observed previously (Figure S2B). However, as neither OTX2 nor SOX1/SOX2 expression levels were able to distinguish between anterior and posterior neural progenitors at D4.5 and D5, we extrapolated the proportions obtained at D4 to quantify their proportions in this case.

14At D2.5 under FGF + CHIR signaling conditions, an additional population of cells with intermediate levels of SOX2 and low expression of BRA, CDX2, SOX1, and TBX6 was detected by the clustering (Figure S2A). As these cells had begun to downregulate SOX2, this indicated that they were transitioning from EPI to CE but had yet to acquire CDX2+ posterior identity. We considered them as a separate group and labeled them transitioning (Tr). Thus, we were able to define clusters, corresponding to the attractors in the landscape, representing epiblast (EPI), anterior neural (AN), CE, posterior neural (PN), and mesoderm (M), and an additional transitioning population (Tr) using a minimal set of 5 marker proteins (BRA, CDX2, SOX1, SOX2, and TBX6) (Figure 2D).

15To test the validity of our cluster assignments and check for missed states, we examined the correlation of protein expression within a cluster and compared our clusters with the cell types found in. Principal component analysis of individual clusters identified distinct associations of proteins with principal components for each cluster (Figure 2E). Each principal component was associated with one or two proteins (for a full analysis see Methods S1 Section 2). To test the hypothesis of multivariate normality, we calculated local protein-protein correlations within each cluster and found smoothly varying patterns (Figure 2F) in the correlation for each cluster. Correlations between protein expression reflect regulatory interactions between the components of the underlying gene regulatory network. Hence, consistency of the correlation structure across each cluster and differences in the structure in different clusters is evidence that the clusters correspond to distinct biological states and that we are not missing any states since, in that case, discontinuities in the correlation structure within clusters would be apparent. Moreover, the boundaries between adjacent clusters were characterized by abrupt changes in this correlation structure (Figure 2F; Methods S1 Section 2). This is consistent with the idea that cell states represent high-dimensional attractors in protein expression space that are separated by saddle points. Previously, hierarchical clustering of transcriptome data during this differentiation produced five major groups representing distinct patterns of gene expression. Even though our choice of clusters was unsupervised, they mapped onto these groups indicating that incorporation of a broader set of genes does not identify any new states.

16The data obtained in our initial experimental series (Figure 3A) indicate that, in the absence of CHIR, EPI cells transition to AN after removal of FGF on D3, suggesting that FGF withdrawal causes a loss of the EPI attractor in such a way that cells escape toward the AN attractor (Figure 3Ca, no CHIR). On the other hand, if CHIR is added on D2 a substantial proportion of cells transition to CE by D3 (Figure 3Ca, CHIR 2–3). This suggests that, with the addition of CHIR, the disappearance of the EPI attractor allows cells to escape toward the CE attractor. The distinct fate outcomes associated with these decisions suggests that each is caused by the bifurcation of the EPI attractor but with escape routes in different directions. In both cases, all cells escape in the same direction toward the same attractor.

17Our approach to model selection used the mathematical classification of the simplest generic 3-attractor landscapes (; Methods S1 Section 3) and the assumption that the system should be part of one of these archetypal landscapes. These can each be modeled by relatively simple two-dimensional dynamical systems that will capture high-dimensional transition dynamics. Only one of these is compatible with the above observations. This is related to Thom’s butterfly catastrophe: the three attractors, AN, EPI, and CE, are separated by two saddle points (Figure 3D) with one of them in the middle, EPI. We call this the binary choice landscape. Changing levels of signal cause bifurcations: cells in the central attractor EPI can transition to either CE or AN (Figures 3Da and 3Dc) but those in states CE and AN are only able to transition to EPI (Figure 3Db). This models an one-or-other/all-or-nothing decision where a population of cells in EPI chooses between the fates CE or AN. If cells commence the transition from EPI to CE and signals change pushing them toward AN, they need to return to EPI first before transitioning to AN. This last feature is observed in our data when CHIR is removed after 24 h (Figure 3Ca).

18Following the transition to CE, the data suggest a second decision associated with cells leaving the CE attractor adopting a mixture of M and PN fates. The ratio of M:PN cells increases with the length of CHIR exposure and is not an all-or-nothing response as in the previous decision (Figure 3Cb). Moreover, the way in which cells leave the CE state appears markedly different from the way that they leave the EPI state (Figure 3Cb, CHIR 2–5). Cells appear to leave the CE state at an approximately constant rate in CHIR induction conditions. This raises the possibility that transitions are caused not by bifurcation but by fluctuation driven escape from the CE attractor basin.

19Among the archetypal landscapes, two contain a part that is compatible with these data and these two parts are equivalent (Figure 3E). We call this the binary flip landscape because it allows signals to flip the escape route of cells leaving the CE state, tipping them into either state M or PN (Figures 3Ea and 3Ec). When combined with stochastic noise, this allows cells leaving CE to be distributed with different ratios between M and PN (Figure 3Eb). Finally, we connect the two landscapes using the common CE attractor (Figure 3F). This is done using a transition function to allow the cell’s trajectory to transition smoothly from one landscape to the other (Methods S1 Section 3).

20To construct models, we took advantage of the normal forms of these two landscapes provided by catastrophe and dynamical systems theory as these minimize the number of parameters and state variables needed to capture the changing geometries. In addition, we postulated a linear relationship between the signals and model parameters and later validated this using the data (Figure 3G). All cells in a sample were assumed to be affected by the same landscape, determined by independent effects of each signal (CHIR and FGF). To fit the data, we found it necessary to add a memory effect to this relationship (Figure 3H). The data showed that removing CHIR at D3 or D3.5 gave different PN:M balances, while the PN:M balance remained near invariant if CHIR was removed at D4 or kept until D5 (Figure 3Cb, CHIR 2–4 and CHIR 2–5). This indicated that the effect of CHIR persisted for some time after its removal. To account for this memory effect, we maintained a CHIR signal in the model for a period of time proportional to the time CHIR had been in the system longer than a time threshold which was estimated in the fitting (Video S2). For FGF, since FGF signaling has been shown to generate a positive feedback and induce expression of FGF ligands, we maintained a level of 90% FGF when exogenous FGF was withdrawn, unless the inhibitor PD0325901 was added, in which case FGF signaling was set to 0.

21The resulting system depends on 20 parameters: 3 × 6 for the relationship from the 3 signal changes. No signaling (no CHIR + PD), addition of CHIR, and addition of FGF: the parameters w in Figure 3G) to the 6 landscape model parameters, a memory threshold and a noise level (Methods S1 Section 3). In order to fit these, we used a training dataset of 7 treatment regimens (see Figures 3A and 4A) measured at 7 time points that distinguish 7 populations: EPI, Tr, AN, CE, PN, M, and UT (a total of 343 measurements of cell fate proportions). We kept separate, for validation, a dataset of 4 experimental conditions with a total of 196 measurements (see Figures 3A and 4C). We took advantage of the ABC framework to generate posterior distributions of the parameters (see Methods S1 Section 4). Assessing correlation between the parameters indicated that there were few redundancies between them.

22By extensive sampling from the resulting posterior parameter distribution for each signal combination, we obtained simulation time series for cells in the landscape using the flow defined by the stochastic dynamical system (Figure 4D). We computed the proportion of cells around each attractor and compared these with the corresponding experimental proportions (see Methods S1 Section 5 for details). For each timepoint and signal combination, this produced a distribution of cell states and resulted in a 2×3 table of landscapes (Figure S3) corresponding to the different signal combinations used in the training experiments: CHIR on/off and FGF on/off or inhibited by PD.

23Even though the parameters come from a distribution, the distribution of cell states was relatively tight about their mean with an average coefficient of variation of 7%. The larger size of some coefficients of variation was due to misidentification of clusters for a small number of experimental conditions. To compare simulated and experimental proportions we, therefore, used the mean simulated proportion. Comparison with the training datasets showed excellent overall agreement (Figure 4A; see Methods S1 for a detailed analysis). The largest differences tended to be associated with situations with a significant number of unclassified cells in the experimental data (UT), which we reasoned were cells transitioning between states. The simulations reproduced subtle features of the data such as the constant rate at which cells leave the CE state in response to CHIR.

24We next examined simulations using signal combinations from the validation datasets computed with the model parameter obtained from the fitting to the training data. Again, there was good overall agreement (Figure 4C). The simulations showed that the cells transitioning at D2.5 were “recaptured” by the EPI attractor after removal of CHIR (Video S1), as observed in the experimental data, where a proportion of cells transitioning to CE at D2.5 return to an epiblast state at D3. Since the data used to fit the model did not contain any conditions in which the EPI attractor was repopulated after cells had left it, these simulations provide a non-trivial validation of the geometry of the landscape and its parameterization.

25Closer inspection, however, revealed some differences between simulations and experimental data. Most notably, there was an underestimation of the CE population in some CHIR conditions (CHIR 2–2.5, CHIR 2–3, CHIR 2–3.5), together with an underestimation of M and overestimation of PN. Taking the results together, there was good agreement between the simulations and experimental data with the model accurately replicating the key decision processes in the differentiation pathway, this included the recapturing of cells by the EPI attractor after short CHIR exposure and the leakage of CE cell into M. Nevertheless, discrepancies in the exact proportions of cells between simulations and experiment prompted us to investigate further refinement to the model fit.

26To test the model and determine the value of refining the fit, we took advantage of the model to design experiments. Simulations testing the effect of alternating periods of CHIR induction with periods of no CHIR suggested that the most informative results were obtained when cells were induced into the transition state by exposure to CHIR for 12 h at D2 and then subjected to either 1 or 2 pulses comprising 3 h with no CHIR followed by 3 h CHIR. After the pulses, cells were kept in continuous CHIR (Figure 5A). Simulations of the single pulse predicted a mixture of anterior and posterior fates at D5 (Figure 5B, Initial Sim). By contrast, simulations of the double pulse predicted approximately 50% of the cells would adopt an AN identity and almost no cells would remain CE at D5 (Figure 5B, Initial Sim). That is, the model predicted that 3 h of no CHIR and no FGF would be sufficient for cells to become AN and that two 3 h pulses of no CHIR conditions would be sufficient for CE cells to become PN. Moreover, the model predicted that the transitioning population (Tr) would be present at day 3 under both pulsing conditions.

27We performed these two experiments (Figure 5B, Test Exp) as part of the test experimental series. As expected, we observed that the transitioning population was present at D3 in both conditions. At later days the experimental data differed from our predictions. First, few cells adopted an AN identity, and second, a substantial CE population remained throughout D3.5 to D5. The continued presence of CE cells at D5 was consistent with the idea that CE represents an attractor in the dynamical landscape and it pointed to the initial fitting underestimating its stability. The discrepancy between simulations and experiments prompted us to refine the model fit.

28To this end, we included the data from the pulsing conditions together with the reference set of conditions to retrain the model (see Methods S1 Section 6 for details). Adjusting these parameters markedly improved the performance of the model (Figures 5B and 5D, Refined Sim). The AN population under the two pulse conditions was almost non-existent and the CE population was increased, matching the experimental data. Data from CHIR 2–2.5 condition that was not used in the training was used as a validation and showed a good agreement with the experimental data. Comparison of the refined landscape (Figure 5E) with the initially fitted landscape (Figure S3) revealed that the CE attractor did not bifurcate under no CHIR + FGF conditions and was very close to bifurcation under no CHIR + No FGF conditions allowing cells to remain in that area for longer. Moreover, the saddle point separating the AN attractor from the rest of the landscape had moved closer to the AN attractor, increasing the amount of time necessary for cells to commit to the AN fate under no CHIR + No FGF conditions.

29We conclude that in all conditions, the CE attractor is close to its bifurcation point. As a result, when it is a genuine attractor, its basin is shallow, hence it is easy for noise driven fluctuations to allow a cell’s trajectory to escape and proceed downhill to the PN or M attractors (Video S2). Thus, we see a gradual loss of CE cells at a roughly constant rate per unit time.

30The fitted parameterized landscape provides insight into the capacity for cells to revert to a previous state: in developmental biology terminology the commitment of cells. The model parameters determined by the fit mean that combinations of FGF and WNT signaling produce a constrained set of landscapes. This applies not just to the signal combinations used in the experiments but other reasonably conceivable combinations. This constraint means that certain transitions are not possible. For example, once cells have transitioned to the AN state, this is irreversible, as this attractor is deep and no bifurcation is allowed that destroys this state, thus cells are committed to AN and can no longer become CE in response to WNT signaling (Video S4). Similarly, the bifurcation of the EPI attractor in response to CHIR provides insight into the inability of CHIR + FGF inhibition with PD to revert epiblast cells to the ESC state (Video S3).

31Moreover, the model predicted a landscape for no CHIR with PD (Figure 5E, top right), a condition not used in the experiments (Figure 7). The model can also be used to predict the effect of different timings of signal changes (Figure 5) and extrapolated to different concentrations of signals (Figure 6).

32To investigate how the balance between M and PN was controlled in the second decision, we performed dose-response experiments with CHIR after D3. We first simulated the model assuming a linear relationship between the effective CHIR signal level and CHIR concentrations (Figure 6). Experimental results suggested a good correspondence for the low concentration (0.5 μM, 10% of saturated), but larger differences were evident for higher concentrations (2.5 μM, 50% of saturated) with a difference of 40% in PN, overestimating the CE and M proportions.

33These results suggested that the model needed to be further refined to account for the dose-response effects of CHIR concentration. We tested whether the addition of a sigmoidal dose-response curve improved the model performance. For the fitting of the dose-response effect, we kept the distribution of all other parameters the same and fit the parameters of the sigmoidal curve. This substantially improved model performance. Although the half-response parameter did not produce a clear distribution, suggesting that additional concentrations would be necessary for a precise estimation of the dose-response curve, the estimation of cell fates in response to both 1.7 μM (30% of saturated) and 2.5 μM were close to the empirical data.

34Having fit, refined, and validated the landscape, we set out to test its accuracy and investigate its implications by identifying predictions for which we had not previously analyzed experimental data.

35First, we noted that, in contrast to the EPI attractor that bifurcates with the saddle connected to the CE attractor in response to FGF + CHIR, the AN attractor was not destroyed by the addition of FGF + CHIR—a saddle remains separating AN and CE (Figure 5E). This predicted that once cells had committed to the AN identity, they would remain anterior and not adopt a CE identity after the addition of CHIR (Figure 7Ba). We tested this hypothesis by adding FGF + CHIR at D4 to cells that had adopted an AN identity. In the standard AN condition (no CHIR + FGF D0-3), cells remained anterior in response to addition of FGF + CHIR at D4 as indicated by the high levels of expression of the marker OTX2 at D5 (Figure 7Bb). These results were consistent with the constructed landscape, confirming that the saddle separating AN from CE does not bifurcate away in response to FGF + CHIR, and indicated that once EPI cells are committed, they remain anterior.

36Second, because the landscapes define a family parametrized by two signals (CHIR and FGF) they include a hypothesized landscape for no CHIR + PD conditions, even though none of the datasets used to construct the model corresponded to this condition (Figure 5E, top right). In this case, the model predicts that, in the absence of WNT, the inhibition of FGF signaling results in a further destabilization of the CE attractor and an expansion of the PN attractor (Figure 7F). We tested this by removing CHIR and adding the FGF inhibitor, PD, to cells that had adopted a CE identity in response to different durations of CHIR exposure. The prediction was that cells in the CE attractor at the time of CHIR removal would transition to PN, while cells that had already adopted an M identity would remain in this state. Consistent with these predictions, the removal of CHIR and inhibition of FGF signaling after 24h (CHIR 2−3 + PD 3–4) resulted in a moderate increase in the PN population compared with control conditions. Maintaining the CHIR for 36h (CHIR 2−3.5 + PD 3.5–4.5) resulted in a marked increase in PN population at the expense of CE identity compared with control conditions (Figure 7D), whereas the M population was less affected. This observation supports the idea that the CE represents a relatively shallow attractor that is maintained by CHIR signaling and destabilized by the removal or inhibition of FGF.

37In the embryo, retinoic acid (RA) emanating from the somites, is proposed to inhibit FGF signaling in differentiating CE cells to promote neural differentiation. We, therefore, compared the effect of RA addition with the effect of FGF inhibition. The addition of RA after 24 h of CHIR (CHIR 2−3 + RA 3–4) or 36 h (CHIR 2−3.5 + RA 3.5–4.5) resulted in almost identical proportions of cell types as obtained for the equivalent experiments with PD (Figure S4). Moreover, we also predicted that if cells were exposed to CHIR from D2-4, the memory effect of CHIR would result in the maintenance of a CE population at D5 even after the addition of PD or RA. Indeed, exposure of cells to RA after 48 h of CHIR induction (CHIR 2−4 + RA 4–5) agreed with model predictions, with a large number of cells adopting M identity by D5. In these experiments a substantial proportion of cells remains CE at 24 h after removal of CHIR and addition of RA, consistent with the memory effect of CHIR hypothesized for the construction of the model.

38Taken together, the results provide a quantitative model for the differentiation of neural and mesodermal cells from pluripotent progenitors that accurately predicts the effect of signals on the proportions of cell types generated. Moreover, the quantified landscape offers an intuitive visualization of the differentiation process that highlights the distinctions between two types of decisions (binary choice and binary flip) that would be obscured in more complex gene-centric models.

39It has become increasingly easy to control cell fate and to generate quantitative data, raising the prospect of predictive engineering of organs and tissues. To support this endeavor, quantitative models that assimilate data and represent the successive steps in development are needed. Models are commonly constructed around genes and their interactions, typically defined with mass action kinetics, but their complexity and the large number of parameters involved can limit the scope of such models and obscure underlying decision mechanisms. The Waddington landscape metaphor is an appealing simplification, yet until recently it lacked both a mathematical foundation and meaningful engagement with data. We demonstrate an approach to construct, parameterize, and analyze dynamical landscapes applicable to any system with discrete quantifiable states (see Methods S1 Section 8 for a description of the general procedure). We show that it makes specific and testable quantitative predictions as well as providing a visualization of the developmental process that offers insight into the underlying mechanisms. We emphasize that the modeling strategy we have developed allows the targeted exploration of specific features of a cellular differentiation process allowing the design of an optimized combination of experimental conditions. This work revealed two distinct strategies by which cells decide between one of two fates. We propose that they represent archetypes for developmental decisions that play a general role in developmental dynamics.

40Gene network models and real cellular systems develop by successive transitions among a small panel of fates, so the mathematics of dynamical systems and catastrophe theory can be applied. Moreover, there is formal mathematical theory for such systems that reflects, extends, and clarifies Waddington’s picture. This theory indicates that even when the state space (gene or protein expression) is high dimensional, developmental transitions comprise the simplest saddles (those with one-dimensional unstable manifolds) and their attractors. Consequently, the topology and behavior of such models can be represented in two dimensions, so that it is conceptually meaningful to work in the low-dimensional space and assert that relevant biology takes place on a two-dimensional manifold. The final crucial step is to make the model quantitative by fitting a rich set of its summary statistics to those of the data. This links the powerful qualitative theories to quantitative models and provides an intuitive but rigorous means to analyze developmental dynamics.

41We applied these ideas to the differentiation of ESCs in response to two signals, WNT and FGF, that together direct differentiation into precursors of neural and mesodermal tissues via intermediates that have received prior study. This system has the virtue that cell identity is controlled by external signals, hence can be modeled without the complexity of cell interactions. Landscape models describe the cells as they transition between valleys in the Waddington metaphor and the most informative data are obtained during transitions. We show that quantitative data using just five molecular markers and a carefully designed array of conditions and times was sufficient to define the cellular states and obtain informative summary statistics: the proportions of cells of a given type at each time point. This information was sufficient to identify the underlying geometrical model. The fixed points of the model were abstract and acquired meaning when gene expression patterns were related to in vivo experiments. The morphogens control the basins of attraction of the fixed points and ultimately the saddle node bifurcations or flips that define cellular decisions. All the relevant cellular decisions could be captured in two dimensions. This strategy for constructing models of cellular differentiation is also compatible with single-cell transcriptome analyses of tissue development that simultaneously profile thousands of genes as long as the data are sufficient to identify the discrete cell states and cell fate transitions of the developmental process.

42The results support the notion that cell states correspond to attractors of a dynamical system. For the EPI, AN, PN, and M states, their assignment as attractors was unambiguous and the differentiation of cells between states involved bifurcations of attractors. For CE, the results were more subtle and revealing. For some morphogen combinations, the CE state appeared as an attractor, albeit a shallow one with a relatively small basin of attraction close to bifurcation, while for other morphogen combinations, the CE state had bifurcated away. The consequence of the shallowness of the attractor is a gradual leakage of cells from the CE state and this, in contrast to the bifurcation of an attractor, explains the transitory nature of the CE state. Cells leave the CE attractor by jumping over a saddle into the basin of a different attractor. This observation not only illustrates the subtlety of the method but also provides an explanation for the progressive elaboration of the anterior-posterior axis of vertebrate embryos. This process involves the production of trunk mesoderm and spinal cord over an extended period of developmental time from bipotential progenitors (CE) located at the tail end of the elongating body. Ensuring the CE state is poised close to bifurcation allows the constant differentiation of a fraction of cells to neural and mesodermal tissue while retaining an uncommitted progenitor population.

43The model we constructed terminates in one of three fates: mesoderm, AN, or PN. Cells in these states are committed, they cannot leave these attractors and remerge on our fate plane, but they are not endpoints as cells continue to differentiate. We envision that longer differentiation trajectories could be constructed, branch by branch, with the modular and hierarchical property of the models facilitating the gradual elaboration of the full differentiation tree. By fitting and refining the saddles and attractors as they move in response to extrinsic signals, such models generate quantitative predictions and suggest complex temporal stimuli that might stabilize or disrupt populations of cells in the dynamical landscape. In this way, landscape models delineate phenomenological features of experimental embryology in quantitative form.

44The landscape that resulted from our analysis comprises two sub-landscapes. These are generic in the sense that they occur naturally in dynamical systems without any special conditions and when they occur, they are structurally stable and robust to misspecification or perturbation of the dynamical system. Each encodes a fundamentally different type of binary decision. For the AN-EPI-CE decision, the data support the notion that the EPI state is destabilized by collision of a saddle point between it and either the CE or AN attractor. After either bifurcation, the saddle separating the EPI attractor from the alternative fate remains intact. The result is an all-or-nothing choice that commits cells to a fate. This agrees with previous experimental observations that epiblast cells commit to an anterior or posterior identity prior to acquiring a neural fate. Since the overall effect of this landscape is that signaling results in all cells in the EPI state choosing and committing to either the CE or AN state, we call this landscape the binary choice landscape.

45The decision for the differentiation of CE to either PN or M is different. The landscape indicates that the proportion of PN and M cells generated in response to different amounts of WNT signaling is a consequence of the location of the unstable manifold of the saddle that separates the CE attractor from PN and M. For prolonged and high CHIR conditions, the escape route connects to the M attractor, whereas in the absence of WNT signaling it flips to PN. For intermediate conditions stochasticity produces a mixed outcome with a progressive increase in the ratio of M to PN as WNT signaling increases. Because the flip in the unstable manifold is the essential dynamical phenomena controlling the allocation of cells, we call this the binary flip landscape. Previously, a pitchfork bifurcation had been postulated for this transition. However, in contrast to our landscapes, this is not generic but only occurs with an additional constraint equivalent to symmetry between M and PN and it is reversible such that cell states are only maintained if signaling is maintained.

46The results help clarify the idea of cell commitment, by explaining when cells that have escaped a basin of a bifurcated attractor can or cannot be recaptured by reversing the signal and reintroducing the attractor. For the EPI to CE transition, recapture is possible for a limited time after which they are committed to the transition. For the EPI to AN transition, once cells have entered the AN state they cannot escape because the AN attractor never bifurcates away. In this view, therefore, commitment becomes a dynamical property of the landscape.

47Theoretical considerations indicate that there are relatively few generic three-attractor landscapes. Thus, the binary choice and binary flip landscapes are likely to be encountered repeatedly in developmental decisions, in particular, in the common situation where a progenitor chooses between two alternative fates. Moreover, as we show here, these landscapes can be linked together to provide an ordered hierarchy that accounts for a differentiation pathway comprising multiple decisions. Together this leads us to suggest that the binary choice and binary flip landscapes represent design archetypes that underpin cellular decisions throughout embryonic development.

48In our approach, we did not fit flow-cytometry data directly to the dynamical variables but instead relied on the assignment of cells to discrete fates. This assignment is non-trivial and may be ambiguous in some systems. We sampled the population frequently in time to capture the transitions among the fates in the topological model and, thus, obtained quantitative predictions about what matters most, the discrete cell states. A future goal will be to develop tools that directly link gene expression measurements and detailed gene centered models to landscape models especially in order to detect the transitions in detail. The minimal spatial organization in the differentiating stem cell colonies, meant only the response of individual cells to exogenous signals needed to be considered. To apply these models to cells of an intact embryo, we would surmise that the cell intrinsic landscape remains the same, but the spatial-temporal behavior of the morphogens would need to be measured or inferred. To probe such models, it will be necessary to perturb signaling pathways with sufficient temporal resolution to catch cell fate transitions, experimental manipulations that do not control time are manifestly less useful. Thus, one can hope to reconstruct development in space-time analogous to the reconstruction of differentiation from single-cell expression data (e.g.,; ) but in a way that provides quantitative predictions and insight into the underlying mechanisms.

49Further information and requests for resources and reagents should be directed to and will be fulfilled by the lead contact, James Briscoe (James.Briscoe@crick.ac.uk).

50This study did not generate new unique reagents

51The XY mouse ES HM1 TetON line was used for all experiments. ES cells were maintained at 37°C with 5% carbon dioxide (CO2) and routinely tested for mycoplasma.

52All mouse ESCs were propagated on mitotically inactivated mouse embryonic fibroblasts (feeders) in DMEM knockout medium supplemented with 1000U/ml LIF (Chemicon), 10% cell-culture validated fetal bovine serum, penicillin/streptomycin, 2mM L-glutamine (GIBCO). To obtain EPI and CE, ESCs were differentiated as previously described with the addition of the porcupine inhibitor LGK974 (Cayman) in all culture medium. Briefly, ESCs were dissociated with 0.05% trypsin, and plated on tissue-culture treated plates for two sequential 20-minute periods in ESC medium to separate them from their feeder layer cells which adhere to the plastic. To start the differentiation, cells remaining in the supernatant were pelleted by centrifugation, counted, and resuspended in N2B27 medium containing 10ng/ml bFGF (Peprotech) + 5μM LGK974, and 50,000 cells per 35mm gelatin-coated CELLBIND dish (Corning) were plated. N2B27 medium contained at 1:1 ratio of DMEM/F12:Neurobasal medium (GIBCO) supplemented with 1xN2 (GIBCO), 1xB27 (GIBCO), 2mM L-glutamine (GIBCO), 40mg/ml BSA (Sigma), penicillin/streptomycin and 0.1mM 2-mercaptoethanol. To generate epiblast (Epi) cells, the cells were grown for 72 hrs in N2B27 + 10ng/ml bFGF + 5μM LGK974. To generate caudal epiblast (CE), cells were cultured with N2B27 + 10ng/ml bFGF + 5μM LGK974 for 48 hrs, then N2B27 + 10ng/ml bFGF + 5μM LGK974 + 5μM CHIR99021 (Axon) (FLC-medium) for a further 24hrs. CE were differentiated to posterior neural progenitors (NP) by removal of bFGF and CHIR from culture medium at 72hrs, and to paraxial mesoderm (M) by removal of bFGF and maintenance of 5μM CHIR from 72hrs onwards. When investigating the role of endogenous FGF production either the MEK inhibitor PD0325901 (500nM) was employed to inhibit downstream FGF signalling or retinoic acid (10nM) was added to inhibit expression of FGF. For all experiments described, cells were cultured for 48hrs before changing medium. Medium changes were then made every 12 hours.

53Cells were washed in PBS and fixed in 4% paraformaldehyde in PBS for 15min at 4°C, followed by two washes in PBS and one wash in PBST (0.1% Triton X-100 diluted in PBS). Primary antibodies were applied overnight at 4°C diluted in filter-sterilized blocking solution (1% BSA diluted in PBST). Cells were washed 3x in PBST and incubated with secondary antibodies at room temperature, for 1hr. Cells were washed 3x in PBST, incubated with DAPI for 5 min in PBS and washed twice before mounting with Prolong Gold (Invitrogen). Cells were imaged on a Zeiss Imager.Z2 microscope using the ApoTome.2 structured illumination platform. Z stacks were acquired and represented as maximum intensity projections using Fiji software. The same settings were applied to all images. Immunofluorescence was performed on a minimum of 2 biological replicates, from independent experiments. Primary antibodies used: SOX2 (Mouse, Santa Cruz sc-365823), OTX2 (Goat, R&D AF1979), CDX2 (Rabbit, Abcam ab76541), TBX6 (Goat, R&D AF4744), FOXC2 (Sheep, R&D AF6989). Secondary antibodies raised in donkey coupled to AlexaFluor- 488, 568, or 647 fluorophores (Molecular Probes) were used at 1:1000 dilution throughout. SOX1 was labelled with AlexaFluor-647 conjugated antibody (Mouse, BD 562224).

54Cells were washed in PBS and dissociated with minimal accutase (GIBCO). Once detached cells were collected into 1.5mL Eppendorf tubes by dissociating in N2B27 and pelleted. Cells were resuspended in PBS, pelleted and resuspended in 4% paraformaldehyde in PBS. Following 15min incubation at 4 C, cells were centrifuged, resuspended in PBS, and stored at 4°C for future analysis. On the day of flow cytometry, cells were counted and equal cell numbers were transferred for staining in v-bottom 96 well plates. Samples were pelleted and resuspended in 50μL FACS block (PBS + 0.2% Triton + 3% BSA). After 10min incubation at room temperature antibodies were added to the sample and incubated overnight at 4°C. Cells were pelleted at 700rcf for 5min and resuspended in 50μL FACS block. One additional wash was performed before acquisition on a Fortessa flow cytometer (BD) equipped with a high-throughput sampler using FACSDiva software. Control samples were labelled with individual fluorophore-conjugated antibodies to compensate for fluorescence spillover. Commercially available antibodies used were SOX2-V450 (Mouse, BD-561610), SOX1-PerCP-Cy.5.5 (Mouse, BD 561549), CDX2-PE (Mouse, BD 563428), BRA-APC (Goat, R&D IC2085A). Anti TBX6 (Goat, R&D AF4744) was conjugated to AlexaFluor-488, and anti FOXC2 (Goat, R&D AF6989) was conjugated to AlexaFluor-647 using Molecular Probes Antibody Labelling Kit (Invitrogen/ ThermoFisher). When labelling OTX2, samples were incubated with unconjugated OTX2 antibody (Goat, R&D AF1979) for 1 hour at room temperature in FACS block. Following 2 washes with FACS block AlexaFluor-488 anti-goat antibody (1:1000) was included in overnight incubations with fluorescently conjugated mouse antibodies.

55Compensation was performed prior to data acquisition in FACSDiva. FCS files were imported into the R programming environment using the ‘flowCore’ package. Individual cells were gated using a standard two-step gating strategy with an additional gate to remove outlier events with especially low fluorescence values as shown in (Figure S5). A custom R script was used to remove outlier events with fluorescent values exceeding a 3-sigma threshold from the mean of the average calculated across all samples of a given experiment. Raw fluorescence-intensity values were used in all downstream analysis and graphical representations of the data.

56See Methods S1

57See Methods S1

58See Methods S1

59See Methods S1

60See Methods S1

61See Methods S1

62See Methods S1

63See Methods S1