(C) PLOS One
This story was originally published by PLOS One and is unaltered.
 . . . . . . . . . .



Dynamic Landscape Analysis of cell fate decisions provides predictive models of neural development from single-cell data [1]

['Marine Fontaine', 'Mathematics Institute', 'University Of Warwick', 'Coventry', 'United Kingdom', 'Zeeman Institute For Systems Biology', 'Infectious Epidemiology Research', 'M. Joaquina Delás', 'The Francis Crick Institute', 'London']

Date: 2026-09 

We estimated model parameters for each SAG concentration using Approximate Bayesian Computation particle filtering [15,18,67–71] which yields a probability distribution over plausible parameter values rather than a single best-fit estimate. This involves simulating each cell’s trajectory using the stochastic model on a finite time interval, calculating the proportions of simulated cells forming a simulated AC at regular intervals, and comparing these with their experimental counterpart (Fig 10C). Parameters are iteratively sampled until convergence to a posterior parameter distribution (PPD) that optimally reproduces the experimental cell state proportions. Details are provided in S2 Appendix, Sect. B4 for the general fitting procedure, Sect. B5 for the fitting to the flow cytometry data, and Sect. B7 for the scRNA-seq data.

Fitting the model to the flow cytometry data for all SAG levels (Fig 11A) and the scRNA-seq data for SAG 500 nM (Fig 11B) gave good agreement with the experimental proportions of the ventral progenitors. FP cells could not be clearly distinguished in the flow cytometry data but the model provided an estimate (12%) of FP proportions at D6 for SAG 500 nM which is consistent with those observed in the scRNA-seq data (7%) in Fig 11B. Additional validation using flow cytometry is given Fig B9 B in S2 Appendix. Importantly, we observed that fitting on a limited set of time points still enabled accurate prediction of cell state proportions at intermediate time points (Fig 11C).

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 .

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.

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.

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.

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.

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.

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.

PNG larger image

TIFF original image Download: Fig 12. Experimental validation of model predictions by signalling perturbations. A. Simulated data showing the attractors and transition routes. Dots represent simulated cells retained at 3 time points (D4, D5, D6). These are aggregated into simulated ACs coloured by the cell state to which they correspond. B. PPDs from model fitting for the 500 nM SAG condition (blue) and the post-high-SAG condition (red) are close to identical, revealing that the two conditions are approximately equivalent. (Left) Projection of the PPDs on the bifurcation set of the binary flip sub-landscape. (Right) Projection of the PPDs on the bifurcation set of the linear sub-landscape. C. Comparison of average simulation results (blue background) for the constant SAG conditions with the corresponding experimental results (yellow background). Results are average over 10,000 simulated proportions using parameters drawn from the fitted PPDs (S2 Appendix Sect. B8). D. As C but for the Up-SAG0/500 condition. The prediction aligns with the experimental data. It shows that delaying SAG application by 24 h commits the cells to differentiate into pMN* rather than p3*. E. As C but for the Down-SAG500/0 condition, modified as follows: (middle) uses the fitted PPDs (500 and 0 nM SAG) obtained from C and assume instantaneous response to signal change, (right) uses the refitted PPD. The latter shows better agreement with the experimental results demonstrating a lack of response to the change of signal. F. As for the first two panels of E for the Down-SAG100/0 condition where the fitted PPDs (100 and 0 nM SAG) obtained from C are used. We observe the same phenomenon as in E (middle panel), despite lower levels of the signal. The data underlying this Figure can be found at https://crick.figshare.com/projects/Neural_Tube_Decision_Landscape/250460. https://doi.org/10.1371/journal.pbio.3003953.g012

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

PNG larger image

TIFF original image Download: Fig 13. Validation of the DP-to-p3 transition. A. Simulated cell trajectories using the PPD from model fitting (500 nM SAG condition). Dots represent simulated cells retained at 4 time points (D4-D7). Cells reaching the p3 AC by D7 are coloured based on their transition history: yellow for cells that follow the intermediate route and green for those that follow the ventral route. Black dots represent cells that adopted alternative end fates (such as MNDiff or FP) or have not yet reached p3 by D7. B. Similar to A, but simulating a delayed pulse of SAG by switching the parameters from 0 to 500 nM SAG after D4. Substantially fewer cells follow the ventral route as most cells had already taken the intermediate route by D4 when the signal changed. C. Proportions of simulated cells in the p3 AC for each day (D5 to D7) computed with respect to the total number of p3 cells at D7 and coloured by the route they originated from. (Left) 500 nM SAG. (Right) Up-SAG0/500 nM SAG condition. The majority of p3 cells under this condition originate from the intermediate route. D. Histograms showing cell counts (y-axis) by the time point (x-axis) at which simulated cells first reach the p3 AC, with colours indicating their route of origin: intermediate (yellow) or ventral route (green). The distributions reflect timing heterogeneity: cells that reach p3 via the intermediate route tend to arrive at late time points (D6-D7). By contrast, cells following the ventral route tend to reach p3 early (D4-D5) while few new cells arrive via this route after D6. E, F. 2D Gene plots (Olig2-Nkx2.2 shown) of the clusters identified in the flow cytometry dataset at each time point. The dots represent individual cells, and contour lines represent the density of cells (plotted only when the cluster contains more than 150 cells). Without Nkx6.1 p0/p1 cannot be distinguished from p2 cells and both are combined into Pax6+ cells. PreNeural, EarlyVentral and FP are combined into Sox2+ cells. E. 500 nM SAG: as early as D4, a group of cells follow the ventral route by up-regulating Nkx2.2 and reach p3 by D5. F. Up-SAG0/500 nM SAG condition: with no SAG for 24 h, cells are forced to take the intermediate route and most of them had up-regulated Pax6 by D4. When SAG is applied, they can only reach p3 by transiting through pMN (D5), then DP (D6), followed by a direct DP-to-p3 transition (D7). G. (Left) 500 nM SAG: experimental proportions of the DP AC and p3Trans (intermediate route) and the Early p3 and p3 ACs (ventral route) at late time points (D5-D7). (Right) 500 nM SAG: simulated proportions of the DP, Early p3 and p3 ACs at late time points (D5-D7). The p3 proportions can be further split into two type of cells depending on which route they followed. H. Same as G for the Up-SAG0/500 nM condition. In the simulation (right), most p3 cells at D7 originate from the intermediate route as shown in B which is consistent with the experiment F. 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.g013

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.

[END]
---
[1] Url: https://journals.plos.org/plosbiology/article?id=10.1371/journal.pbio.3003953

Published and (C) by PLOS One
Content appears here under this condition or license: Creative Commons - Attribution BY 4.0.

via Magical.Fish Gopher News Feeds:
gopher://magical.fish/1/feeds/news/plosone/