Quantitative, real-time, single cell analysis in tissue reveals expression dynamics of neurogenesis

Despite single cell transcriptomics, how cells make transitions within tissues in real time, is not understood. Here, we use single cell live imaging of an endogenous HES5 reporter and absolute quantification to gain a dynamic view of neurogenesis in the embryonic mammalian spinal cord. We show that HES5 protein expression fluctuates in dividing neural progenitors and becomes more frequently periodic in the transition towards differentiation, creating transient oscillations with higher fold-changes. This dynamic behavior means that the HES5 population heterogeneity observed between cells at a fixed time-point, is a composite of short-term and longer-term dynamics.


Introduction
During embryogenesis cells balance proliferation with differentiation to make cell state transitions that lead to the formation of functional organs. This is exemplified by the development of the central nervous system, which requires the balance of neural progenitor maintenance with differentiation during multiple waves of differentiation in to neuronal and glial cell-types 1 . In the dorso-ventral (D-V) axis of the spinal cord elegant experiments have shown that fate decisions require integration of a wide range of signals over time, many in the form of morphogen gradients, resulting in downstream gene expression changes 2,3 . Single-cell transcriptomics have greatly enhanced our understanding of these gene expression changes, the gene networks involved in fate decisions and of the bifurcation points where decisions are made [4][5][6][7] .
Transcriptomic studies are powerful in revealing cohorts of up-regulated or downregulated genes and in defining sub-states and branching points. However advances in single-cell live imaging of gene expression has shown that it is often highly dynamic, suggesting that the control of cell state transitions is more complex [8][9][10] . Rather than being in an "on" or "off" state, a handful of transcription factors have been shown to oscillate with periodicity of a few hours 9,11 . While such oscillations have been long described in somitogenesis 12 , in the development of neural tissues they are a relatively recent discovery. This is because unlike somitogenesis where oscillations are synchronous within each somite, they tend to be asynchronous in neural cells and so required unstable reporters and single cell imaging to be discovered 13 . Thus, it is not only changes in gene expression levels that are important, but the short term dynamics of gene expression can also carry important information for cell state transitions. Indeed, there is experimental and theoretical evidence that cell fate transitions may be controlled by a change in the dynamic pattern of gene expression, which could be from oscillatory to stable expression, or to oscillatory with different characteristics 9,14,15 . In the case of the transcriptional repressor HES1, a key target of Notch signaling, it has been known that oscillatory expression is driven by transcriptional autorepression coupled with delays, instability of mRNA and protein and non-linearity of reactions, common principles of many biological oscillators 16,17 . Like HES1, HES5 is a Notch target bHLH transcription factor (TF) which is highly expressed by neural progenitor cells and decreases in expression as differentiation proceeds 18,19 . Knockout mice and over-expression studies have shown that HES5 functions to maintain the undifferentiated progenitor state through repression of proneural genes, such as Neurog2 and Atoh1 that promote neuronal differentiation [20][21][22] . Like HES1, HES5 has been shown to oscillate in neural progenitors in vitro 9 .
More recently, it was shown that a change in HES dynamics is mediated by a change of the parameters or initial conditions of the oscillator, most likely represented by a change in mRNA stability or protein translation under the influence of a microRNA, miR-9 [23][24][25] . Other theoretical studies provide additional support for the importance of a change in dynamics by showing that gene expression networks in the D-V dimension of the spinal cord can generate multi-way switches (stable or oscillatory) 26 .
An additional revelation of single-cell live imaging studies is that gene expression is characterised by varying degrees of noise due to the stochastic nature of transcription [27][28][29] . Current ideas for the role of such embedded stochasticity include cases where it would be an advantage 30,31 or conversely, an impediment for cell fate decisions 32,33 and mechanisms to suppress noise after a fate-decision 34 .
However, although these studies have shed new light into the problem of cell-state transition, how cells make decision in the context of multicellular tissue is poorly understood. This is because both single-cell transcriptomics and live imaging data are routinely performed in single cells taken out of the tissue environment. Existing studies of oscillatory expression in the mouse brain and spinal cord lack the statistical power needed to give a comprehensive understanding of the dynamics in the tissue 11,35 . A study using electroporation of a promoter reporter of Hes5-1 in chicken spinal cord tissue reported activation of Notch signaling throughout the progenitor cell cycle but most frequently before mitosis 36 . This approach suffered from plasmid loss and varying degrees of plasmid transfection and did not report on endogenous HES5.
Here, we have taken advantage of a Venus::HES5 knock-in reporter to study the dynamics of HES5 in the context of a tissue, with single cell resolution, over time. We have developed an ex-vivo organotypic slice system of embryonic mouse spinal cord (E10.5) where HES5 dynamics can be continuously imaged over at least 12 hours. We have also developed statistical tools to analyse the behaviour of hundreds of HES5 expressing cells over time. Coupled with absolute quantitation of the HES5 fluorescent fusion protein, we were able to create a finely resolved map that shows precisely the heterogeneity that exists in the tissue and furthermore gives insight into how it is generated.
We report that HES5 expression has a 10-fold range between cells in a single expression domain that can be accounted for by short-term fluctuations and longerterm trends of decreasing HES5. We use hierarchical clustering to define distinct clusters of HES5 expression dynamics and then use cell position and division properties to infer cell state. Surprisingly, we find that most of the oscillatory behaviour is observed in cells that transition towards differentiation where it is coupled with an overall decrease in HES5 expression. By contrast, dividing neural progenitor cells are less frequently periodic but significantly more noisy in their HES5 expression. Computational modelling with stochastic differential delay equations, parameterised using experimental values and Bayesian inference, suggest that in the spinal cord tissue environment the Hes5 genetic oscillator operates close to a bifurcation point where noise can tip it from aperiodic to periodic expression. Taken together, our findings suggest that single progenitor cells in a tissue environment are "noisy" and are thus primed to enter a transient oscillatory phase as the cells differentiate. Our work reveals for the first time the single-cell dynamics during cell state transitions in a tissue context. Additionally, our study reveals that tissue level single-cell heterogeneity has a complex origin in both short and long term dynamics.

Venus::HES5 reporter in the ventral embryonic spinal cord recapitulates endogenous features
We characterised the Venus::HES5 knock-in mouse 9 to ensure that it is a faithful reporter of the un-tagged gene. In transverse sections of the spinal cord at E10.5 Venus::HES5 showed a broad ventral and a smaller dorsal domain (Fig.1a). The ventral domain, which is the focus of this study, encompasses mainly ventral interneuron (p0/p2) and some ventral motorneuron progenitors (p2/pMN) as shown by mapping to appropriate regional (Supp. Fig. 1a and b) and neural/neuronal markers (Sox1/2+ progenitors, NeuN+ neurons, Fig 1b) and consistent with reports of endogenous HES5 7 .
Both mRNA and protein half-lives of Venus::HES5 are unstable with similar values to untagged HES5 (approximately 30 mins for the mRNA and 80-90 mins for the protein). These findings confirm that the Venus::HES5 fusion protein is a faithful reporter of endogenous un-tagged HES5 expression (Supp. Fig.1c-f).

Direct imaging and fluorescence correlation spectroscopy (FCS) of Venus::HES5 protein quantifies the range and level of HES5 expression in single cells
Dynamic expression can give rise to tissue level single-cell heterogeneity which may be masked by population averaging on non-quantitative methods. Here, we have imaged nuclear Venus::HES5 intensity directly, to avoid signal amplification associated with immunohistochemistry, and characterised the true heterogeneity in expression within the ventral domain in tissue slices (Fig.1b). We then obtained absolute quantitation of Venus::HES5 molecules at the single cell level by FCS of neural progenitor cells in live homozygous Venus::HES5 E10.5 embryo slices (Fig.1c,d Supp. Fig.2a,b,c). FCS is an absolute quantification method that records fluorescence fluctuations as molecules diffuse through the confocal volume 37 . Therefore, it provides a means of distinguishing true signal from auto-fluorescence and allows for a comparison between different experiments. Single cells showed a 10-fold range of nuclear Venus::HES5 protein expression within the ventral Venus::HES5 expression domain, from 26nM to 319nM. (Fig.1d). The mean Venus::HES5 nuclear concentration was calculated as 140nM, or 46,250 molecules per nucleus. Heterozygous embryos showed lower mean protein expression as could be expected by monitoring the expression of one allele (Supp. Fig.2d) but similar variability between cells. These findings show a high degree of variability in Venus::HES5 expression between cells which is similar in homozygous and heterozygous embryos suggesting that integrating the expression from 2 alleles does not diminish the variability that cells experience.

Absolute protein quantitation reveals spatial patterns of single-cell heterogeneity throughout the HES5 expression domain
FCS can be performed for a limited number of live cells in the tissue, while an intensity map based on the Venus signal can be obtained for all cells from snapshot images. By combining the two approaches we were able to obtain a quantitative map of all expressing cells along the D/V axis of the spinal cord 38 . We plotted the distribution of single-cell Venus::HES5 intensities from manual segmentation of nuclei in a single slice (Fig.1g) against the distribution of single-cell FCS protein concentration (Fig.1h) over multiple slices and experiments. The resulting quantilequantile (Q-Q) plot is linear and only deviates from linearity at the very high and low values (Fig.1i). Therefore we can turn intensity in an image into protein concentration (Fig.1j) by scaling the intensity value by the gradient of the linear Q-Q plot. Once the Venus::HES5 protein concentration distribution has been obtained by FCS it can be applied to multiple images to generate more quantitative maps without the need to repeat the FCS (see Supp. Fig. 2e for another example).
We then used this quantitative map to look at the global and local pattern of HES5 concentration. We split the ventral domain into 2 regions due to the difference in the width of the ventricular zone along the D-V axis (as indicated by boxes in Fig.1j) We observed a non-linear global reduction of Venus::HES5 concentration with increasing distance from the ventricle (Fig. 1k). The shoulder-point corresponded to around 50µm and 30µm in the dorsal-most (1) and ventral-most (2) regions respectively, suggesting that at this distance, cells start to decrease HES5. At any given distance there is a large degree of cell-to-cell variability in Venus::HES5 concentration. For example, within the first 20µm (roughly 2 nuclei lengths) of the ventricular zone in region 1, some cells have as low as 130nM of Venus::HES5 while others have as high as 250 nM of HES5 expression, a 1.9-fold difference. The concentration difference between a cell and its nearest neighbour (Supp. Fig.2f) increased further away from the ventricle, reaching a maximum of 191nM, a 4.5-fold difference (Fig.  1l). This trend was confirmed in embryos that had not undergone intensity:concentration scaling (Fig.1m). Thus, a global reduction of Venus::HES5 expression further from the ventricle is accompanied by a high degree of fine-grained heterogeneity that increases with increasing distance from the ventricle.

Hierarchical clustering reveals distinct Venus::HES5 expression dynamics
Single cell expression heterogeneity, as described above, may be the result of multiple possibilities: i) fluctuating expression alone (Fig.2ai), which could be periodic and asynchronous ii) distinct but stable cell-state subpopulations (Fig.2aii) or iii) an expression decline as cells transition from one stable state to another (Fig.2aiii). Hypothesis (i) implies HES5 dynamic expression satisfies ergodicity, i.e. variability in a single cell over time can recapitulate the tissue level heterogeneity 39 . To resolve these possibilities for the first time, we performed live imaging of Venus::HES5 expression dynamics at single cell level, in the tissue context of ex-vivo slices in Sox1+ neural progenitor cells (Fig. 2b,c).
We observed multiple types of single-cell Venus::HES5 dynamic behaviours in heterozygous neural progenitor cells (Fig. 2d) over a time period of 12-15 hours. Hierarchical clustering of the standardised Venus::HES5 intensity timeseries suggested 4 clusters of different types of long-term Venus::HES5 expression dynamics ( Fig.2e and Supp. Fig. 5). Cells in cluster 1 and 2 show fluctuating expression around a stable mean whereas cells in clusters 3 and 4 show gradually decreasing and fluctuating HES5 expression (Fig. 2e). The non-standardised mean expression of cells in each cluster maintain this trend (Fig.2f) which is further exemplified by single cell traces (Fig.2g).
The coefficient of variation (C.O.V) of Venus::HES5 over time in single neural progenitor cells (standard deviation of intensity divided by the mean intensity) increased when cells were tracked over 4,8,12,14.25 and 17.25 hours (Fig.2h). By 8-12 hours multiple cells in clusters 3 and 4 had reached similar or higher levels of variation as the variation observed between cells at a single snapshot (Fig.2h) suggesting that declining expression is a major contributor to the tissue heterogeneity. In contrast, cells in clusters 1 and 2 rarely, if ever, reached tissuelevels of variation between cells at a single point in time, suggesting that short-term dynamics have a lesser contribution to overall tissue heterogeneity (and excluding scenario i).
Thus, there exists a limited number of HES5 expression dynamics in neural progenitor cells, represented by clusters, and heterogeneity is generated by a mix of declining expression (long-term trends, scenario Fig.2a iii,) and dynamic fluctuations (short term dynamics, scenario Fig.2a ii,) around a slowly varying mean.

Distinct patterns of Venus::HES5 expression dynamics correlate with cellstates and decision processes.
We hypothesise that the different clusters of Venus::HES5 expression, revealed by hierarchical clustering, may represent different cell-states present in the Sox1+ progenitor population. It is well known that proliferating neural progenitors (Sox1+/2+) are found apically in the ventricular zone, undergo inter-kinetic nuclear migration (INM) and divide at the apical surface of the ventricle 1,40 . Newly born cells fated towards neuronal differentiation migrate basally away from the apical surface in to the mantle zone and exit the cell cycle, turning on markers of differentiation (Tuj1 and NeuN) 40 . We therefore sought to infer the cell state by position, motility and division, recorded from individually tracked nuclei with the ventricle as a reference point. The average position of cells in cluster 1 was significantly closer to the ventricle than those in cluster 3 (Fig.3c). However, the main difference between the positions of the clusters is evident in a zone greater than 50µm from the ventricle, where very few cells of cluster 1 reside and cells in cluster 3 are more abundant (Fig.3d). By contrast, the zone within the first 50µm of the ventricle is equally occupied by cells in clusters 1-4 (Fig.3d).
Nuclei of cells in cluster 1 moved both apically and basally, consistent with interkinetic nuclear migration (INM) but had the shortest displacement as they returned to the ventricular zone. Meanwhile nuclei of cells in cluster 3 and 4 had a larger displacement which was unidirectional towards the basal side ( Fig. 3a,b,e,f and Supp. Fig.7a). The average position of nuclei in cluster 3 tends to be located further away from the ventricle, suggesting they are on their way to differentiation ( Fig.3c and Supp. Fig.7a). Measurement of the size and position of the Sox2+ domain by immunostaining of slice cultures showed that many cells in cluster 3 and 4 moved out from the Sox2+ zone into the mantle zone with concurrent decreasing Venus::HES5 ( Fig.3a,b, Supp. Fig.7a,b,c). In addition, the percentage of cells that divided in the 12-hour monitoring window was significantly higher in cluster 1 and 2; indeed, very few cells in cluster 3 or 4, (if any) were observed to divide (Fig.3g). Given these findings, we inferred that cells in cluster 1 and 2 are proliferating progenitors and cells in cluster 3 and 4 are transitioning towards differentiation.
We confirmed our interpretation of cell-state based on cellular behaviour by treating the spinal cord slice cultures with the Notch inhibitor DBZ to promote differentiation 7 . Spinal cord slice tissue treated with 2µM DBZ showed significantly lower mean Venus::HES5 intensity than control DMSO treated slices (Fig.3h, Supp.Fig. 7d) and an increase in the early neuronal marker B3-tubulin especially in apical regions (Fig.3i). The disorganisation of the neural tube in DBZ treated slices is similar to Hes KO phenotypes 41 . The average position of single cells in DBZ treated slices was further from the ventricle (Supp. Fig. 7e) and they showed significantly increased apico-basal displacement confirming that the Notch inhibition had pushed cells towards basal migration and differentiation (Fig.3j). Hierarchical clustering of standardised Venus::HES5 single-cell intensities showed that 98% of cells in the DBZ treated slices were found in clusters 3 and 4 ( Fig. 3k,l) while the distribution of DMSO Venus::HES5cells recapitulated the presence of all 4 clusters (Supp. Fig.7f,g).
We conclude that cells characterised by a temporally fluctuating Venus::HES5 expression pattern around a high mean (cluster 1 and 2) are proliferating neural progenitors maintained by Notch signalling, while cells with decreasing Venus::HES5 levels over time (clusters 3 and 4) are non-dividing progenitors undergoing cell state transition to differentiation.
6. Cells undergoing state transition to neurons are more frequently oscillatory while proliferative progenitors are more noisy HES5 expression has been reported to be periodic in proliferating mouse neural stem cells derived from the embryonic cortex 9 . Therefore we sought to determine whether any of the clusters showed oscillatory expression. HES5 traces show a high degree of variability in period and amplitude. Detecting oscillatory gene expression in such noisy timeseries, is very challenging and we have previously developed a statistical approach for noisy bioluminescent data, whereby single cell periodicity can be inferred and a statistical confidence level is imposed at population level 42 . Here, we have developed an extension of this method to take into account that fluorescence intensity timeseries from tissue are inherently more noisy (Supp. Methods A). Further, we introduced robust procedures by which globally-informed priors are used to facilitate parameter inference from short timeseries (see Supp. Methods A).
To analyse oscillations, we first corrected for long-term changes in level (trend) caused by HES5 downregulation (Fig. 4a). We then analysed the detrended data with an oscillatory covariance model and inferred the period, amplitude and lengthscale (Fig. 4b). The first two parameters characterise the overall behaviour of the oscillator while the latter accounts for variability in the peaks over time. We compared the oscillatory (alternative) model fit with that obtained from an aperiodic (null) covariance model fit using the log-likelihood ratio (LLR) statistic, which is high for oscillators (Supp. Fig.8b) and low for non-oscillators (Supp. Fig. 8c). Finally, we identified oscillatory cells in each experiment using a strict false-discovery rate criteria set at 3% (Supp. Fig.8e).
We found that overall only 47% of cells showed oscillatory Venus::HES5 expression dynamics (Supp. Fig.8a), while the rest were fluctuating and aperiodic. The estimated period of Venus::HES5 oscillations showed a distribution with a mean of 3.3±0.3 hours (Fig.4d) while H2B::mCherry expression from the ROSA26 locus collected from the same nuclei was aperiodic (Supp. Fig.8a). Surprisingly, oscillations were not restricted to proliferating progenitor cells. Instead, there was a tendency for more progenitor cells on their way to differentiation (clusters 3 and 4) to pass oscillatory tests than dividing progenitors in clusters 1 and 2 (Fig.4c). By contrast, proliferating progenitors in cluster 1 had significantly greater noise than differentiating cells in cluster 3, (more cells with a higher squared-standard deviation of the de-trended Venus::HES5 signal (Fig.4e)). In agreement with this, the likelihood of a cell to have oscillatory Venus::HES5 significantly increased with an increasing average distance from the ventricle (Fig. 4f), whereas noise decreased ( Fig.4g) Given that progenitor cells close to the ventricle (cluster 1 & 2) must turn into the transitory and differentiating cells in cluster 3 and 4, we conclude that progenitor cells have high, dynamic and noisy Venus::HES5 expression which evolves in to a more oscillatory signal as Venus::HES5 decreases and the cells undergo differentiation. Although our observational time window is relatively short, data collected from a few cells in cluster 1 demonstrate this noisy to oscillatory transition in Venus::HES5 expression, supporting this view (Fig.4h, Supp. Fig.9a, Supp. Movie 2).

Mathematical modelling shows the Hes5 oscillator is poised at aperiodic to oscillatory transition.
To understand how the HES5 dynamics of clusters 1 and 2 are generated and how they may transition from aperiodic to periodic expression, we used a stochastic delay differential equation model of an auto-negative feedback network ( Fig.5a and Supp. Methods B) 30,[43][44][45] We parameterized the model using protein and mRNA half-lives (Supp. Fig.1c,d) and Approximate Bayesian Computation (ABC) 46 to infer parameters that are not directly experimentally known (see Supp. Methods B). ABC has advantages over commonly-used point estimates because it provides a probability distribution for estimated parameters thus quantifying parameter uncertainty. We used ABC to search for parameters that give rise to experimentally observed summary statistics of HES5 expression (see Supp. Methods B) and we found that the experimentally measured distribution of oscillation periods and relative standard deviation values in clusters 1 and 2 (Supp. Fig. 10a and Supp. Fig. 10b respectively) are consistent with the predictions from these parameters (Fig. 5b,c).
HES5 expression simulated from inferred parameters can be aperiodic (Fig.5di) or oscillatory (Fig. 5dii,iii) depending on the parameters, as illustrated qualitatively by a sharpening of the peak in the power spectrum and expressed quantitatively by coherence 30 (Supp. Methods B). At unique combinations of parameter values the stochastic model can generate different proportions of aperiodic and oscillatory HES5 expression, across traces and within the same trace. This is consistent with our experimental observations where less than half of cells pass oscillatory tests and we can observe changes in expression dynamics.
We investigated how HES5 expression may transition from aperiodic to oscillatory in two steps. Firstly, we investigated how oscillation coherence varies in response to a single parameter change, in this case the protein degradation rate across parameter space using Bayesian inference ( Fig.5e where each curve corresponds to one possible parameter combination). The experimentally measured protein degradation rate (protein half-life of 90 minutes, blue-line Fig.5e) defines a transition point where the range of possible coherence values changes sharply.
Secondly, we determined the expected coherence in relation to the protein and mRNA degradation rates for the full stochastic model (Fig.5f) and the deterministic model (Fig. 5g). The experimentally measured mRNA and protein degradation rates were located in a region of parameter space where oscillations are expected in the stochastic model, but not in the deterministic model. This is consistent with the results of a full Bayesian comparison between the stochastic and deterministic model (Supp. Methods B) where the likelihood of the deterministic model to describe the HES5 expression statistics is more than 160 times smaller than that of the stochastic model (Supp. Fig.10d). Oscillations were observed in the stochastic system and our experimentally measured degradation rates placed the stochastic system at the boundary of high and low coherence.
Taken together, our modelling suggests that the HES5 oscillator in spinal cord neural progenitor cells is enabled by noise 30,47 and operates very close to the boundary between aperiodic and oscillatory model dynamics, where small parameter changes can cause a transition between non-oscillatory (low coherence) and oscillatory (high coherence) expression.

Oscillations on a downward trend increase the fold-changes experienced by target genes
Detecting oscillations in differentiating cells (cluster 3 and 4) was unexpected based on previous studies. To understand the significance of this finding, we have characterised the dynamic of cells in cluster 3 and 4 in greater detail. We used the Hilbert transform technique for reconstructing instantaneous amplitude (Fig 6a and Supp. Fig 8g) and instantaneous phase (Supp. Fig. 8h) from detrended data. Phase information was used to identify peaks and troughs in the signal to extract fold changes in amplitude (see Supplementary Methods A).
The gradual decrease in mean Venus::HES5 expression in cluster 3 and 4 cells is accompanied by decreasing amplitude of oscillations over time (Fig.6a,b and Supp. Fig 8g) suggesting that Venus::HES5 undergoes amplitude death as the expression is terminated, while the period of oscillators appears to be unaffected (Supp. Fig. 8h). Surprisingly, although amplitude decreases with differentiation, the maximal peak-totrough fold change in Venus::HES5 expression in single cells was significantly higher in differentiating cells in clusters 3 and 4 than proliferating progenitors in cluster 1 (Fig.6c). Furthermore, within cluster 3, oscillatory cells have a higher mean peak-totrough fold change than aperiodically fluctuating cells (Fig.6d,e). Taken together, these findings suggest that oscillations on top of a long-term decreasing signal transiently promote larger fold-changes and thus, may impart greater changes in HES5 molecule number to downstream target genes than either one alone.

Discussion
In this paper, we have investigated how individual Sox1+ neural progenitor cells make cell state transitions in the multi-cellular environment of the embryonic spinal cord. We have also developed statistical tools and computational models to analyse and interpret stochastic dynamics. Our main findings are twofold: firstly, we found that oscillatory expression of HES5 with a periodicity of a few hours is observed in neural progenitor cells in their native tissue environment and occurs more frequently and with higher fold change in cells that are in transition to a differentiated state (Fig.  6f). Secondly, we show that cell-to-cell heterogeneity in HES5 at the tissue level is a composite of long term dynamics (decline in expression) and short term dynamics (fluctuations in a short time scale), the contribution of each was precisely identified using absolute protein quantification.
The first finding supports the existing hypothesis that HES5 oscillates in neural progenitor cells and that HES5 expression declines as cells differentiate. These findings support the view that changes in expression dynamics correlate with transitions in cell state 9 . However, contrary to expectations, we show that only about 30-40% of the dividing progenitors located close to the ventricle show oscillations that pass a statistical test we developed for stochastic oscillators. Our findings also contrast with the current view that only neural progenitors show gene expression oscillations 13 . Instead, we have uncovered a previously unknown state in progenitor cells further away from the ventricle, transitory to differentiation, where oscillations tend to be more frequently observed. By contrast, dividing progenitors that are closer to the ventricle, have a higher level of HES5 expression noise and have a higher proportion of cells that show aperiodic fluctuations. Therefore we observe both oscillatory and non-oscillatory dynamics within 2 defined sub-states -proliferative progenitors and transitory differentiating cells, with the oscillatory dynamics being more likely in differentiating cells.
Computational modelling helped to determine how the dynamics of HES5 expression are generated. A stochastic model of genetic auto-repression recapitulated statistics of the experimentally measured period and maximum amplitude of HES5 protein expression and indicated that HES5 oscillations are enabled by stochastic amplification 47 . According to our mathematical modelling, the HES5 auto-repression network operates near a bifurcation boundary where small changes in parameter values, such as protein degradation rate, can cause cells to switch between aperiodic and periodic expression. Thus the stochastic dynamical model of HES5 supports the existence of transitions between oscillatory and non-oscillatory regimes in the same cells.
In our modeling, we have included the effects of intrinsic stochastic noise, since this choice of noise does not introduce further model parameters and it is inevitably associated with any rate process. Phillips et al. suggested that the low HES1 molecule number leads to stochastic oscillations of HES1 through a finite number effect. By contrast, HES5 molecule number is not low, (approximately 40-60k molecules per nucleus for HES5 versus while 2-3K per nucleus for HES1 30 ). Thus other sources of noise may need to be considered. We hypothesise that noise from stochastic activation of Notch cell-cell signalling in the densely packed tissue as cells undergo INM may actually interfere with the periodicity and coherence of HES5 oscillatory dynamics. An alternative and not mutually exclusive possibility is that additional noise is introduced by cell division 48 or the cell cycle, which would fit our observation that dividing progenitors of cluster 1 have the highest proportion of noisy cells. While there is a view that noise is an undesirable feature of biological systems that needs to be mitigated, there are several cases where noise and the stochasticity that derives from it, may in fact carry some benefit to the decision making process 30 . In the case reported here, the benefit of noise may be to "prime" HES5 expression such that it is poised to become oscillatory.
What is the impact of HES5 periodic and/or aperiodic fluctuations on the control of neurogenesis? The HES5 oscillator operates around a high mean with low peak-totrough amplitude in dividing progenitor cells and it is unclear whether such low amplitude noisy oscillations of HES5 have a role beyond a priming function. The small differences in peak and trough levels may be difficult to differentially decode by downstream genes. Most likely, these oscillations are a by-product of an "active" negative feedback loop that is required for maintaining the HES5 level around a high mean, thus repressing pro-neural genes in most apical progenitors. By contrast, oscillations in the transition to differentiation are coupled with an overall declining trend, and thus generate larger fold differences, which may be easier for downstream targets to decode. This is analogous to a ball bouncing down steps and undergoing greater height drops (oscillatory expresion) than a ball rolling down a ramp (aperiodic expression). Since HES proteins are transcriptional repressors for pro-neural genes, such as Neurog2 and Atoh1 21,49 we predict that the larger fold-changes generated by oscillatory decline in HES5 induces an oscillatory onset of downstream proneural genes 11,35 . We would argue that coupling HES5 oscillations with a declining trend is an ingenious biological way for the cells to be able to decode what is normally a very shallow HES5 oscillator and importantly, to couple it with the process of differentiation.
The second main contribution of this paper is to increase the depth of our understanding of the degree and origin of cellular heterogeneity in gene expression in a tissue environment. While heterogeneity is a common emerging theme during mammalian development, it's origins remains poorly understood. We conclude that HES5 expression in the spinal cord is not an ergodic system since tissue level variability cannot be explained from short-term single cell variability but through a combination of cell sub-states co-existing in the tissue (which can be resolved spatially and dynamically) and transitions between these sub-states. Specifically, hierarchical clustering identified 4 clusters based on their short term and long term dynamics. We have found that a progenitor zone close to the ventricle (<50µm) shows maximum heterogeneity in cell-states, as all 4 dynamic expression clusters are equally represented in this zone, but minimum cell-to-cell heterogeneity in HES5 expression levels. By contrast, in the progenitor zone further from the ventricle there is minimum heterogeneity in cell-states, as it occupied mainly by cells in clusters 3 and 4, and maximum cell-to-cell heterogeneity in HES5 expression levels, approaching a 10-fold range in HES5. Furthermore, single cells undergoing differentiation start to down-regulate Venus::HES5 at any point between 20-50µms away from the ventricle indicating that cells can make the cell fate decision at any point along the apical-basal dimension of the progenitor zone. This contrasts with the widespread schematic view that cell fate is controlled deterministically at global tissue level through signalling gradients. Instead we find that progenitor cells make stochastic fate decisions through a complex and yet unresolved integration between global and local cell-cell signalling.
Our findings highlight the importance of integrating gene expression dynamics with spatio-temporal cell behavior to understand cell state transitions in real time in a multicellular tissue. Challenges for the future include understanding the local spatial cell-to-cell pattern in HES5 expression and extending this type of analysis to observing the expression dynamics of several genes simultaneously, for which a step change in live imaging capabilities/image analysis will be needed.

Animal models
Animal experiments were performed under UK Home Office project licenses (PPL70/8858) within the conditions of the Animal (Scientific Procedures) Act 1986. Animals were only handled by personal license holders. Venus::HES5 knock-in mice (ICR.Cg-Hes5<tm1(venus)Imayo>) 9 were obtained from Riken Biological Resource Centre, Japan and maintained as a homozygous line. In these mice the mVenus fluorescent protein is fused to the N-terminus of endogenous HES5. Sox1Cre:ERT2 mice (Sox1tm3(cre/ERT2)Vep 50 were obtained from James Briscoe with the permission of Robin Lovell-Badge. R26R-H2B::mCherry mice 51 were obtained as frozen embryos from Riken Centre for Life Science Technologies, Japan and C57Bl6 mice were used as surrogates. Sox1Cre:ERT2 and R26R-H2B::mCherry were crossed to generate a double transgenic line homozygous for R26R-H2B::mCherry and heterozygous for Sox1Cre:ERT2.

Embryo slicing
Homozygous Venus::HES5 knock-in females were mated with R26R-H2B::mCherry Sox1Cre:ERT2 males and E0.5 was considered as midday on the day a plug was detected. Intra-peritoneal injection of pregnant females with 2.5 mg Tamoxifen (Sigma) was performed 18 hours prior to embryo dissection. Whole embryos were screened for H2B::mCherry expression using Fluar 10x/0.5 objective on a Zeiss LSM880 confocal microscope and the trunks of positive embryos were embedded in 4% low-gelling temperature agarose (Sigma) containing 5mg/ml glucose (Sigma).
200µm transverse slices of the trunk around the forelimb region were obtained with the Leica VT1000S vibratome and released from the agarose. Embryo and slice manipulation was performed in phenol-red free L-15 media (ThermoFisher Scientific) on ice and the vibratome slicing was performed in chilled 1xPBS (ThermoFisher Scientific).

Fluorescence Correlation Spectroscopy (FCS)
E10.5 transverse spinal cord slices heterozygous or homozygous for Venus::HES5 were stained on ice for 1.5 hours with 50µM Draq5 (ThermoFisher Scientific) diluted in phenol-red free L-15 (ThermoFisher Scientific) media. FCS experiments and snapshot images of whole spinal cord were carried out using a Zeiss LSM880 microscope with a C-Apochromat 40x 1.2 NA water objective on slices placed directly on a glass-bottomed dish (Greiner BioOne) kept at 37°C and 5%CO 2 . FCS signals were collected inside single nuclei in either the ventral region alone or both dorsal and ventral regions for tissue experiments. Venus (EYFP) fluorescence was excited with 514 nm laser light and emission collected between 517 and 570nm. Data from individual cell nuclei was collected using 5 x 2 s runs at 0.15 to 0.3% laser power which gave <10% bleaching and a suitable count rate ~1 kHZ counts per molecule (CPM). To obtain molecule number, autocorrelation curves were fit to a two-component diffusion model with triplet state using the Levenberg-Marquardt algorithm in MATLAB optimization toolbox with initial conditions assuming a 'fast' diffusion component 10x faster than the 'slow' component as previously described 52 . Measurements collected from cells exhibiting large spikes/drops in count rate or with low CPM (<0.5 kHz), high triplet state (>50%), or high bleaching (>10%) were excluded from the final results. Number and brightness analysis of the count rate 53 showed a high correlation with molecule number obtained from autocorrelation curve fitting. The effective confocal volume had been previously determined with mean 0.57fL ± 11 fL (S.D.) allowing conversion from molecule number to concentration 54 . Single-cell data of number of molecules in the cell nucleus was obtained by adjusting concentration to the average volumetric ratio between nuclear volume and confocal volume. Mean nuclear volume of 523 fL was estimated using H2BmCherry intensity and 3D reconstruction from z-stack images in Imaris (Bitplane).

Generating a quantitative expression map
Individual Draq5+ nuclei in a tile-scan image of a transverse slice of the whole E10.5 spinal cord were manually segmented as ellipses using ImageJ and the Venus background intensity subtracted. A quantile-quantile plot was generated for the distribution of nuclear Venus::HES5 intensities from manual segmentation of a single image and the distribution of nuclear Venus::HES5 concentrations from FCS of cells throughout the E10.5 spinal cord from multiple slices and experiments. Linear regression was used to generate a calibration curve between Venus::HES5 intensity and Venus::HES5 concentration over the middle 90% of the range. The gradient of the line was used as a scaling factor and applied to the pixel intensity values in the segmented image to transform intensity to concentration.

Analysis of variability in Venus::HES5 in snapshot images
The centroids of the manually segmented cells from a quantitative expression map were used to measure distance from the ventricle and perpendicular to the D/V axis. Neighbours were ranked based on distance from the centroid of the cell of interest and the nearest neighbours were classified as the cells in the first rank (Supp. Fig.  2f). Coefficient of variation of Venus::HES5 intensity was measured by manual segmentation of Draq-5 stained transverse slices of whole E10.5 spinal cord in ImageJ.

Embryo slice culture and live imaging
E10.5 spinal cord slices for live timelapse microscopy were placed on a 12mm Millicell cell culture insert (MerckMillipore) in a 35mm glass-bottomed dish (Greiner BioOne) incubated at 37°C and 5%CO 2 . The length of the legs of the cell culture insert were reduced to decrease the distance from the glass to the tissue. 1.5mls of DMEM F-12 (ThermoFisher Scientific) media containing 4.5mg/ml glucose, 1x MEM nonessential amino acids (ThermoFisher Scientific), 120ug/ml Bovine Album Fraction V (ThermoFisher Scientific), 55µM 2-mercaptoethanol, 1x GlutaMAX (ThermoFisher Scientific), 0.5x B27 and 0.5x N2 was added. Movies were acquired using Zeiss LSM880 microscope and GaAsP detectors with a Plan-Apochromat 20x 0.8 NA objective with a pinhole of 5AU over a z-stack of roughly 70 um every 15 mins for 18-20 hours. DMSO or 2µM DBZ (Tocris) was added to media immediately before imaging.

Image analysis and cell tracking
Briefly, single cells were tracked using the H2B::mCherry channel. Single-cell Venus and mCherry expression were normalised to the whole tissue mean for the relevant channel to account for any possible photobleaching. For hierarchical clustering single-cell Venus::HES5 expression from 12-hour tracks was standardized by subtracting the mean and dividing by the standard deviation of the single-cell signal.
Single neural progenitor cells in E10.5 spinal cord slices were tracked in Imaris on the H2BmCherry channel using the 'Spots' and 'Track over time' function. Spot detection algorithm used background subtraction and tracking used the Brownian motion algorithm. All tracks were manually curated to ensure accurate single-cell tracking. A reference frame was applied to the movie along the dorso-ventral and apico-basal axes of the spinal cord to allow the distance from the ventricle to be calculated. To account for any photobleaching and allow comparison of intensities between movies the mean intensity of mCherry and Venus in each spot was normalised to the mean intensity of mCherry or Venus in the whole tissue. The whole tissue volume was tracked using the 'Surfaces' and 'Track over time' function.
There was no correlation in Venus::HES5 and H2BmCherry expression suggesting the Venus::HES5 dynamics were not a result of global changes in transcription or translation in the cell or microscope anomalies (Supp . Fig 3a-d and further examples in Supp. Fig. 6). We also investigated the relationship between Venus::HES5 and zposition in the tissue (Supp Fig. 4a-d). As expected from imaging through tissue there was a small negative correlation (r = -0.24) between Venus::HES5 intensity and z-position when all cells and time-points were plotted (Supp. Fig. 4b). However the range of z-positions in a single cell 12-hour track was rarely greater than 25µm, therefore it is unlikely the fluctuations and oscillations in Venus::HES5 are a result in changes in z-position (Supp. Fig. 4c). Further at the single-cell level there is no difference in the correlation coefficient between z-position and Venus::HES5 intensity when comparing oscillatory and non-oscillatory cells (Supp. Fig. 4d).

Hierarchical clustering
Prior to analysis, timeseries of single cell Venus::HES5 expression were normalised to tissue mean to account for bleaching per independent experiment and in addition standardised (z-score calculation) by subtracting the mean of the timeseries from each timepoint and dividing by the standard deviation of the timeseries. Standardising the data enables clustering on relative expression changes rather than absolute expression levels. Standardized single cell timeseries were then subject to hierarchical clustering using Euclidean distance and Ward's linkage in RStudio (R Project). Experiments were clustered separately and each clustergram independently identified 4 clusters per experiment. The elbow method to look at the variance explained as a function of number of clusters (nbclust package, R), suggested 4-6 clusters as the optimal cluster number however 5 and 6 clusters were not favoured by silhouette method (nbclust package, R) so we chose 4 clusters (data not shown). Cluster relationships varied between experiments thus for annotation between experiments corresponding clusters labels were determined by calculating average single-cell coefficient of variation (COV) in Venus::HES5 over time for each cluster and comparing to results of clustering experiment 1 (Supp. Fig. 5b). Thus, four clusters with the same COV profile are reproducibly identified in each experiment.
For DBZ-treated cells, data could not be corrected for photobleaching since Venus::HES5 downregulation is induced at tissue level causing a significant drop in tissue mean and masking effects from bleaching. Prior to analysis both DMSO and DBZ timeseries were standardised by z-scoring. To enable comparison between DBZ-treated and negative control DMSO-treated cells, experimental data from both treatment conditions were clustered together (Fig. 3k) as well as clustering DMSO independently of DBZ (Supp. Fig. 7) yielding similar cluster profiles to untreated cells (Supp. Fig. 5a).

Detection and analysis of oscillations -short term dynamics
We analysed oscillations in the HES5 timeseries using the Gaussian Processes approach in Phillips et al. 42 . Data was de-trended to remove long term behaviour such as down-regulation and to recover the oscillatory signal with zero mean (see Supplementary Methods Section A). We used maximum likelihood estimation to fit the de-trended data timeseries with two competing models: a fluctuating aperiodic one (null model) and an oscillatory one (alternative model). We used the loglikelihood ratio statistic to compare the likelihood of data being oscillatory or nonoscillatory and determined the oscillators based on a false discovery rate of 3% independently per experiment. Additional procedures for dynamic parameter inference were introduced and are detailed in Supplementary Methods Section A (see also Supp. Fig. 8d,f).

Fold change detection in signal amplitude
We used an unsupervised technique of amplitude reconstruction based on the Hilbert transform. This allows us to robustly identify peaks and troughs in noisy timeseries and compute fold changes as peak-to-trough variations in the signal containing the long term trend. This procedure is described in Supplementary Methods.

Stochastic model of genetic auto-repression
We model protein expression dynamics emerging from genetic autorepression and transcriptional delay using an established mathematical model 45 . This model includes the effects of transcription and translation as well as degradation of protein and mRNA. We adjust the model to include noise due to intrinsic stochasticity 43,55 . We parameterized the model using the experimentally measured mRNA and protein halflives (30 and 90 mins respectively, Supp. Fig.1c,d,e,f) and approximate Bayesian computation (ABC) 46,56 (Supp. Fig.10c,d,e, Table 1) to search for parameter combinations which modelled a mean HES5 expression of 55000-65000 protein molecules (Supp. Fig.10a) and standard deviation over time of greater than 5% (Supp. Fig.10b). ABC provides probability distributions of model parameters given experimentally observed data; ABC is a standard method for parameterising stochastic mathematical models. We limit the parameter inference to biophysically realistic parameter ranges (Supp. Table 1&2). Details of the model implementation and parameter inference are available in the supplementary methods (Supplementary Methods B).

Statistical testing
Statistical tests were performed in GraphPad Prism 7. Data was tested for normality with D'Agostino-Pearson test. The relevant parametric or non-parametric test was then performed. If necessary outlier removal was performed using ROUT method (GraphPad). Coefficient of variation is defined as standard deviation (SD) over the mean.
Bar plots (simple or stacked) and discrete scatter plots show mean or mean±SD where multiple independent experiments are analysed. Statistical significance between 2 datasets was tested with either Student t-test (parametric) or Mann-Whitney test (non-parametric). Statistical significance (p<0.05) for 2+ datasets was tested by Kruskall-Wallis with Dunn's multiple comparison correction. All tests were 2-sided. Multiple comparison testing involved comparing all pairs of data columns. Correlations were analysed using Spearman rank correlation coefficient. Sample sizes, experiment numbers, p values<0.05 and correlation coefficients are reported in each figure legend.

Code and Data availability
Data fitting for detections of oscillations has been implemented in Matlab R2015a using the GPML toolbox (Rassmussen and Hannes 2010) and custom designed routines available at http://gaussianprocess.org/gpml/code/matlab/doc/. Code for stochastic model of genetic auto-repression and Bayesian inference available online under https://github.com/kursawe/hesdynamics. Matlab custom designed routines for analysis of FCS available on request. Single-cell Venus::HES5, H2B::mCherry intensity and positional information available on request from the corresponding authors. Mean=61a.u SD=39a.u h) Nuclear Venus::HES5 concentration in homozygous E10.5 Venus::HES5 knock-in embryos across entire spinal cord measured by tissue FCS. 442 cells, 4 experiments. Mean is 148nM, SD=58nM. i) Quantile-quantile plot of nuclear Venus::HES5 concentration (h) vs nuclear Venus::HES5 intensity (g) for E10.5 homozygous embryos. Red line -linear fit over middle 90% range. j) Quantitative map of nuclear Venus::HES5 concentration in whole live E10.5 spinal cord. Intensity values scaled according to linear fit of Q-Q plot in i). Scale bar 50µm. k) Nuclear Venus::HES5 concentration by distance from ventricle in region 1 (upper box in j) and region 2 (lower box in j) l) Concentration difference between a cell and its nearest neighbours for cells less than or greater than 50µm (region 1) from the ventricle (n=154, n=73 respectively. p=0.0007 in Mann-Whitney test), or 30µm (region 2) from the ventricle (n=91, n=135 respectively. p<0.0001 in Mann-Whitney test). m) Coefficient of variation in Venus::HES5 intensity between cells less than or greater than 50µm from the ventricle in ventral domain in E10.5 Venus::HES5 embryos. (n=4 embryos, at least 24 cells per embryo, 2 experiments, p=0.04 in paired t-test.)     Fig.9b). d) Ten example traces generated using the model are shown at three different parameter points in the panels i), ii), and iii). The power spectrum does not have a dominant non-zero peak in i) whereas the power spectra in ii) and iii) do have a dominant nonzero peak with decreasing width from ii) to iii) showing increasing coherence. Parameter values are (i) α m =0.64/min, α p =17.32/min, P 0 =88,288.6, τ=34min, n=5.59 (ii) α m =39.93/min, α p =21.56/min, P 0 =24,201.01, τ=33min, n=4.78 (iii) α m =44.9/min, α p =3.13/min, P 0 =35,080.2, τ=40min, n=5.62. The half-lives of the protein and mRNA are set to 90 and 30 minutes, respectively. e) Response curves in coherence when changing the protein degradation rate. The black line is located at the degradation rate corresponding to a 90 minutes HES5 protein half-life. f-g) Variation in expected coherence for the stochastic model (f) and the deterministic model (g) of HES5 expression as protein and mRNA degradation rates are changed. The blue dots mark experimentally measured values for the protein and mRNA degradation rates, corresponding to a 90 and 30 minute half-life, respectively. Experimentally measured degradation rates are located on the slope of increasing coherence values in the stochastic model, and in a region of no expected oscillations in the deterministic model.