Box 1.– A dynamical-systems view of cellular decision-making: attractors and saddles. Terms explained: attractors, attractor clusters, saddles, unstable manifold, transition route, near-bifurcation escape, fold bifurcation, basins of attraction, microheterogeneity, landscape.

We interpret cell fate decisions using dynamical systems theory applied to gene regulatory networks (GRNs). Molecular interactions define a stochastic dynamical system governing cell state evolution, while external signals reshape the available states and transitions by modulating system parameters. The central structures are attractors, saddles, and the transition routes connecting them, together forming a dynamical landscape.

The GRNs considered here are strongly dissipative, so trajectories converge to stable fixed points corresponding to cell fates; we do not consider oscillatory states. Although gene expression is noisy, the noise is usually insufficient to destroy the underlying deterministic structure. We therefore treat cellular dynamics as a stochastic perturbation of a deterministic system: the deterministic component defines the possible states and connections, while stochasticity influences the timing and route of transitions.

Attractors and cell fates. Stable cell fates correspond to deterministic attractors (Fig 1A). Noise causes cells to fluctuate around them, producing clouds of gene-expression profiles we term attractor clusters (ACs) (Fig 1B).

We identify ACs through properties expected of attractors: incoming transition routes, persistence across time points, approximately unimodal expression distributions, loss of stability via characteristic bifurcations, and transitions to neighbouring ACs.

Saddles and transition routes. Certain saddles act as “mountain passes” enabling transitions between fates. We focus on saddles with a one-dimensional unstable manifold, which act as passes in the landscape. Transitions occur either through near-bifurcation escape, where stochastic fluctuations drive escape, or when a signal change causes the saddle and attractor to collide in a fold bifurcation, destroying the attractor. In both cases the most likely route over the saddle follows the saddle’s unstable manifold—the curve of points whose infinite-backward-time limit is the saddle. Such a well-defined route also exists when the attractor is destroyed in a fold.

Because the system is stochastic, cells only approximately follow this unstable manifold and in data this single path is replaced by a transition route that approximates it.

The proximity of an attractor and its neighbouring saddle determines the attractor’s depth. Cells in deep attractors require large perturbations to escape and therefore represent robust fates, whereas cells in shallow attractors are easily destabilised. Signals alter depth by changing GRN parameters, allowing stable fates to become unstable.

Stable manifolds and basins of attraction. The complementary stable manifolds (green curves in Fig 1A) partition the phase space into basins of attraction, regions whose trajectories converge to the same attractor (coloured regions in Fig 1A). Crossing a stable manifold corresponds to a commitment event.

Because cells in a given cell state have a slightly different molecular state, cell populations contain microheterogeneous landscapes. Apparently similar cells may therefore adopt different fates not only through stochastic fluctuations but also because they occupy different sides of basin boundaries or bifurcation curves (Box 3). The extra stochasticity introduced by microheterogeneity further smears the ACs and transition routes.

Head attractor clusters. At branching decisions there is always an upstream AC from which cells can access multiple downstream fates; we term this the head AC. Microheterogeneity and signalling differences within the head AC generate distinct escape trajectories, producing branching behaviour (Fig 1B).

Landscapes and topology. The attractors, saddles, and connecting unstable manifolds together define the dynamical landscape (Fig 1A). Its graph structure—with attractors as nodes and transition routes as edges—can possess nontrivial topology.

Fig 2. Attractor clusters in single-cell RNA sequencing data reveal discrete cell states during neural development.

A. Schematic diagram of neural tube patterning regulated by the Shh gradient. B. Schematic of in vitro differentiation of mESCs and timing of the scRNA-seq data collection (D, Day). C. Example of the outward clustering methodology illustrated by analysis of the late time points (D5-D7). (Left) 2D PCA projection of cells that express one of the marker genes Olig2, Nkx2.2 and Shh. (Middle) PCA projection from d-dimensional effective gene space using the d most significantly differentially expressed genes between the 3 samples of cells (here d = 112). Cells are clustered using Leiden clustering (a graph-based algorithm for identifying densely connected groups of cells) with fine resolution and coloured accordingly. (Right) The fine Leiden clusters are merged based on their molecular identities. D. The most significantly differentially expressed genes distinguishing the 3 samples of cells shown in C (left panel). E. Extraction of the group of adjacent ACs, Early p3, p3 and FP, visualised in 2D LDA space. Histograms for each linear discriminant score (i.e., the dot product of the data vector with the ith LDA eigenvector) show that ACs (such as FP) exhibit a well-defined approximately Gaussian structure in LDA space. F. Temporal progression of Early p3, p3, and FP identifies Early p3 as the ("head") state from which cells transition to either p3 or FP. Grey dots represent all cells in the group (D5-D7), coloured dots indicate cells present at the indicated time points with the colour indicating their cluster assignment. Contour lines indicate the density of cells at each time point. G–H. Identification of early ACs (D3-D4) and transitioning cells between them labelled -Trans. G. PCA projection of cells from D3-D4 samples showing the early ACs. H. Dot plot indicating the expression of markers used to identify the early ACs. The size of the circle represents the percentage of cells in the AC and the intensity of shading indicates the mean expression of the gene. I. The proportion of cells (in percent) in each AC and transitioning clusters at the early time points (D3-D4). J–L. Identification of late ACs (D4-D7) and transitioning cells between them, labelled as -Trans. Similar to G–I but using cells from the D4-D7 samples. The D8 sample is not displayed; it predominantly contains MNs, V3 neurones, MNDiff cells, and p3 cells. The data underlying this Figure can be found at https://doi.org/10.5281/zenodo.15584010.

https://doi.org/10.1371/journal.pbio.3003953.g002

Fig 3. Neural tube decision landscape.

Schematic diagram summarising all the ACs associated with the cell states found in the scRNA-seq dataset at 500 nM SAG across all time points (D3–D7). This is included to introduce the complete set of states in the landscape and the possible transitions between them. These will be discussed in later sections.

https://doi.org/10.1371/journal.pbio.3003953.g003

Outward clustering identifies cell states through iterative gene selection.

We developed an iterative gene selection and clustering approach that we term outward clustering: beginning from a small set of known marker genes, the method iteratively expands—or “clusters outward”—by identifying genes that are differentially expressed between the clusters defined at each step, progressively building up the effective gene space without presupposing any cell states in advance. The approach is based on the principle that genes relevant to the GRN should show differential expression between cell states or transition pathways (detailed in Section A1.1 in S1 Appendix). Since cell states are regarded as initially unknown, we bootstrapped the analysis using a small number of marker genes (typically three or four) that define known cell types, regions or pathways. For early time points (D3-D4), we used TBXT/BRA, Foxc2, Foxa2, and Olig2, which identify neuromesodermal progenitors, mesodermal tissues, and early neural lineages [36]. We identified subsets of cells expressing just one of these markers and performed differential gene expression analysis between these subsets. This process generated gene modules of approximately 100 genes, determined by setting statistical significance thresholds. For later time points, we used Olig2, Nkx2.2, and Shh which are markers of pMN, p3, and FP progenitors and performed an equivalent analysis (Fig 2C, 2D).

Although the early and late gene modules produced in this way define relatively high-dimensional data spaces (which we call effective gene spaces), they are reduced in dimension by a factor of about 30 compared to the dimensionality of all highly variable genes. This allowed almost all analyses to be carried out in gene space, as advocated in [45], without resorting to low-dimensional projections, except for visualisation. This, for example, allowed us to compute cell fate differentiation paths in the effective gene spaces as opposed to tracing them on a 2D or 3D visualisation. The only exception was for the fine-resolution Leiden clustering for which we projected the data into a lower-dimensional PCA space [46]. Following this, we merged fine-scale clusters into cell state clusters based on the molecular identities that were identified (Fig 2C).

We confirmed that we had captured the genes that are differentially expressed between each cluster and the others. We checked for missed clusters by examining unassigned cells for potentially missed states and testing whether additional genes showed significant differential expression between cluster pairs. We also confirmed robustness by substituting alternative correlated marker genes and verifying that cluster structures remained consistent. After analysing transition pathways (described below), we expanded the gene set to include genes showing significant variation in transitions, along with expert-curated transcription factors. None of these additions substantially affected subsequent analyses.

Validation of attractor clusters.

In a deterministic system, a stable gene expression state (i.e., an attractor) would be a single fixed point: a precise combination of gene activity levels maintained by regulatory feedback. The stochastic nature of molecular dynamics and the data modality causes cells states to fluctuate around this point, replacing it with a cloud of similar expression profiles. We call this cloud an attractor cluster (AC). It has key properties that can be used to validate its status: it captures nearby states, has a well-defined unimodal statistical structure, will either be stable or bifurcate in characteristic ways, and will connect via saddle points to one or more upstream and downstream ACs.

For this validation we used multiple criteria. Firstly, we tracked each AC over time to assess its stability at successive time points, identify any loss of stability, and characterise resulting transitions to adjacent ACs. Secondly, in order to clarify the topological relationship between ACs and more clearly see the relationship between an AC and neighbouring ACs connected by transition routes, we extracted small groups of neighbouring ACs and their associated transitioning clusters – which we term sub-landscapes– and projected them into 2- or 3-dimensions using linear discriminant analysis (LDA) [47]. LDA maximises separation between clusters and faithfully represents them because, although LDA rescales axes, it is linear and therefore does not introduce artificial bends, folds, or discontinuities and faithfully maintains the topological relationships beetween states (and how routes connect them). To illustrate this approach, consider the previously unappreciated cell state represented by the Early p3 AC together with the p3 and FP ACs (Fig 2E). These together with the two transition routes from Early p3 make up such a sub-landscape where the Early p3 AC loses stability causing cells to transition to either p3 or FP, consistent with published observations [48] (Fig 2F).

Thirdly, Wilcoxon rank-sum tests confirmed that known progenitor marker genes distinguished ACs. Fourthly, we measured Kullback-Leibler divergence (a statistical measure of how well-separated two clusters are) in LDA space (Table A1 in S1 Appendix), which provides a lower bound for the same divergence measured in the effective gene space. Fifthly, cell distributions within each AC were checked to be unimodal along linear discriminant components, consistent with approximately linear dynamics near attractors (Fig 2E). Finally, we examined gene-gene correlations within and between ACs, with differing correlations providing evidence of distinct regulatory interactions in each cluster (Section A1.3 in S1 Appendix).

An important aspect of our analysis is the reproducibility of our determination of ACs. We have data across multiple timepoints and, for the flow cytometry data, which is discussed below, across multiple experimental conditions. Directly comparing the different ACs across these conditions indicated that they are reproducible, as can be appreciated by comparing data at a given time point and/or experimental condition with the position of the dataset consisting of all times and conditions. To visualise this when projecting into a lower dimensional space we always show the data of current interest against a grey plot of all the data in that sub-landscape. This reproducibility provides additional evidence to determine the persistence of ACs and the location of transition routes. These reproducibility checks are performed in gene expression space, LDA projections and on the standard 2D plots of expression levels of pairs of key genes. As an example consider Fig 5A, 5B below, showing the branching decision from NMP. This shows the overall structure, the emptying of the NMP state and filling of the Early Mesoderm and PreNeural ACs and the stability of their position over four time points.

This method identified all expected cell states plus several previously uncharacterised intermediates which span early mesoderm specification to ventral neural progenitors [34](Figs 2G – 2L and 3).

At early time points, ACs included NMPs [49] (TBXT/Bra), mesodermal cells [50] (Foxc2), PreNeural cells [43,51] (Nkx1.2, Msx1, Sox2), transitioning NMPs (NMPTrans), and Early Mesoderm [52,53] (Tbx6). As differentiation progressed, ventral neural progenitor states emerged: pMN [54] (Olig2), p3 [55] (Nkx2.2), and FP [56] (Foxa2, Arx, Shh). The p0/p1 (Pax6, Irx3) and p2 (Pax6, Irx3, Nkx6.1) ACs appeared at low numbers at D3.8 and D4, consistent with high SAG concentrations.

Additional ACs/cells states and genes that characterise them included MNDiff (Olig2, Neurog2), representing pMN progenitors undergoing motor neurone differentiation [44,57], and DP (Olig2, Nkx2.9, Hes1), representing an intermediate between p3 and pMN identities [58]. The EarlyVentral population (Nkx6.1, Foxa2) at D3.8-D4 represents a common precursor of FP and p3 [36]. At D5, we identified Early p3 (Nkx2.9), distinct from conventional p3 (Nkx2.2, Ascl1) observed at D6-D7 [48,59,60], and V3 neurones (Sim1, Nkx2.2), the differentiated p3 progeny [55].

Overall, this analysis demonstrates that outward clustering robustly identifies all known cell states and reveals developmental intermediates.

Direct transitions and binary decisions.

When a near-bifurcation escape takes place, the simplest outcome is that all cells leaving an attractor follow the same path to the same downstream attractor. These direct transitions involve no branching and commit cells to a single fate. The next case is a choice between two potential transition routes, whose minimal architecture is the binary flip landscape defined in Box 2. In such a landscape, when cells escape the head attractor the choice between fates depends on the signals received and the cell’s state, and changes in either can redirect the route to the alternative downstream attractor. Populations typically exit the head attractor via both routes, splitting into two groups with distinct gene-expression trajectories (Fig 5A). We identified the NMP, PreNeural, and Early p3 ACs as head ACs of binary flip sub-landscapes, while the p2 attractor heads a related but slightly more complex architecture (discussed below). The binary decision from NMP to either neural or mesodermal fates [34,35] is a canonical example: LDA projection shows clear separation between the NMP, PreNeural, Early Mesoderm, and Mesoderm attractors, with transitioning cells forming bridges between populations (Fig 5A). Between D3 and D3.6, cells diverge from NMP along two trajectories—one leading to PreNeural and the other to Early Mesoderm (D3.6), which then transitions directly to Mesoderm at D4 (Fig 5B). The PreNeural decision behaves analogously: cells with high Foxa2 expression take the ventral route while those with high Pax6 and low Foxa2 expression take the intermediate route [36] (Fig A5 in S1 Appendix).

Transition analysis.

To validate the bifurcation structures and identify the genes involved in the transitions, we developed a gene expression analysis using tubular neighbourhoods along splines that approximate the escape-trajectory manifolds described above (Section A2 and Fig A4C, A4D in S1 Appendix). Rather than relying on a single path connecting the source AC to a downstream AC, we repeatably subsample the transitioning cells, fit a spline approximating the mean pathway for the subsample, and compute the mean expression of the significantly varying genes within the corresponding tubular neighbourhoods along the spline. We then confirm that the resulting ensemble of time series produced in this way is reproducible across subsamples (Fig 5E) and in the case where branching routes are expected we demand these ensembles are internally consistent and clearly distinct from each other. We do this because cells are microheterogeneous (near-identical cells in or transitioning from a single AC are unlikely to have identical landscape) and therefore we require that all reasonable candidates for the unstable manifold have similar gene expression time series that are also clearly distinct from other routes. This approach detects genes that vary significantly during transitions - including genes that are not differentially expressed between ACs but have nontrivial transition dynamics, thereby identifying potential GRN components. In addition to the genes defining the effective gene space, we tracked approximately 150 curated transcription factors relevant to neural development.

To analyse branching cell state transitions from a given AC, we restrict attention to the cells within that AC and the neighbouring ACs to which they connect. As an illustration, the cells in Fig 5A suggest a potential branching from the NMP AC towards the PreNeural and Mesoderm ACs.

We apply two criteria to validate coherent temporal gene-expression patterns and to distinguish the proposed bifurcation from the alternatives. For the first criterion, we applied the subsampling and spline procedure described above, then verified that the resulting expression time series fell into two well-defined classes to which cells could be reliably assigned. For the NMP example, almost all relevant cells followed one of two well-defined pathways with clearly distinct temporal expression patterns (Fig 5C–5E). We also assessed whether the variation along the splines was consistent with the trajectories of a stochastic dynamical system: that is, whether the fluctuations were structured rather than dominated by excessive noise or high-frequency fluctuations. More details are given in Section A2.1, A2.3 in S1 Appendix.

For the second criterion we assess geometry of departure from the head AC (Box 2) to determine the architecture of the landscpae. For example, while in a binary flip, cells leave along a shared initial direction before diverging, while in a binary choice, they depart in different directions from the outset. In practice, we assessed this by examining whether the gene-expression profiles of transitioning cells at the earliest exit stage form a single cluster (flip), or two separate clusters (choice) (Box 2).

To these criteria we later add further validation of the proposed bifurcation structure by fitting a quantitative model to the data.

This approach was used to validate binary flip architectures, including escape from a head attractor immediately before branching. For the NMP head attractor, it revealed upregulation of Nkx1.2 and Irx3 for escape along the PreNeural pathway and of Tbx6 and Foxc2 for the Mesoderm pathway (Fig 5E), providing support for a binary flip in this case.

The nature of the transition is further revealed by the temporal progression of cells along the unstable manifolds of the flip landscape (Fig 5F), which shows how the distributions of cells within the tubular neighbourhood advances along the manifold over time. Here we observe a pulse of cells leaving an AC seemingly in unison. The plot would look different if cells escape at a roughly constant rate via stochastic fluctuations: the distribution would then remain stationary rather than progressing, providing a way to distinguish the two cases. It also gives timing information: here the broader variance at later stages (e.g., at D3.4) supports the idea that cells experience delays and do not transition simultaneously.

The critical transition from ACs was further validated using the transition index of Mojtahedi and colleagues [16,61,62], which tends to rise sharply at exit (e.g., Fig 5G for exit from NMP). We applied the same analyses to the other branching sub-landscapes (Section A2.4 in S1 Appendix).

Fig 6. Flow cytometry analysis of cell states reveals SAG concentration-dependent neural progenitor specification.

A. Schematic of the signalling conditions used for flow cytometry. B. 2D PCA projection of the ACs found using flow cytometry data. Data for all SAG concentrations and time points are combined. C. Dot plot indicating the expression of markers used to characterise the ACs in B (for the other clusters see Fig A8 B in S1 Appendix). The Floor Plate cluster (FP) was not identifiable in this dataset (see discussion in Section A4.2 in S1 Appendix). D. Percentage occupancy of the ACs as a function of Day (D) and SAG concentration (0, 10, 100, and 500 nM SAG). This shows the average proportions over 3 experimental replicates for 10 nM and 100 nM and 5 replicates for 0 and 500 nM (Section A3 in S1 Appendix). Blank spaces in bar plots correspond to the proportions of transitioning and unclassified cells, these are shown in Fig A8 C in S1 Appendix for all replicates. E. 3D representation in LDA coordinates of the clusters at all time points (unclassified cells excluded). At 0 nM SAG cells remain confined to the p0/p1 AC but start transitioning to p2 and pMN at 10 nM SAG. The ventral cell branch (green/cyan dots) that includes EarlyVentral, Early p3, and p3 only appears at high SAG concentrations (100, 500 nM). F. Schematic diagrams of attractors present at different SAG levels and the connections between them. This represents the consensus state because a landscape is associated with an individual cell, which are heterogeneous in state and signal. G, H. Temporal evolution of selected ACs (D4-D6) in LDA coordinates as a function of time and SAG concentration. Grey dots represent all cells in the group, coloured dots indicate cells present at the indicated time points with the colours indicating their cluster assignment. Contour lines indicate the density of cells at each time point. G. Binary flip decision from PreNeural to either p0/p1 or EarlyVentral. The figure shows 0, 500 nM SAG (see Fig A13 in S1 Appendix for other concentrations). (Top) The p0/p1 AC is stable when SAG is absent (0 nM) and there is no transition to EarlyVentral. (Bottom) The transition to EarlyVentral only appears at high SAG concentrations (100, 500 nM). PreNeural, p0/p1, and EarlyVentral empty over three days at the highest SAG level (500 nM). The attractors have bifurcated, but the unstable manifolds from PreNeural to pMN and Early p3 still pass through the regions in gene space formerly occupied by the bifurcated attractors. H. Direct transitions from p2 to pMN and from pMN to MNDiff. The figure shows 10, 500 nM SAG (see Fig A14 A in S1 Appendix for other concentrations). (Top) At 10 nM SAG, the p2 AC is relatively stable but shallow as cells can escape and transition to pMN over time. (Bottom) The same transition at 500 nM SAG where p2 has already bifurcated showing more rapid occupation of the pMN AC. The data underlying this Figure can be found at https://crick.figshare.com/projects/Neural_Tube_Decision_Landscape/250460 and at https://doi.org/10.5281/zenodo.15584010.

https://doi.org/10.1371/journal.pbio.3003953.g006

Fig 7. Sub-landscapes identification.

A. Schematic diagram of all the possible transitions between cell states. A sub-landscape is defined as a group of neighbouring ACs and transitioning cells, projected in LDA space and highlighted on the diagram. A mathematical model should reproduce all the possible transitions. B. Distinct dynamical systems that model the transitions in the intermediate sublandscape. In Panel 1, p2 and pMN are close to their saddle and cells transition easily towards MNDiff. In Panel 2, the unstable manifold is directed towards DP and cells transition to DP via pMN. In Panel 3, the pMN attractor is bifurcated and DP is close to its right saddle, allowing a DP-to-p3 transition. C. The ACs MNDiff, pMN, DP and p3 form a smooth structure in LDA space, supporting our choice of using a linear landscape as in B to model the transitions from pMN to either MNDiff or DP. The data underlying this Figure can be found at https://doi.org/10.5281/zenodo.15584010.

https://doi.org/10.1371/journal.pbio.3003953.g007

Connecting the sub-landscapes to obtain a global model.

The normal form models for the sub-landscapes are then connected using the simple topological relationships deduced from the data, with connections occurring along transition routes that avoid attractors and bifurcation-critical regions (Figs 9 centre and 10A). This design enables the sub-landscapes to be joined using ideas from differential topology [66] (Section B2.4 in S2 Appendix) for joining local systems together to produce a global one. The result is a dynamical system that depends smoothly on its parameters (Section B2.5 in S2 Appendix). The model qualitatively captures the key bifurcations observed in the experimental data, including fold and flip bifurcations.

Dynamical stochasticity and cell-dependent dynamics.

To model the stochastic nature of the dynamics due to intrinsic noise arising from the randomness of biochemical reactions within a cell we include a spatially dependent diffusion term setting the noise level as an additional parameter. We allowed the noise to vary between the various regions of phase space, accounting for transition-specific noise variability (Section B3.1 in S2 Appendix). Finally, velocity parameters are included that allow the speed along trajectories to vary from region to region. This provides a highly flexible but minimally parameterised model that can accurately capture the behaviour in each local sub-landscape.

As mentioned above when discussing microheterogeneity, nongenetic variability, arising from within-cell processes such as differences in cell size, ribosome number etc, play a role in determining cell fates. To account for this nongenetic variability, we allow for variation in each cell’s parameters producing a distinct dynamic for each cell (Fig 10B, 10C).

When we combine these various types of stochasticity we observe dynamics where cells starting in identical states can traverse widely different developmental trajectories (Fig 10B). Moreover, the timing with which cells follow the same trajectories tends to be highly variable because cells can remain trapped in shallow attractors for extended periods before escaping due to stochastic fluctuations, and then progressing relatively linearly along the unstable manifolds (Fig 10C). We quantify this variation in timing below in Fig 13D when we discuss the DP to p3 transition, showing that it can be very large and represents a challenge to the notion of pseudotime.

Fig 11. Parameter estimation and bifurcation analysis of the dynamic landscape model.

A. Experimental proportions (yellow background) compared to the simulated proportions (blue background) after parameter fitting for each SAG concentration using the flow cytometry data. The experimental proportions were obtained by averaging the proportions for 3 individual experimental series for 10 nM SAG, 100 nM SAG and for 5 experimental series for 0 and 500 nM SAG (S1 Appendix Sect. A3). The predicted FP proportions (12%) obtained from simulations at 500 nM SAG are consistent with those (7%) observed in the scRNA-seq at the same SAG concentration (see B below). B. Similar to A for the scRNA-seq data at 500 nM SAG. C. (Left) Experimental proportions for flow cytometry 100 nM SAG. (Middle) Fitting the model to only the D4 and D6 experimental proportions accurately predicts the intermediate D5 proportions. (Right) Fitting the model to only the D5 experimental proportions predicts the proportions at D4 an D6. D–G. Dots represent parameters accepted by the Approximate Bayesian Computation particle filter. Contour lines are used to estimate the PPDs projected on the corresponding sub-landscape parameter space (blue: 500 nM SAG scRNA-seq, red: 500 nM SAG flow cytometry, orange: 100 nM SAG flow cytometry, yellow: 10 nM SAG flow cytometry, purple: 0 nM flow cytometry). The other curves show the bifurcation set where fold and flip bifurcations occur for the different sub-landscapes. Panels (E–G.) contain the sub-landscape phase portraits that arise for parameters chosen within each region of bifurcation set. D. (Left) p2 sub-landscape bifurcation set. (Right) PreNeural sub-landscape bifurcation set. The flow cytometry and scRNA-seq PPDs at 500 nM SAG fall in the same regions of parameter space, delimited by fold and flip curves. The associated dynamics and decisions are thus qualitatively identical. E. Intermediate sub-landscape bifurcation set with the flow cytometry and scRNA-seq PPDs at 500 nM SAG. The shaded region indicates the set of parameters for which the DP attractor is present. As the blue PPD overlaps the DP fold bifurcations curve, we observe heterogeneity: some cells remain trapped in the DP attractor whereas others rapidly transition to p3. The flow cytometry PPD (red) remains largely within the shaded region, where the DP attractor is present, although the fraction of cells with parameters close to the DP fold curve still transition from DP to p3 due to stochastic fluctuations. F. Bifurcation set of the sub-landscape controlling the p2 transition. The yellow fold curve represents the parameter values where the p2 attractor bifurcates. At 10 nM SAG, the PPD lies on both sides of the p2 fold curve, suggesting that the p2 attractor is present, but shallow, for half of the cells and bifurcated for the other half in which case cells transition to a downstream attractor. At higher SAG concentrations (100, 500 nM), the PPDs lie outside the shaded region suggesting that the p2 attractor is bifurcated for all cells. The 0 nM PPD (purple) is not shown, as the parameters did not converge due to the absence of p2 cells at this SAG concentration (the parameters are distributed across the entire region). G. Bifurcation set of the sub-landscape governing the binary decision from PreNeural to intermediate or ventral fates. Parameters on either side of the grey flip curve yield distinct sub-landscapes: one with an unstable manifold (blue) connecting PreNeural to an intermediate fate (p0/p1, p2, pMN), and another (red) to a ventral fate (EarlyVentral, Early p3, p3, FP). Crossing the purple fold curve causes the PreNeural attractor to bifurcate, with fate determined by the basin of attraction the cell is in. At low SAG concentrations (0–10 nM), PPDs lie left of the flip curve, indicating transitions primarily towards intermediate fates. At higher concentrations (100–500 nM), PPDs are near the fold–flip point, allowing heterogeneous transitions to both intermediate and ventral fates. The data underlying this Figure can be found at https://crick.figshare.com/projects/Neural_Tube_Decision_Landscape/250460 and at https://doi.org/10.5281/zenodo.15584010.

https://doi.org/10.1371/journal.pbio.3003953.g011

Comparison between flow cytometry and scRNA-seq.

Given longstanding questions about the relationship between expression levels of mRNA and protein measurements, it was unclear whether model parameters fitted to both types of data would agree. We therefore compared the PPDs obtained by fitting the model to the scRNA-seq data with those of the flow cytometry at 500 nM SAG by projecting them onto the bifurcation sets for each sub-landscape to reveal how signalling conditions control cell fate. We define two PPDs as qualitatively the same if they occupy the same region of the bifurcation set delimited by fold and flip curves. For most sub-landscapes, the SAG 500 nM PPDs were qualitatively the same (Fig 11D), This agreement additionally validates the robustness of the landscape structure across measurement modalities.

There was, however, one notable difference. The scRNA-seq PPD crossed the DP fold curve while the cytometry PPD did not (Fig 11E), forcing DP cells to rapidly transition to p3. We hypothesised that this discrepancy reflected the substantially increased number of p3 cells in the scRNA-seq data on D7 (Fig 11B), as these include the V3 neurones that become abundant at D7. In comparison, the flow cytometry data extended only to D6 and the numbers of p3 cells were lower and did not include V3 neurones as only Sox2-expressing cells were selected. Simulations revealed that the substantial p3 population at late time points cannot be achieved solely through the PreNeural-to-p3 pathway via the ventral sub-landscape. This suggested that the PPD crossing the DP fold curve (Fig 11E) reflects an alternative route: direct transition from DP to p3. The model predicts that cells can reach the p3 fate through two distinct developmental pathways, a prediction we validate below.

Bifurcation sets determine how molecular heterogeneity translates into decision heterogeneity.

As discussed in Box 3, bifurcation sets fundamentally shape how molecular heterogeneity affects cell fate outcomes. When noise levels are low, two cells, whose landscape parameters are in the same region of the bifurcation set, will make the same decisions with high probability. Conversely, nearly identical cells can make reliably different fate choices if a bifurcation curve divides them into distinct regions. This implies that heterogeneity in differentiation outcomes among molecularly similar cells arises from a cell population straddling bifurcation boundaries rather than from intrinsic noise (Fig 1B). Understanding bifurcation set structure and how estimated parameters map onto it is therefore essential for interpreting decision heterogeneity in single-cell data. We analysed the four PPDs obtained by fitting the model to flow cytometry data at constant SAG concentrations.

The bifurcation set controlling the direct transition from p2 contains a fold bifurcation curve (Fig 11F). Crossing this curve destabilises the p2 attractor, forcing a cell to transition to pMN or DP. At 10 nM SAG, parameters fall on both sides of the fold curve, creating heterogeneity where some cells remain in p2 while others transition. This reproduces the bistability in the p2 and the pMN cells population observed at D6 in Fig 6H. At higher concentrations (100, 500 nM), the PPDs lie below the fold curve, indicating that the p2 attractor has bifurcated for all cells and driving uniform differentiation to pMN. Even at low SAG concentration (10 nM), this in vitro protocol fails to stabilise the p2 state. The model predicts that producing a persistent population of p2 cells is not possible by modulating Shh signalling alone.

In the PreNeural sub-landscape bifurcation set, a flip and a fold curve intersect at a fold-flip point (Fig 11G). At low SAG (0, 10 nM), the PPDs localise to the left of the flip bifurcation curve, biasing cells to intermediate fates. At high SAG (100, 500 nM), the PPDs concentrate near the fold-flip point, generating heterogeneous outcomes between the two fates.

Other informative examples showing the relationship between bifurcation and PPDs mappings and validating agreement between scRNA-seq and flow cytometry data are discussed in S2 Appendix Sect. B5.1 and Sect. B7.3).

Overall, the fitted model accurately captured the experimental dynamics (S1–S3 Videos), and successfully reproduced the observed bifurcation behaviours. This framework revealed how individual bifurcations combine to generate complex decision structures and provides a qualitative and quantitative explanation for Shh signalling dependent patterning of ventral neural tissue and makes testable predictions.

Model predictions of signalling induced landscape changes.

We used the fitted model to develop and test various hypotheses. We first asked whether temporal changes in Shh signalling resulted in instantaneous or delayed landscape reconfigurations. We focussed on three principal decision points: the transition from PreNeural to either the intermediate or ventral sub-landscapes; the transition from p0/p1 to p2 that leads to pMN; and the connection from DP to p3 responsible for the circular topology of the landscape.

To examine the first two decisions, we recast the model into a simplified version that focussed on these transitions (Fig 12A). This reduced model aggregated several related cell states (combining NMP/PreNeural, EarlyVentral/Early p3/p3, and pMN/MN/DP) while preserving the essential topology of a binary flip landscape followed by a linear transition pathway (S2 Appendix Sect. B8). In this framework, the ventral sub-landscape consists of a single state p3*, and the intermediate sub-landscape comprises two states: p2 and the combined pMN/MNDiff/DP state, relabelled pMN*. The model accurately reproduced the proportions of cell states observed in various SAG concentrations (Fig 12C). The full global landscape model produces equivalent predictions for these signalling perturbation experiments; details are provided in S2 Appendix Sect. B6. The simplified version is presented here.

Timed pulses. We used the model to simulate the effect of timed pulses of Shh signalling on neural progenitor proportions. First, we simulated the effect of 24 h with no SAG followed by 48 h with 500 nM SAG. Under the hypothesis of instantaneous landscape responses, the increase in SAG would trigger immediate bifurcations, destabilising dorsal states (p0-p2). This predicted that, at D4, after 24 h with 0 nM SAG, cells would adopt a p0/p1 state, at which point an abrupt increase to 500 nM SAG would destabilise this state, driving the cells to transition further into the intermediate sub-landscape and not into the usual ventral route to p3. Since the p2 state is absent at 500 nM SAG, cells would rapidly adopt a pMN* identity. Importantly, because most cells are predicted to have left the PreNeural attractor and committed to an intermediate fate by D4, few PreNeural cells would be left to transition to the ventral sub-landscape. Consequently, under this hypothesis, the model predicts substantially fewer p3* cells and a purer population of pMN* at D5-D6 compared to constant 500 nM SAG exposure.

We experimentally tested this hypothesis by culturing cells in the absence of SAG for 24 h, from D3 to D4, and then adding 500 nM SAG for the subsequent 48 h (D4-D6). We analysed the proportions of cell types at D4, D5, and D6 using flow cytometry. The experimental data were in good agreement with the simulations, with the majority of cells adopting a pMN* identity at D6 (Fig 12D). These results indicated that the response to raising levels of SAG resulted in fast predictable changes to the landscape.

We next simulated the opposite pulsing regime—500 nM SAG for 24 h followed by 0 nM SAG for 48 h. to explore the impact on the landscape of an abrupt removal of SAG at D4. Assuming cells respond immediately to the removal of the signal, triggering instantaneous bifurcations in the landscape, the model predicts the reestablishment of the p0/p1 state and cells that had not reached the p2 state by D4 would remain trapped in p0/p1. Consequently, we would expect a higher proportion of p0/p1 cells to persist at D6 compared to continuous exposure to 500 nM SAG.

Experimental results showed a markedly different outcome from this prediction: pMN* cells dominated at D5 and D6 and p0/p1 cells were largely absent as in the constant exposure case (Fig 12E). The almost complete absence of p0/p1 indicated an irreversible commitment to pMN* fates and that the cells did not respond to an instantaneous signal withdrawal. This is consistent with previous experimental observations of bistability and hysteresis in neural progenitor specification [72].

The similarity in cell states proportions between constant 500 nM exposure and SAG removal after 24 h rules out the assumption of an instantaneous cellular response to signal withdrawal. To distinguish whether the response to withdrawing SAG is delayed or entirely absent, we fitted the model directly to the experimental data (Fig 12E) and compared the resulting PPD ‘post-500 nM-SAG’ with that from continuous 500 nM SAG exposure. The PPDs were identical (Fig 12B), indicating that the landscape remains unchanged after SAG withdrawal and the dorsal states (p0-p2) are not reestablished. Similar results from 100 to 0 nM SAG transitions demonstrate that 100 nM SAG is also sufficient for irreversible p0/p1 destabilisation (Fig 12F). Further analysis, including an experiment using a shorter 12 h pulse of 500 nM SAG, confirmed that the sustained landscape effect persists throughout the post-withdrawal period (S2 Appendix Sect. B6).

Experimental validation of convergent topology. Finally, the global model predicts that, at 500 nM SAG, a substantial proportion of cells reach p3 via the intermediate DP route, arriving later than cells taking the direct ventral pathway. Consistent with this prediction, scRNA-seq data showed concurrent DP depletion and p3 accumulation between D6 and D7 (Fig 11B).

To test this prediction, we simulated two signalling regimes: constant 500 nM SAG versus 24 h at 0 nM followed by 72 h at 500 nM SAG (Fig 13A, 13B). We tracked which route cells followed to reach p3. Under the second signalling regime, the vast majority of p3 cells took the intermediate DP-to-p3 route rather than the ventral route via EarlyVentral and Early p3 (Fig 13C). Entry time analysis revealed distinct kinetics for the two pathways. Cells reached p3 earlier via the ventral route, whereas the intermediate route resulted in delayed arrival due to retention in shallow intermediate states (Fig 13D).

To test these predictions, we performed experiments extending the time course to D7 (Fig 13E, 13F) and analysed the effect of adding 500 nM SAG at D4 after 24 h without SAG. Under these conditions, the model predicts that PreNeural cells are forced to take the intermediate route as almost all cells initially transition to p0/p1. Subsequently, elevating SAG levels to 500 nM SAG destabilises the p0/p1 and p2 attractors, making p3 accessible via the DP saddle. Flow cytometry data at D7 revealed that cells had transitioned from DP to p3, characterised by decreased Olig2 expression while maintaining Nkx2.2 (Fig 13F). This is particularly clear at D7. In conditions in which SAG exposure is delayed by 24 h, p3 cells emerge despite the near absence of Early p3 at earlier time points, ruling out the possibility of a ventral route.

We compared the simulated proportions of DP, Early p3 and p3 ACs for both experimental conditions (Fig 13G, 13H), using the PPDs previously obtained for 0 and 500 nM (this experiment was not included in the fitting of the model). The model captured the increased proportion of DP cells and the passage from DP to p3 via the intermediate route (Fig 13H) (S3–S4 Videos).

Together, these results further validate the model and demonstrate the complex relationship between signal dynamics and cell fate allocation. The analysis indicates that transient high-level Shh signalling induces irreversible commitment to ventral fates, while also providing support for the dual pathways to p3 identity. These findings establish how cells integrate and retain signalling information to make robust fate decisions, with the mathematical model providing a framework for understanding the underlying dynamical principles.

Fig 14. Analysis of human neural progenitor development validates the circular topology of the decision landscape.

A. Summary of the outward clustering methodology for the human organoid data. B. The most significantly differentially expressed genes found using the 4 samples defined in A. C. The four cell type domains including neural progenitors were extracted, annotated and separately projected using LDA from the effective gene space defined by the gene modules in B. D. 3D LDA projection of the neural progenitor ACs for all time points and experimental conditions. This reveals the circular topology of the transition pathways. E. Simulated neural progenitor ACs produced by the landscape model of this system after fitting in a similar fashion to that described above. Dots represent simulated cells retained at 3 time points corresponding to D3, D5, D7 using the parameters from the fitted PPDs for the 18h-Delay condition. Grey dots represent transitioning cells that were not assigned to any simulated AC. F. Dot plot indicating the expression of markers used to identify the tissue domains. The intensity of shading indicates the mean expression of each gene in each domain. G. Proportions of cells in each domain of C and F by time point in the two experimental conditions. H. Dot plot indicating the expression of markers used to identify the ACs and transitioning clusters of the Neural Progenitors domain in D. I, J. The PreNeural sub-landscape extracted from the human dataset (I) and the mouse datasets (J) and transition analysis for marker genes that characterise the binary decision as cells exit PreNeural. Unstable manifolds defining routes between pairs of ACs were computed in effective gene space defined by the genes in B and projected onto 2D LDA space. Marker genes were analysed along both routes. Exit from PreNeural is marked by down-regulation of Cdx2 and Nkx1.2. A branching event occurs (around index 10), followed by up-regulation of markers characteristic of ventral states along the red route (leading to EarlyVentral), and markers such as Pax6 characteristic of the intermediate states p0/p1, and p2 along the blue route. K. (Top) Comparison of the experimental and simulated proportions after parameter fitting for the 18 h-Delay condition. (Bottom) Similar for the 24 h-Delay condition. This shows an increase of ventral states such as p3 and FP. L. (Top) Bifurcation set of the PreNeural sub-landscape shows the qualitative difference between the 18 h- and the 24 h-Delay conditions. The 18h-Delay PPD (blue) appears on the left of the flip curve and above the PreNeural fold curve. Under this condition, the PreNeural attractor is shallow for most cells and is connected to an intermediate fate. By contrast, the 24 h-Delay PPD (red) lies to the right of the flip curve, favouring cell trajectories from PreNeural to a ventral fate and promoting the generation of more p3 and FP cells. (Bottom) Typical configurations of attractors, saddles and unstable manifolds for the cell landscape models for the PreNeural sub-landscape for parameters in each of the PPDs. In the 18 h-Delay, cells are most likely to transition via the intermediate route towards p2 and to Early p3 via the ventral route in the 24 h-Delay. As indicated by the PPDs in M, the attractors p0/p1 and EarlyVentral are bifurcated for all cells in both experimental conditions. Hence, the connections from PreNeural lead directly to p2 (blue unstable manifold) or to Early p3 (red unstable manifold). M. Bifurcation set of the sub-landscape controlling the p0/p1 and EarlyVentral fold bifurcations. Under both conditions, the PPDs are in the region where the p0/p1 and the EarlyVentral attractors are bifurcated for all cells. N. Bifurcation set of the p2 transition shows that, under both conditions, the p2 AC is shallow as most parameters cluster near (and above) the p2 fold curve in the shaded region. The data underlying this Figure can be found at https://doi.org/10.5281/zenodo.15584010.

https://doi.org/10.1371/journal.pbio.3003953.g014

Conservation of the decision landscape topology.

Despite differences in developmental timing, gene expression patterns characterising the ACs (Fig 14H) and individual analysis of each sub-landscape showed striking conservation between mouse and human systems (Fig 14I, 14J for the PreNeural sub-landscape and Figs A16-A17 in S1 Appendix for all the other sub-landscapes). We therefore fitted our landscape model to the two datasets 18h-delay and 24h-delay focussing on the neural progenitor landscape. We anticipated that, since the 18h-delay condition contains less notochord and therefore less Shh signalling than the 24h-delay condition, the 18h-delay dataset would behave as if it had been exposed to lower SAG concentrations than the 24h-delay dataset. Cells in the D3 PreNeural AC were used as initial conditions and comparisons were made to the proportions of cells in each simulated AC with the experimental cell states proportions. The model successfully reproduced the experimental data (Fig 14K). Details of the fitting and the simulation can be found Section B9 in S2 Appendix.

Using bifurcations to quantify the effect of notochord-derived Shh on the decision landscape.

Analysis of the bifurcation sets revealed key differences between the 18h-delay and 24h-delay conditions. The 24h-delay generated more notochordal cells (Fig 14G), resulting in higher proportions of floor plate and p3 progenitors and fewer pMN cells by D7 than the 18h-delay. Importantly, the bifurcation analysis showed that these differences primarily affected the PreNeural flip bifurcation (Fig 14L), while other bifurcations remained unchanged between conditions (Fig 14M, 14N). In both experimental conditions, the p0/p1 attractors were unstable across all cell landscapes, making the p0/p1 AC transient - mirroring the behaviour observed at high SAG concentrations in the mouse experiments (Fig 14M). Most p2 attractors remained shallow, with the corresponding PPDs clustering just within the p2 bifurcation fold curve (Fig 14N). This suggests that the increased number of notochordal cells in the 24-delay condition, and consequently the higher level of Shh signalling, primarily influenced the first neural progenitor cell fate decision, directing more cells towards the ventral sub-landscape, providing quantitative insight into how notochord-derived Shh signalling affects ventral patterning.

These findings demonstrate that the landscape model successfully captures neural patterning dynamics across species while revealing features of human development. The analysis identified how cell state transitions are affected by notochord-derived Shh: while the PreNeural flip bifurcation showed high sensitivity, other transitions displayed qualitative robustness to signalling variations. This raises important questions about how dorsal fates are established in the presence of strong ventralising signals from the notochord, suggesting the necessity of additional regulatory mechanisms. Future studies will need to explore how combinations of signals shape the early neural tube landscape.