text:
Thank you for visiting nature.com. You are using a browser version with limited support for CSS. To obtain the best experience, we recommend you use a more up to date browser (or turn off compatibility mode in Internet Explorer). In the meantime, to ensure continued support, we are displaying the site without styles and JavaScript.

# Dynamics of lineage commitment revealed by single-cell transcriptomics of differentiating embryonic stem cells

### Subjects

Gene expression heterogeneity in the pluripotent state of mouse embryonic stem cells (mESCs) has been increasingly well-characterized. In contrast, exit from pluripotency and lineage commitment have not been studied systematically at the single-cell level. Here we measure the gene expression dynamics of retinoic acid driven mESC differentiation from pluripotency to lineage commitment, using an unbiased single-cell transcriptomics approach. We find that the exit from pluripotency marks the start of a lineage transition as well as a transient phase of increased susceptibility to lineage specifying signals. Our study reveals several transcriptional signatures of this phase, including a sharp increase of gene expression variability and sequential expression of two classes of transcriptional regulators. In summary, we provide a comprehensive analysis of the exit from pluripotency and lineage commitment at the single cell level, a potential stepping stone to improved lineage manipulation through timing of differentiation cues.

## Introduction

In vitro differentiation is a key technology to enable the use of embryonic and induced pluripotent stem cells as disease models and for therapeutic applications1, 2. Existing directed differentiation protocols, which have been gleaned from in vivo development, are laborious and produce heterogeneous cell populations3. Protocol optimization typically requires costly and time-consuming trial-and-error experiments. To be able to design more efficient and specific differentiation regimens in a systematic way it will be necessary to gain a better understanding of the decision-making process that underlies the generation of cell type diversity4.

Lineage decision-making is fundamentally a single-cell process5 and the response to lineage specifying signals depends on the state of the individual cell. A substantial body of work has revealed lineage biases related to, for example, cell cycle phase or pre-existing subpopulations in the pluripotent state4, 6,7,8. The commitment of pluripotent cells to a particular lineage, on the other hand, has not yet been studied systematically at the single-cell level. We consider a cell to be committed, if its state cannot be reverted by removal of the lineage specifying signal.

Here we set out to characterize the single-cell gene expression dynamics of differentiation, from exit from pluripotency to lineage commitment. Using single-cell transcriptomics we find that retinoic acid drives the differentiation of mouse embryonic stem cells to neuroectoderm—and extraembryonic endoderm—like cells. Between 24 h and 48 h of retinoic acid exposure, cells exit from pluripotency and their gene expression profiles gradually diverge. By pseudotime ordering we reveal a transient post-implantation epiblast-like state. We also study the influence of the external signaling environment and identify a phase of high susceptibility to MAPK/Erk signaling around the exit from pluripotency. We employ a minimal gene regulatory network model to recapitulate the dynamics of the lineage response to signaling inputs. Finally, we identify two classes of transcription factors which have likely distinct roles in the lineage decision-making process.

## Results

### Retinoic acid driven lineage transition

Mouse embryonic stem cells (mESCs) are a well-characterized model system to study in vitro differentiation. Here, we focused on mESC differentiation driven by all-trans retinoic acid (RA), which is widely used in in vitro differentiation assays9 and has important functions in embryonic development10. E14 mESCs were grown feeder free in 2i medium11 plus LIF (2i/L) for several passages to minimize heterogeneity before differentiation in the basal medium (N2B27 medium) and RA (Fig. 1a). Within 96 h the cells underwent a profound change in morphology from tight, round, homogeneous colonies to strongly adherent, morphologically heterogeneous cells (Fig. 1a). To characterize the differentiation process at the population level we first measured gene expression by bulk RNA-seq at 10 time points during 96 h of continuous RA exposure (Supplementary Fig. 1). Genes that are absent in the pluripotent state but upregulated during differentiation can reveal the identity of differentiated cell types. To find such genes we clustered all genes by their temporal gene expression profiles using k-means clustering (Methods, Supplementary Fig. 1a). By testing for reproducibility through repeated clustering (stability analysis12, see Methods) we determined that there were 6 robust gene clusters. The two clusters that showed a continuous increase in expression over the time course (clusters 5 and 6 in Supplementary Fig. 1a), were enriched with genes that have functions in development and differentiation (Supplementary Fig. 1b). In particular, established neuroectoderm and extraembryonic endoderm (XEN) markers belonged to these clusters. Mesodermal markers, on the other hand, were not up-regulated. (Supplementary Fig. 1c, d). This observation is in agreement with earlier reports showing that RA induces neuroectodermal and XEN lineages while suppressing mesodermal gene expression10, 13, 14.

We next set out to identify the final cell types present after 96 h of RA exposure. The up-regulation of both ectodermal and XEN markers seemed to indicate that cells adopted these two fates. Since population level measurements are not able to resolve population heterogeneity, we turned to the recently developed Single Cell RNA Barcoding and Sequencing method15 (SCRB-seq, Supplementary Fig. 2). We quantified the transcriptional profiles of over 2000 single cells, sampled at 9 time points during differentiation, typically spaced 12 h apart. To visualize the heterogeneity of gene expression profiles and find subpopulations that emerge during differentiation, we used principal component analysis (PCA) and k-means clustering of the cells (Fig. 1b, Supplementary Fig. 3a, b). Repeated k-means clustering of the cells (stability analysis12, see Methods) indicated that the population was homogeneous at 0 h and two robust clusters were present at the end of the differentiation time course (96 h). To reveal the identity of the two observed clusters, we turned to the composition of the first two principal components. The first principal component (PC 1) was primarily composed of established markers for the XEN lineage (Sparc, Col4a1, Lama1, Dab2), while PC 2 comprised markers of neuro-ectodermal development (Prtg, Mdk, Fabp5, Cd24) (Supplementary Fig. 3a, b). Accordingly, we identified one cluster as XEN-like and the other one as as ectoderm-like (Fig. 1b). Hierarchical clustering supported our interpretation of the PCA results (Supplementary Fig. 3c). In particular, we observed that genes from gene cluster 5 (Supplementary Fig. 1a), which includes ectoderm markers, were more broadly expressed in the ectoderm-like cells. By contrast, genes from cluster 6, which includes XEN markers, were largely restricted to XEN-like cells.

To confirm the existence of two cell types by an independent method, we next sought to find surface markers that would allow us to identify and purify the cell types. Cd24, which is among the genes with the highest loadings in PC2, is an established marker for neuroectodermal lineages16. Pdgfra is the earliest known marker of the primitive endoderm lineage in vivo17. Antibody staining of these two markers showed two well-separated subpopulations at 96 h (Supplementary Fig. 4a): an ectoderm-like subpopulation (CD24 + /PDGFRA-) and a XEN-like subpopulation (CD24−/PDGFRA + ). The frequencies of these two subpopulations were robust across multiple biological replicates (Supplementary Fig. 4b) and in accordance with the single-cell RNA-seq results. We then purified ectoderm-like and XEN-like cells after 96 h of RA exposure and cultured them in the same medium (N2B27 supplemented with EGF and FGF2). After continued culture, the two subpopulations showed markedly different morphologies (Supplementary Fig. 4c) and distinct gene expression patterns, as measured by bulk RNA-seq (Supplementary Fig. 4d, f). Ectoderm-like cells expressed neuro-ectodermal and neural crest markers and were similar in their expression profile to neural progenitor cells and neural crest cells in vivo. XEN-like cells expressed primitive endoderm markers and resembled an embryo-derived XEN cell line and yolk sac tissue. Taken together, these results provide evidence that the observed cell clusters corresponded to stable neuroectoderm-like and XEN-like cell types with likely in vivo correlates.

### Exit from pluripotency between 24 h and 48 h of RA exposure

Having established the identity of the differentiated cell types we next sought to study the exit from pluripotency in detail. At the population level, we detected a gene expression response to differentiation conditions within only 6 h, as well as a second wave of gene expression changes between 24 h and 36 h (Supplementary Fig. 1e). While the immediate response was a direct effect of the switch to RA containing media, as evident from the upregulation of direct RA targets, we hypothesized that the second wave of changes indicated the exit from pluripotency. In support of this hypothesis we found that pluripotency markers were strongly down-regulated between 24 h and 48 h (Supplementary Fig. 1c, d). Cell morphology and cell cycle phase lengths (Supplementary Fig. 5a–c) also changed significantly during the same time interval, in agreement with the observed expression dynamics. As a functional assay we used replating of the cells at clonal density in 2i/L medium. 90% of the cells could not grow in this selective medium anymore by 36 h of RA exposure (Supplementary Fig. 5d). Taken together, our population level gene expression measurements and functional assays suggested that cells exited pluripotency between 24 h and 48 h of RA exposure.

### Gradual divergence of gene expression profiles

To visualize gene expression dynamics around the exit from pluripotency at the single–cell level we used t-distributed stochastic neighbor embedding18 (t-SNE) of our SCRB-seq data set. t-SNE maps gene expression profiles to a low-dimensional space and places similar expression profiles in proximity to each other. Here we used t-SNE to map the expression profiles of individual cells throughout the time course on a single axis (Fig. 1c). We assessed the pluripotency status of individual cells by the expression level of the established pluripotency marker Rex1 19. t-SNE showed that gene expression changed homogeneously throughout the population for the first 12 h of RA exposure, which was likely a direct effect of the RA containing medium. At this stage Rex1 expression was high throughout the population. The subsequent steep increase in single-cell variability of gene expression at 24 h (Fig. 1d) indicated that gene expression profiles started to become more heterogeneous during the exit from pluripotency. Simultaneously, Rex1 expression started to decline in a subset of cells, confirming the exit from pluripotency at the single-cell level. To pinpoint the time when distinguishable cell types first appeared during the differentiation time course, we calculated gene expression variability for individual cell clusters formed by k-means clustering (Fig. 1d), instead of the whole population. Starting at 48 h, within-cluster variability using 2 clusters was reduced compared to population variability, signifying the emergence of the two cell types. Clustering into 3 or 4 clusters did not reduce the variability much further. Taken together, t-SNE mapping and variability analysis showed that cells exited pluripotency and started to diverge in gene expression between 24 h and 48 h of RA exposure.

To further quantify the divergence of gene expression profiles we classified cells based on their similarity (Pearson correlation) with the average profiles of either mESCs at 0 h or the two differentiated cell types at 96 h (Fig. 2a–c). Cells which were more similar to a differentiated cell type than to mESCs first appeared between 24 h and 48 h of RA exposure, which matched the dynamics visible in the t-SNE map (Fig. 1c). Importantly, average expression profiles of the three classes were similar around the exit from pluripotency and only diverged more quickly afterwards (Fig. 2d). These observations suggested that the cells adopted the final cell fates only gradually, potentially via distinct transitory states.

### Initial differentiation into post-implantation epiblast

We next wanted to zoom in further on the initial lineage decision, right after the exit from pluripotency, to reveal potential intermediate cell states. To achieve this goal, we had to remove possible obfuscating effects related to the asynchrony of differentiation. The transient coexistence of all three classes of cells (Fig. 2c) and the heterogeneous expression of Rex1 (Fig. 1c) around the exit from pluripotency had indicated that differentiation was indeed asynchronous. Confounding effects due to asynchronous differentiation can be mitigated with the help of pseudo-temporal ordering of cells20. Here we defined a pseudo-time based on the Pearson correlation with mESCs or the differentiated cell types at 96 h (Fig. 3a, b). This pseudo-time thus reflects the progress of differentiation of an individual cell along the ectoderm- or XEN-like lineage. Pseudo-temporal ordering reduced the co-existence of cell types to a small period in pseudo-time (Fig. 3c) and was thereby able to clarify expression dynamics. Furthermore, it revealed that pluripotency factors were down-regulated already before the branch point, where differentiated cell types could first be distinguished (Fig. 3d, e). During the same time, markers of post-implantation epiblast21 (e.g. Pou3f1, Fgf5) were up-regulated. This intermediate period might represent a phase of homogeneous lineage priming or subtle population heterogeneity that we cannot resolve given the technical noise of our single-cell RNA-seq method. After the branch point, several neuroectodermal markers (like Pax6, Sox11 or Nes) were up-regulated in the ectoderm-like branch. Established XEN markers (e.g. Gata6, Dab2), on the other hand, were restricted to the XEN-like branch, as to be expected.

Gene expression dynamics in pseudo-time seemed to suggest a transient state in which the cells resembled the post-implantation epiblast. To further clarify the relationship of RA differentiation with in vivo development we used PCA to compare our data set to RNA-seq measurements of pre- and peri-implantation tissues21 (Fig. 4a and Supplementary Fig. 6). This analysis revealed that mESCs were most similar to pre-implantation epiblast (E4.5), as has been shown previously22. During differentiation the cells first moved closer to the E5.5 epiblast around 48 h before separating into two subpopulations (Fig. 4a). At 96 h, the XEN-like subpopulation was closest to E4.5 primitive endoderm. The occurrence of these XEN-like cells is thus likely due to a trans-differentiation from E4.5 or even E5.5 epiblast–like cells. The initial lineage decision in our system is therefore between continued differentiation along the epiblast lineage and trans-differentiation to a primitive endoderm-like state.

To confirm the single-cell RNA-seq results with an independent method we sorted cells based on PDGFRA and CD24 expression at 48 h, 72 h and 96 h and profiled the expression of the sorted subpopulations by bulk RNA-seq (Fig. 4b). At 48 h only few cells expressed PDGFRA but the majority expressed CD24. Most importantly, in PDGFRA negative cells the expression of post-implantation epiblast markers increased with CD24 expression. By 96 h the expression of post-implantation epiblast markers had largely disappeared. XEN markers, on the other hand, were expressed exclusively in PDGFRA positive cells at 72 h and 96 h. To determine cell identities in the bulk expression data set in an unbiased way we used the KeyGenes algorithm23 together with pre- and peri-implantation tissues21 as training set (Fig. 4c). KeyGenes identified mESCs as E4.5 epiblast, in agreement with our PCA (Fig. 4a) and previous results22. Notably, at 48 h PDGFRA negative/CD24 low cells were classified as E4.5 epiblast, while PDGFRA negative/CD24 high cells were identified as E5.5 epiblast. CD24 thus indicated the adoption of a post-implantation epiblast-like state, in agreement with previous findings24. PDGFRA positive cells, on the other hand, were consistently identified as E4.5 primitive endoderm. Bulk RNA-seq of sorted subpopulations and KeyGenes analysis thus confirmed that cells either continued to differentiate along the epiblast lineage or adopted a XEN-like cell type.

### Regulation by the external signaling environment

Having characterized the gene expression dynamics of the exit from pluripotency and the subsequent lineage transition, we next wanted to identify effectors of the lineage decision. Notably, mESCs lost their ability to differentiate into a XEN-like lineage when they were cultured, prior to differentiation, in serum and LIF conditions (without feeders) instead of 2i/L (Supplementary Fig. 7a, b). The ability of RA to drive ectodermal differentiation seemed unaffected under these conditions, as reported before25. Since culture conditions had such a strong impact on the developmental potential of mESCs we wanted to explore the contribution of specific signaling pathways on the cellular decision. We differentiated mESCs with RA in the presence of a MEK inhibitor (MEKi, PD0325901), which abrogates MAPK/Erk signaling; a GSK3 inhibitor, which effectively stimulates Wnt signaling (GSK3i, CHIR99021), LIF, which activates the JAK/Stat pathway or an FGF receptor inhibitor (FGFRi, PD173074). (Supplementary Fig. 7c–h). The first 2 of these molecules are components of the defined 2i medium and are known to prevent differentiation while stabilizing the pluripotent state. The presence of GSK3i or LIF led to an overall reduction of differentiated cells (Supplementary Fig. 7c), consistent with their role in stabilizing pluripotency. Addition of MEKi alone, however, led to a specific reduction of the XEN-like subpopulation (Supplementary Fig. 7c–e), in agreement with previous results26, 27. This effect was unlikely due to interference with RA signaling since increasing RA concentration did not reverse the effect (Supplementary Fig. 7f). In contrast to the MEK inhibitor, the FGF receptor inhibitor not only suppressed the XEN-like population but also greatly reduced the ectoderm-like population (Supplementary Fig. 7g, h). This observation is in agreement with earlier studies that reported a requirement for FGF signaling in mESC differentiation28 and lineage segregation in the early mouse blastocyst29. Taken together these experiments clearly demonstrate that RA driven XEN-specification requires the same signaling pathways as other differentiation regimens and XEN-specification in vivo, despite the pleiotropic nature of RA.

### Phase of high susceptibility to external signal inputs

We next wanted to establish when mESCs are sensitive to RA signaling and how long the signal would have to be applied to drive a complete lineage transition. Having observed that gene expression responds to differentiation conditions within 6 h (Supplementary Fig. 1e), we hypothesized that a short pulse of RA might be sufficient to induce XEN specification. To test this hypothesis we applied a precisely defined pulse of RA by first exposing the cells to RA for a defined period of time and then switching to a highly potent pan-RA receptor antagonist30 (Fig. 5a). These experiments showed that, contrary to our expectation, RA had to be applied for at least 24 h for XEN-like cells to appear. Longer pulses resulted in a gradual increase of the XEN-like fraction. A 36 h long pulse of RA resulted in 20% XEN-like cells at the 96 h time point, roughly half of what we found after uninterrupted RA exposure (Fig. 5a). This indicated that even after 36 h of RA exposure and significant down-regulation of the pluripotency network XEN specification continued to depend on RA-signaling. Timed abrogation of MAPK/Erk signaling by MEKi resulted in a similar response as an RA pulse (Fig. 5b). At least 24 h of uninterrupted MAPK/Erk signaling was necessary for XEN-like cells to occur. Longer durations of MAPK/Erk signaling resulted in an increase in the XEN-like subpopulation. This effect plateaued after 48 h, which suggested that XEN-like cells then became independent of MAPK/Erk signaling and thus stably committed. We also wanted to establish when cells lost their ability to respond to RA signaling. To this end we first differentiated the cells in basal (N2B27) medium and started RA exposure after a defined time period (Fig. 5c). When RA exposure was delayed by up to 12 h, we did not observe any difference in the lineage distribution at the 96 h time point. For longer delays of RA exposure, we found that the fraction of XEN-like cells declined. This observation demonstrated that the cells quickly lost their susceptibility to RA after the exit from pluripotency. Taken together, these signaling experiments revealed a short transient phase after the exit from pluripotency, during which cells were maximally susceptible to external signaling cues to inform their lineage decision.

### A minimal gene regulatory network of lineage bias

Interestingly, our experiments revealed a difference in the lineage response dynamics between the RA pulse and RA delay. While cells abruptly lost their ability to become XEN-like after only 12 h in N2B27 (Fig. 5c), the RA pulse had to be applied for at least 24 h to cause XEN-specification and longer pulses elicited a gradually increasing response (Fig. 5a). This asymmetry could be related to the fact that N2B27 on its own drives differentiation towards the neuroectoderm lineage31. Correspondingly, we consistently found that the majority of cells became ectoderm-like when there was no RA present (Fig. 5a–c). To explore the role of an intrinsic epiblast or ectoderm bias we developed a simple phenomenological model based on a minimal gene regulatory network (GRN)27, 32. Briefly, the GRN is comprised of two lineage-specific, auto-activating expression programs that mutually repress each other (Fig. 5d, Supplementary Fig. 8a). This GRN can produce two stable attractors that correspond to two differentiated cell types. Here, we added repression of both lineages by the pluripotency network to model the pluripotent state. Consistent with our data, we assumed that the pluripotency program is turned off after 12 h. To model the ectoderm bias we assumed that auto-activation of the ectoderm program was stronger than auto-activation of the XEN program in the absence of RA. In the presence of RA auto-activation of both programs was taken to be equal. Due to the great importance of gene expression noise in lineage decision-making5, 4, we also incorporated noise in our model (Methods). Stochastic simulations of the 3-state GRN reproduced the asymmetry between RA pulses and delays (Fig. 5e). The frequency of XEN-like cells decreased sharply when the delay in RA signaling was extended beyond the exit from pluripotency. The RA pulse, on the other hand, had to be applied for a longer period of time to cause XEN-specification and the response was more gradual. This behavior can be explained by the fact that in the absence of RA cells are quickly drawn to the ectoderm attractor after the exit from pluripotency. When RA is added after a delay, the cells are already in the proximity of the ectoderm attractor and cannot escape it anymore, which causes the lack of XEN cells. Notably, the asymmetry between the response curves was reduced by gene expression noise. Noise allowed the cells to switch between the basins of attraction of the two attractors (Supplementary Fig. 8b), thereby equalizing the intrinsic difference between the two attractors. Taken together, our stochastic simulations showed that an intrinsic ectoderm bias can explain the difference in the response dynamics between an RA signal delay and an RA pulse.

### Two classes of transcriptional regulators

Having revealed a highly dynamic susceptibility to signaling cues, we were wondering if the expression of transcriptional regulators was equally dynamic. To that end we focused on transcriptional regulators that show lineage specific expression when the two lineages can be first discerned robustly, around 48 h (see Methods for the list of GO terms used to define transcriptional regulators). Since these regulators were typically lowly expressed, they were not well-represented in the SCRB-seq data set. Therefore, we collected another single-cell RNA-seq data set using SMART-seq233 at four early RA differentiation time points (0 h, 12 h, 24 h and 48 h). We first identified XEN-like and ectoderm-like cells at the 48 h time point (Supplementary Fig. 9a). The remaining cells were likely mostly undifferentiated cells as several pluripotency factors were differentially expressed in this population (Supplementary Fig. 9b). In the cells classified as XEN- or ectoderm-like we found 50 transcriptional regulators to be differentially expressed between the two lineages (Fig. 6a, Supplementary Fig. 9b). 22 of those genes (dubbed “early”) were present already in mESCs. These early regulators were broadly co-expressed in individual cells at the beginning of the time course (Fig. 6b and Supplementary Fig. 9c). Compared to canonical pluripotency factors, early regulators showed a smaller level of co-expression with each other in the pluripotent state, in particular if they belonged to different lineages (Fig. 6b, c). Individual mESCs thus expressed varying ratios of XEN and ectoderm specific early regulators. Over time, co-expression of XEN and ectoderm specific early regulators declined but they never became completely mutually exclusive (Supplementary Fig. 9c). Hence, we speculated that other transcriptional regulators might be up-regulated in lineage biased cells and take over lineage specification from the early regulators. Indeed, 28 of the identified differentially expressed regulators (dubbed “late”) were, by definition, not significantly expressed at the beginning of the time course (Fig. 6a, Supplementary Fig. 9b). These late regulators were overall positively correlated with early regulators of the same lineage and anti-correlated with regulators of the opposing lineage (Fig. 6d and Supplementary Fig. 9c). This correlation pattern suggested that early regulators might have a role in lineage biasing, whereas late factors could be involved in lineage commitment.

To confirm the sequential expression of early and late regulators, we next focused on four transcription factors, chosen based on their reported function for the specification of ectoderm (Gbx2 34 (early), Pax6 35 (late)) and extraembryonic endoderm (Tbx3 36 (early), Gata6 37 (late)). Notably, Tbx3 and likely also Gbx2 are direct targets of RA38, 39. In agreement with their reported roles we found these 4 factors to be differentially expressed in ectoderm-like and XEN-like cells, respectively, in our SCRB-seq data set (Supplementary Fig. 9d). To quantify correlation patterns with high precision we used single-molecule FISH (smFISH40) due to its superior sensitivity and precision compared to single-cell RNA-seq (Supplementary Fig. 10a). We measured the expression of the early factors (Fig. 7a, c) or the late factors (Fig. 7b, d) together with the pluripotency factor Nanog and quantified co-expression at all time points (Supplementary Fig. 10b–d). In agreement with the SMART-seq2 data, early factors were broadly co-expressed in the pluripotent state and a smaller subpopulation of co-expressing cells persisted during differentiation (Fig. 7c, Supplementary Fig. 10c). Importantly, mESCs expressed the early factors at highly variable ratios: 30% of mESCs did not express the early ectoderm factor Gbx2 at a significant level, while almost all cells expressed the early XEN factor Tbx3 (Fig. 7e). smFISH further confirmed that late factors were only sporadically expressed before the exit from pluripotency but strongly up-regulated in separate subpopulations thereafter. These subpopulations likely corresponded to lineage-committed cell states (Fig. 7d and Supplementary Fig. 10d, e). Interestingly, a simultaneous measurement of the early ectoderm factor Gbx2 and the late ectoderm factor Pax6 revealed their positive correlation throughout the time course, even before the exit from pluripotency (Fig. 7f). A possible explanation for such a correlation might be a lineage-biasing role for Gbx2. All in all, the smFISH measurements clearly confirmed differences in the expression dynamics and correlation patterns of early and late transcriptional regulators.

## Discussion

In summary, we leveraged a recently developed high-throughput single-cell transcriptomics method to dissect the exit from pluripotency and dynamics of lineage commitment in RA driven differentiation of mESCs with high temporal resolution. We characterized the influence of the external signaling environment and explained the dynamics of the signaling response with a minimal gene regulatory network. We finally identified potential transcriptional regulators of lineage decision and commitment.

In particular, we showed that after 96 h of RA exposure mESCs had differentiated into neuroectoderm-like and XEN-like cells. By purification and continued culture we showed that these cell types are stable and not just transient expression fluctuations. In agreement with previous results22 we found mESCs cultured in 2i/L to be transcriptionally most similar to E4.5 epiblast in vivo (Fig. 4a–c). At E4.5 the lineage decision between primitive endoderm and epiblast has already occurred, so a priori it would not be expected that mESCs should be able to generate XEN cells. The potential to create XEN-like cells could be explained by a subpopulation of cells in the pluripotent state that resembles an earlier developmental stage. In our single-cell RNA-seq data set we could not find evidence for such pre-existing heterogeneity (Fig. 4a). Alternatively, RA might have caused the dedifferentiation of the whole mESC population to an earlier developmental stage after which the cells could follow the in vivo bifurcation between E4.5 epiblast and primitive endoderm. While the whole population indeed initially moved closer to the E3.5 inner cell mass during the first 24 h, cells then moved towards E5.5 epiblast before discernible XEN-like cells appeared (Supplementary Fig. 6). Hence, most likely XEN-like cells are created by trans-differentiation from E4.5 or E5.5 epiblast-like cells and mESCs initially decide between progression along the epiblast lineage and the XEN-like cell type right after the exit from pluripotency. The epiblast lineage then further develops to neuroectoderm-like cells by 96 h. A recently published study by Klein et al. used single-cell RNA-seq to characterize mESC differentiation by LIF withdrawal41 and also found a small XEN-like subpopulation. That and other studies42, 27 show that XEN-like cells occur more generally in in vitro differentiation of mESCs and are not an idiosyncratic artefact of exposure to RA. We also found that mESCs grown in 2i/L (but not in serum and LIF) efficiently generate XEN cells under RA exposure (Supplementary Fig. 7a, b). Similarly, Schröter et al. have observed, for a different differentiation assay, that pre-culture in 2i/L greatly increases the number of XEN-like cells generated from mESCs27. Together with those results our observations thus support a model in which mESCs grown in 2i/L functionally correspond to a slightly earlier developmental stage than mESCs grown in serum and LIF42.

Despite the artificial nature of the lineage transition described here, we observed several similarities with the epiblast/primitive endoderm bifurcation in vivo. A recent study by Saiz et al. in the mouse embryo43 showed that epiblast and primitive endoderm are specified asynchronously from a pool of progenitor cells, which also happened in our experiments (Fig. 2c). Furthermore, in the absence of primitive endoderm inducing signals, all cells of the inner cell mass become epiblast-like in vivo29, 43,44,45. In our experiments the majority of mESCs adopted the epiblast/ectodermal lineage in the absence of RA, in agreement with the literature31, 46. Saiz et al. observed that MEKi prevented the specification of primitive endoderm, in agreement with an earlier report by Nichols et al47. The experiments by Saiz et al. also revealed that the susceptibility to MEKi disappeared gradually between E3.5 and E4.5. Our experiments with MEKi showed similar dynamics (Fig. 5b). Thus, both in vivo and in vitro, cells seem to gain competence to respond to XEN specifying signals over time. In vitro, the susceptibility to signaling inputs is thought to be contingent on the down-regulation of pluripotency factors and the exit from pluripotency46, 48, a notion which is supported by our study (Fig. 5a–c). Our results thus clearly reveal a window of opportunity right after the exit from pluripotency, which might be exploited to guide lineage decisions with maximal efficacy.

Using a minimal GRN to model the lineage decision, we also showed that an inherent epiblast/ectoderm bias can cause the observed asymmetry between an RA delay and an RA pulse (Fig. 5d). A similar GRN has been used successfully before in a report by Schröter et al., studying induction of the XEN lineage by exogenous Gata4 expression27. Importantly, our model does not strictly require an ectoderm bias. An initial bias for progression along the epiblast lineage (and continued differentiation to ectoderm under RA) would be sufficient. This interpretation is in line with previous results that reported expression of non-ectodermal markers during early stages of differentiation in N2B2749. Notably, in our model, gene expression noise was able to reduce the asymmetry between the RA regimens, because gene expression trajectories could switch more easily between the basins of attraction of the two lineage attractors (Supplementary Fig. 8b). The impact of noise in the context of lineage decisions was recently addressed in a publication by Marco et al50. In that study the authors focused on the ability of noise to destabilize committed cell states. Here we showed that noise can also impact commitment dynamics and even mask an intrinsic lineage bias. This result suggests that gene expression noise could be exploited to influence lineage decision-making in vitro.

Our study further identified early-expressed lineage specific transcriptional regulators that are heterogeneously expressed in the pluripotent state and thus have a potential role in biasing the lineage decision. Importantly, the two factors we studied in detail, Gbx2 and Tbx3, were previously determined to be part of an essential pluripotency network51,52,53,54,55. It has been suggested before that some pluripotency genes are also involved in lineage specification25, 48, 49. Thomson et al. showed that Sox2 and Oct4 promote the neuroectodermal and mesendodermal lineage, respectively48. Malleshaiah et al. reported similar functions for Nac1 and Tcf3, respectively25. Future research will have to show whether Gbx2 and Tbx3 have similar roles for the epiblast/neuroectoderm and XEN lineage, respectively. In fact, for Tbx3 Lu et al. recently demonstrated a dual function in self-renewal and XEN specification36. The observed correlation between Gbx2 and Pax6 suggests a function of Gbx2 in epiblast or neuro-ectoderm specification. The long-tail distribution of Gbx2 in mESCs hints at infrequent transcriptional bursting and possibly distinct subpopulations56. The causal relationship between Gbx2 and Pax6 and the functional relevance of the Gbx2 high subpopulations will be explored in a future study. Late-expressed lineage specific transcription factors, like Pax6 and Gata6, which were not expressed in the pluripotent state, have a role in lineage commitment. They can thus serve as bona fide lineage markers.

Transient phases of susceptibility to lineage cues, such as the one characterized in this study, might be valuable windows of opportunity for the control of lineage decisions. We speculate that exit from a pluripotent cell state necessarily coincides with a phase of instability and increased gene expression variability, as demonstrated recently for lineage decisions in the hematopoietic system57, 58. Based on our results we would like to propose tentative transcriptional signatures of such phases (Fig. 8): 1. down-regulation of pluripotency factors (Fig. 1c), 2. a sudden increase in single-cell gene expression variability (Fig. 1d), 3. slowly diverging lineage specific expression patterns (Fig. 2d), 4. co-expression of early-expressed (thus potentially lineage-biasing) transcriptional regulators (Fig. 6b), 5. sporadic expression of late-expressed (thus potentially lineage-committing) transcriptional regulators (Fig. 7d). We hope that these results will be a stepping stone towards finding more efficient ways to guide lineage decisions.

## Methods

### Cell culture

All cell lines were grown routinely in modified 2i medium11 plus LIF (2i/L): DMEM/F12 (Life technologies) supplemented with 0.5x N2 supplement, 0.5x B27 supplement, 0.5mM L-glutamine (Gibco), 20 µg/ml human insulin (Sigma-Aldrich), 1 × 100U/ml penicillin/streptomycin (Gibco), 0.5x MEM Non-Essential Amino Acids (Gibco), 0.1 mM 2-Mercaptoethanol (Sigma-Aldrich), 1 µM MEK inhibitor (PD0325901, Stemgent), 3 µM GSK3 inhibitor (CHIR99021, Stemgent), 1000 U/ml mouse LIF (ESGRO). Cells were passaged every other day with Accutase (Life technologies) and replated on gelatin coated tissue culture plates (Cellstar, Greiner bio-one).

E14 cells were provided by A. van O., V6.5 cells were provided by R.J. Both cell lines were regularly tested for mycoplasma infection.

### Differentiation

Prior to differentiation cells were grown in 2i/L for at least 2 passages, with the exception of the experiment shown in Supplementary Fig. 7a,b: here cells were grown in knockout DMEM (Thermofisher) supplemented with 10% ES cell screened FBS (Sigma), 1 × 100U/ml penicillin/streptomycin (Gibco), 0.1 mM 2-Mercaptoethanol (Sigma-Aldrich) and 1000 U/ml mouse LIF (ESGRO) for 3 passages prior to differentiation. For all differentiation experiments cells were seeded at a density of 2.5 × 105 cells per 10 cm dish and grown over night (12 h). After washing cells twice with PBS, differentiation was carried out in basal N2B27 medium (2i/L medium without the inhibitors, LIF and the additional insulin) supplemented with all-trans retinoic acid (RA, Sigma-Aldrich). RA concentration was 0.25 µM unless stated otherwise. Spent medium was exchanged with fresh medium after 48 h.

For the RA pulse experiments (Fig. 5a) cells were first differentiated with 0.25 µM RA for the indicated amounts of time, washed three times with PBS and cultured in basal medium with 2.5 µM of the RA receptor antagonist AGN 193109 (sc-210768, Santa Cruz Biotechnology). At this concentration this antagonist completely inhibits signaling through all-trans retinoic acid30.

For the differentiation under perturbation of various signaling pathways (Supplementary Fig. 7c) we used the MEK inhibitor PD0325901 (Stemgent, standard concentration 1 µM or dilutions thereof), GSK3 inhibitor CHIR99021 (Stemgent, standard concentration 3 µM or dilutions thereof) or mouse LIF (ESGRO, 1000 U/ml). For the experiments with MEK inhibition shown in Fig. 5b and Supplementary Fig. 7d,e we used PD0325901 at a concentration of 0.5 µM. For differentiation under inhibition of FGF signaling, shown in Supplementary Fig. 7g,h we used the FGF receptor inhibitor PD173074 (Sigma-Aldrich) at a concentration of 1 µM.

Multiple biological replicates of the differentiation of E14 cells with RA were performed, where replicates were characterized with different methods to cross-validate the results: SCRB-seq (1 replicate), SMART-seq2 (1 replicate), smFISH (3 replicates where 2 replicates used the same probe set), antibody staining (3 replicates). Morphologies similar to the ones shown in the representative images in Fig. 1a and Supplementary Fig. 5b were observed in at least 5 independent biological replicates of the experiment.

### Long-term culture of differentiated cells

Cells that were differentiated for 96 h with RA were sorted into ectoderm-like (CD24A + /PDGFRA-) and XEN-like (PDGFRA + /CD24-) and replated on poly-D-lysine and laminin coated tissue culture dishes in basal (N2B27) medium complemented with 20 ng/ml mouse EGF (E5160, Sigma) and 10 ng/ml mouse FGF2 (SRP4038-50UG, Sigma). Ectoderm-like cells were propagated by dissociation with Accutase (Life Technologies) and replating under identical conditions every 3–4 days. Floating aggregates of XEN-like cells were propagated in suspension in uncoated plastic petri dishes. Aggregates were not dissociated but the medium was refreshed typically every 4 days. Morphologies similar to the ones shown in the representative images in Supplementary Fig. 4c were observed in 3 independent biological replicates of the experiment.

### Antibody staining and FACS sorting

We used the following antibodies: APC Rat Anti-Mouse CD24 (BD Bioscience, 562349), PE Rat Anti-Mouse CD24 (BD Bioscience, 553262), Anti-Mouse CD140a (PDGFRA) FITC (eBioscience,17-1407), Anti-Mouse CD140a (PDGFRA) APC (eBioscience,17-1401), all at a dilution of 1:1000. Cells growing in 6-well plates were washed once with PBS and then incubated in a volume of 500 µl of basal (N2B27) medium with antibodies for 30 min at 37 °C, in the dark. Subsequently, cells were washed once with PBS, 300 µl Accutase (Life Technologies) was added and cells were gently dissociated by pipetting up and down. After adding 600 µl of basal medium the cell suspension was loaded on a flow cytometer (LSR II, BD Bioscience) or cell sorter (FACSAria III, BD Bioscience). Cells growing in 10 cm dishes were first dissociated and incubated in 1 ml medium with the same incubation conditions and antibody concentrations as for adherent cells. After staining in solution, cells were spun down, the supernatant was removed and cells were resuspended in 1 ml of basal medium before flow cytometry or sorting.

Sorting gates for positive and negative populations were set by comparison to the signal measured in undifferentiated mESCs. For the experiments shown in Fig. 4b, c cells were sorted according to quartiles of CD24 signal at 48 h or terciles of CD24 signal at 72 h.

### Colony formation assay

Cells were differentiated with or without RA as described above for various amounts of time and then replated at a density of 5 × 104 cells/well in a gelatinized 6-well tissue culture plate in 2i/L. Colonies were grown for 2 additional days, washed twice with PBS and then imaged in PBS. Remaining colonies were counted automatically by a custom made image analysis script written in MATLAB. The number of surviving colonies was normalized to the first data point (replating of untreated cells growing in 2i/L).

### Measurement of cell cycle phases

Cells growing on gelatinized tissue culture dishes were washed twice with PBS, detached with Accutase (Life technologies) and resuspended in full medium. Formaldehyde was added to the cell suspension to a final concentration of 4%. Cells were incubated for 12 min at room temperature while being rotated and then spun down for 3 min at 90 x g. Subsequently cells were permeabilized at least over night in 70% ethanol. Cells were stained with Hoechst 33342 in PBS for 1 h and fluorescence measured on a flow cytometer (LSR II, BD Biosciences). The Dean-Jet-Fox model59 was fit to histograms of the fluorescence signal to determine the relative lengths of the cell cycle phases reported in Supplementary Fig. 5c.

### Single cell isolation for SCRB-seq

For each differentiation time point cells were harvested and medium removed by spinning for 5 min at 90 x g. RNA was stabilized by immediately resuspending the pelleted cells in RNAprotect Cell Reagent (Qiagen) and RNaseOUT Recombinant Ribonuclease Inhibitor (Life Technologies) at a 1:1000 dilution. Just prior to fluorescence-actived cell sorting (FACS), the cells were diluted in PBS and stained for viability using Hoechst 33342 (Life Technologies). 384-well SBS capture plates were filled with 5 μl of a 1:500 dilution of Phusion HF buffer (New England Biolabs) in water and individual cells were then sorted into each well using a FACSAria II flow cytometer (BD Biosciences) based on Hoechst DNA staining. After sorting, the plates were immediately sealed, spun down, cooled on dry ice and then stored at −80°C.

### SCRB-Seq of isolated single cells

Frozen cells were thawed for 5 min at room temperature and cell lysis was enhanced by a treatment with proteinase K (200 μg/mL;Ambion) followed by RNA desiccation to inactivate the proteinase K and simultaneously reduce the reaction volume (50 °C for 15 min in sealed plate, then 95 °C for 10 min with seal removed).

To start, diluted ERCC RNA Spike-In Mix (1 μl of 1:107; Life Technologies) was added to each well and the template switching reverse transcription reaction was carried out using Maxima H Minus Reverse Transcriptase (Thermo Scientific), our universal adapter E5V6NEXT (1 pmol, Eurogentec):

5′-iCiGiCACACTCTTTCCCTACACGACGCrGrGrG-3′

where iC: iso-dC, iG: iso-dG, rG: RNA G, and our barcoded adapter E3V6NEXT (1 pmol, Integrated DNA Technologies):

5′-/5Biosg/ACACTCTTTCCCTACACGACGCTCTTCCGATCT[BC6]N10T30VN-3′

where 5Biosg = 5′ biotin, [BC6] = 6 bp barcode specific to each cell/well, N10 = Unique Molecular Identifiers. Following the template switching reaction, cDNA from 384 wells was pooled together, and then purified and concentrated using a single DNA Clean & Concentrator-5 column (Zymo Research). Pooled cDNAs were treated with Exonuclease I (New England Biolabs) and then amplified by single primer PCR using the Advantage 2 Polymerase Mix (Clontech) and our SINGV6 primer (10 pmol, Integrated DNA Technologies):

5′-/5Biosg/ACACTCTTTCCCTACACGACGC-3′

Full length cDNAs were purified with Agencourt AMPure XP magnetic beads (0.6x, Beckman Coulter) and quantified on the Qubit 2.0 Flurometer using the dsDNA HS Assay (Life Technologies). Full-length cDNA was then used as input to the Nextera XT library preparation kit (Illumina) according to the manufacturer’s protocol, with the exception that the i5 primer was replaced by our P5NEXTPT5 primer (5 μM, Integrated DNA Technologies):

5′-AATGATACGGCGACCACCGAGATCTACACTCTTTCCCTACACGACGCTCTTCCG *A*T*C*T*-3′

where * = phosphorothioate bonds.

The resulting sequencing library was purified with Agencourt AMPure XP magnetic beads (0.6x, Beckman Coulter), size selected (300–800 bp) on a E-Gel EX Gel, 2% (Life Technologies), purified using the QIAquick Gel Extraction Kit (Qiagen) and quantified on the Qubit 2.0 Flurometer using the dsDNA HS Assay (Life Technologies). Libraries were sequenced on Illumina Hiseq paired-end flow cells with 17 cycles on the first read to decode the well barcode and UMI, a 9 cycle index read to decode the i7 Nextera barcode and finally a 46 cycle second read to sequence the cDNA.

### RNA-seq on bulk samples

Bulk RNA-seq samples comprise complete populations at 10 time points during RA differentiation (Supplementary Fig. 1) as well as various sorted subpopulations (Fig. 4b, c) and long term cultured ectoderm- and XEN-like cells (Supplementary Fig. 4d, e). Cells were collected in RNAprotect, lysed in QIAzol (Qiagen) and total RNA was extracted and purified using Direct-zol RNA MiniPrep (Zymo Research). DGE libraries were prepared from 10 ng of extracted total RNA, using the protocol described above for SCRB-seq with the exception of using more concentrated E3V6NEXT and E5V6NEXT (10 pmol).

### SCRB-seq and bulk RNA-seq read alignment

All second sequence reads were aligned to a reference database consisting of all mouse RefSeq mRNA sequences (obtained from the UCSC Genome Browser mm10 reference set: http://genome.ucsc.edu/), the mouse mm10 mitochondrial reference sequence and the ERCC RNA spike-in reference sequences using bwa version 0.7.4 with non-default parameter “-l 24”. Read pairs for which the second read aligned to a mouse RefSeq gene were kept for further analysis if 1) the initial six bases of the first read all had quality scores of at least 10 and corresponded exactly to a designed well-barcode and 2) the next ten bases of the first read (the UMI) all had quality scores of at least 30. Digital gene expression (DGE) profiles were then generated by counting, for each microplate well and RefSeq gene, the number of unique UMIs associated with that gene in that well. Python scripts implementing the alignment and DGE derivation are available from the authors upon request.

### SMART-seq sample preparation and read alignment

The single-cell SMART-seq2 libraries were prepared according to the SMART-seq2 protocol33, 60 with some modifications61. Briefly, total RNA from single cells sorted in lysis buffer was purified using RNA-SPRI beads. Poly(A) + mRNA from each single cell was converted to cDNA which was then amplified. cDNA was subjected to transposon-based fragmentation that used dual-indexing to barcode each fragment of each converted transcript with a combination of barcodes specific to each single cell. Barcoded cDNA fragments were then pooled prior to sequencing. Sequencing was carried out as paired-end 2 × 25 bp with 8 additional cycles for each index. Alignment of the reads and calculation of gene expression was done through the Tuxedo pipeline (Tophat, Cuffquand, Cuffnorm)62. Gene expression was expressed as reads per kilobase exon model per million mapped reads (RPKM).

### Computational analysis bulk RNA-seq experiments

The bulk RNA-seq results were normalized by the total amount of reads per time point. Only those genes with non-zero mean were considered for further analysis. For k-means clustering of the temporal profiles we first determined the number of robust clusters. Stability analysis12 indicated that there were 6 robust clusters (Supplementary Fig. 1a). We then performed gene ontology enrichment analysis using the DAVID bioinformatics resource63 the results of which are summarized in Supplementary Fig. 1b. Only the clusters of monotonically upregulated genes (clusters 5 and 6) showed significant enrichment for GO terms related to development, morphogenesis and differentiation. The heat maps of bulk RNA-seq data depict expression relative to Gapdh expression (Supplementary Fig. 1c). To quantify global changes in gene expression we calculated the L2 norm (Euclidean norm) for individual time points including all genes with non-zero average expression across all time points. Differences in the L2 norm between time points are reported in Supplementary Fig. 1e.

To reveal the identity of sorted subpopulations (Fig. 4c) the KeyGenes algorithm23 was used with a panel of pre-/peri-implantation tissues21 as training set. Since there were 3 replicates per tissue in the training set, leave-one-out cross-validation had to be used instead of 10-fold cross-validation.

Expression in the long term cultured ectoderm- and XEN-like cells was compared to these tissue expression data sets from the literature: neural progenitor cells64, neural crest cells65, yolk sac66 and a XEN cell line67.

Differential expression between mESCs and ectoderm-like or XEN-like cells at 96 h (Supplementary Fig. 4e) was identified by an MA-plot based method using biological replicates for all three conditions68.

### Computational analysis SCRB-seq experiments

A histogram of the total number of UMIs detected per cell is shown in Supplementary Fig. 2a. To reduce the influence of technical noise we discarded cells with less than 2000 UMIs (red vertical line in Supplementary Fig. 2). This cutoff nearly minimized the upper bound of the counting error per gene (Supplementary Fig. 2b) estimated by

$$\varepsilon = \frac{1}{{\sqrt { < {\rm{UMI}} >} }} \cdot \frac{1}{{\sqrt {\# {\rm{cells}}} }}$$

while not significantly reducing the number of detected genes (13,720, Supplementary Fig. 2c)—defined as the number of genes, which had more than one UMI in more than one cell. Due to this cutoff 2451 out of 3456 measured cells were used for further analysis (Supplementary Fig. 2e), where these are the numbers of cells analyzed at each time point:

In individual cells with more than 2000 total UMIs 850 genes were detected on average.

For all further analyses, except the calculation of Fano factors, the data was normalized in the following way to account for differences in efficiency of transcript recovery between wells: UMI counts were divided by the total number of UMI counts per cell and then multiplied by the median of total UMI counts across all cells growing in 2i/L. For the calculation of Fano factors (Fig. 1d) UMI counts were down-sampled to 2000 UMI counts per cell. This down-sampling procedure ensured that the contribution of counting error to the Fano factors was equal for all cells from all time points. To include only those genes, which exhibited significant, biological variability, we considered the coefficient of variation (CV) of individual genes over all time points with respect to the mean expression level as well as the CVs of ERCC spike-ins with known abundance (Supplementary Fig. 2f). The increase in variability with decreased average expression reflected higher technical and counting noise for lowly expressed genes. We used the 829 genes, which had the 5% highest ratios of CV and the moving average of the CV for principal component analysis, k-means clustering and t-SNE mapping (see below).

To further characterize the performance of SCRB-seq we first compared SCRB-seq data averaged over cells for individual time points with bulk RNA-seq and found them to be strongly correlated (Supplementary Fig. 2g, Pearson correlation ρ = 0.75). We compared 100 randomly selected pairs of cells growing in 2i/L and found that SCRB-seq measurements of individual cells were strongly correlated (Supplementary Fig. 2h, Pearson correlation ρ = 0.63). By analysis of UMI counts of ERCC spike-in RNA we determined that UMI counts scaled approximately linearly with the spiked-in transcripts—the slope of a linear fit to the log-log plot of spike ins vs. UMI counts was 0.78. The efficiency of transcript recovery as determined from the offset of that linear fit was about 0.9% (Supplementary Fig. 2i).

For principal component analysis (PCA) we considered genes, which belonged to the upregulated clusters (clusters 5 and 6, Supplementary Fig. 1a) and were among the most variable genes (Supplementary Fig. 2f). Prior to PCA expression profiles of individual genes were converted to z-scores using the average expression over all time points and the moving average of the coefficient of variation (Supplementary Fig. 2f) to preserve biological variability. PCA was performed with all cells across all time points and expression profiles of individual cells were then projected on the principal components thusly determined. The genes with the highest loadings in the first two principal components are listed in Supplementary Fig. 3a and their loadings are represented graphically in Supplementary Fig. 3b.

To discover clusters of cells we used k-means clustering including all 829 most variable genes with (1—Pearson correlation) as the distance metric. Cluster-wise assessment of stability12 was used to determine the robustness of clusters. In particular, we calculated the Jaccard similarities between clusters found in bootstrapped samples. Clusterings resulting in Jaccard similarities close to 0.5 were considered unstable. In this way two stable clusters were found for the 96 h time point. For earlier time point cells were classified according to similarity with the clusters found at 96 h or mESCs at 0 h. In particular, we first calculated the mean expression profiles of mESCs, as well as the XEN-like and ectoderm-like subpopulations at 96 h. Then Pearson correlation was calculated between those average profiles and expression profiles of individual cells at earlier time points. A cell was classified as a particular cell type when the correlation with this particular cell type exceeded the correlation with all other cell types.

Gene expression of individual genes in the SCRB-seq data set was represented in color by normalizing to the maximum expression per time point, linear histogram stretching (1st to 99th percentile) and subsequent linear mapping to a custom colormap (Supplementary Fig. 9d).

For t-distributed stochastic neighbor embedding (t-SNE) we considered genes, which were among the most variable genes (Supplementary Fig. 2f). Prior to t-SNE mapping profiles of individual genes were converted to z-scores using the average expression over all time points and the moving average of the coefficient of variation (Supplementary Fig. 2f) to preserve biological variability. One-dimensional t-SNE maps were computed using the MATLAB Toolbox for Dimensionality Reduction (v0.8.1—March 2013) (18, L.J.P. van der Maaten, http://lvdmaaten.github.io/drtoolbox/). Expression of Rex1 was represented in color by normalizing to the maximum expression, linear histogram stretching (0th to 95th percentile) and subsequent linear mapping to a custom colormap.

The Fano factor reported in Fig. 1d is the Fano factor of individual genes averaged over all significantly variable genes. For this calculation down-sampled SCRB-seq data was used (see above). To determine the Fano factor of possible subpopulations, cells were first clustered by k-means clustering. Then the Fano factor was calculated separately for each cluster and averaged over all clusters. This procedure was carried out repeatedly and the resulting Fano factors were again averaged.

Hierarchical clustering of the SCRB-seq data (Supplementary Fig. 3) was performed using standard MATLAB routines. The particular clustering method was complete-linkage clustering using (1—Pearson correlation) as the distance metric.

For pseudotime ordering of cells a correlation-based pseudotime was defined by τ = Rpluri—0.5*(Rect + Rxen), where are Rpluri, Rect and Rxen are the Pearson correlations of an individual expression profile with the average expression of mESCs, ectoderm-like cells at 96 h and XEN-like cells at 96 h, respectively.

For comparison with expression in in vivo tissue21 we performed PCA on standardized in vivo data using the 829 most variable genes defined above. We then projected our standardized SCRB-seq data on the plane spanned by the first two principal components.

### Computational analysis SMART-seq2 experiments

Only cells with at least 200000 reads per cell were used, resulting in the following numbers of cells analyzed at the respective time points:

For all further analyses the data was normalized in the following way to account for differences in the total number of reads between samples: RPKM for individual genes were divided by the total number of RPKM per cell and then multiplied by the median of total RPKM across all cells growing in 2i/L. Cells with high expression of Cd24 or Pdgfra were classified as shown in Supplementary Fig. 9a. Out of the 82 cells measured by SMART-seq2 at 48 h, 10 were considered XEN-like (Pdgfra high) and 29 ectoderm-like (Cd24 high). To compute p-values for gene expression differences in these subpopulations we used a null model that assumes that all cells were essentially identical and gene expression differences were only due to biological and technical noise. We repeatedly sampled 10 or 29 cells, respectively, with replacement from the pool of cells which did not express Pdgfra or Cd24 and calculated the average expression level for each gene. The distribution of average expression levels for each gene thusly obtained was then fit with a normal distribution. The p-value was then calculated using this normal distribution and the average expression level observed in the Cd24 or Pdgfra high cells.

To account for multiple hypothesis testing we used the Benjamini-Hochberg procedure and set the false discovery rate to 0.05. Additionally, we required a minimal fold-change of 2 and an absolute expression level bigger than 1 normalized RPKM for a gene to be accepted as differentially expressed Finally, we considered only genes which were defined as transcriptional regulators by gene ontology (GO) term annotation (GO:0003700, GO:0044212, GO:0045944, GO:0006355, GO:0000981). We considered a gene to be expressed / not expressed in the pluripotent state when it was robustly expressed (normalized RPKM > 5) in at least 50% / less than 5% of the cells at 0 h.

We combined the transcriptional regulators identified in this way with pluripotency network factors51 to arrive at a set of transcription factors which are likely relevant for the lineage decision studied here. For the calculation of co-expression (Fig. 6b, c) we considered a gene to be expressed at normalized RPKM values over 1.

### Single-molecule FISH

Cells growing in gelatinized tissue culture dishes were washed twice with PBS, detached with Accutase (Life technologies) and resuspended in full medium. Formaldehyde was added to the cell suspension to a final concentration of 4%. Cells were incubated for 12 min at room temperature while being rotated and then spun down for 3 min at 90 x g. Subsequently cells were permeabilized at least over night in 70% ethanol. For hybridization and imaging cells were attached to chambered cover slides (Nunc Lab-Tek) coated with poly-l-lysine.

In the case of intact colonies, adherent cells were fixed for 15 min with 4% formaldehyde by adding formaldehyde to the growth medium and subsequently permeabilized in 70% ethanol.

Oligonucleotide libraries with 20-nt probes for Nanog, Gbx2, Tbx3, Gata6 and Pax6 were designed and fluorescently labeled as previously described40. Briefly, a home-made MATLAB script was used to design probes with close to 45% GC content. The probes were then checked for low-complexity sequences and binding to other than the desired transcript by BLAST (https://blast.ncbi.nlm.nih.gov/Blast.cgi). Suitable probes were ordered from Biosearch Technologies and fluorescently labeled using a 3’-amine modification of the oligos and amine reactive dyes (GE healthcare). The hybridization buffer used for smFISH contained 2 × SSC buffer, 25 or 40% formamide, 10% Dextran Sulphate (Sigma), E. coli tRNA (Sigma), Bovine Serum Albumin (Ambion) and Ribonucleoside Vanadyl Complex (New England biolabs). 50–75 ng of the desired probes were used per 100 µl of hybridization buffer. (The mass refers only to pooled oligonucleotides, excluding fluorophores, and is based on absorbance measurements at 260 nm). Probes were hybridized for 16–18 h at 30 °C, after which we washed cells twice for 30 min at 30 °C in wash buffer (2 × SSC, 25% formamide (for all probes except Gbx2 and Tbx3) or 40% formamide (for Gbx2 and Tbx3)), supplemented with Hoechst 33,342. For microscopy, we filled the hybridization chamber with a mounting solution containing 1 x PBS, 0.4% Glucose, 100 μg/ml Catalase, 37 μg/ml Glucose Oxidase, and 2 mM Trolox. Imaging was done exactly as described previously69. Images were taken on a NIKON Ti-E inverted fluorescence microscope equipped with a Roper scientific PIXIS 1024B camera and a 100x oil immersion objective (numerical aperture 1.49). Custom filters (Omega Optical) were used for imaging TMR and Alexa 594 and a standard filter (Chroma) for Cy5. Exposure times ranges between 1 and 3 s and the distance between planes in a z-stack was 0.3 μm Home-made MATLAB scripts were used for image analysis. Cells positive for one of the assayed genes were classified as shown in Supplementary Fig. 10b.

### Quantification of the flow cytometry experiments

The distribution of cells in the space of CD24 and PDGFRA expression was modeled by the sum of 4 bivariate normal distributions. This model has in principle 19 free parameters (8 for the means, 8 for the standard deviations and 3 for the size of the relative contributions). To ensure robust fitting to the date we reduced the number of parameters to 9 by keeping the standard deviations constant and only allowing 4 different values for the means.

$$\begin{array}{l}\\ p\left( {x,y;{f_1},{f_2},{f_3},\mu _x^{low},\mu _y^{low},\mu _x^{high},\mu _y^{high},{\sigma _x},{\sigma _y}} \right)\\ = {f_1} \cdot N\left( {x,y;\mu _x^{low},\mu _y^{low},{\sigma _x},{\sigma _y}} \right) + {f_2} \cdot N\left( {x,y;\mu _x^{low},\mu _y^{high},{\sigma _x},{\sigma _y}} \right)\\ + {f_3} \cdot N\left( {x,y;\mu _x^{high},\mu _y^{low},{\sigma _x},{\sigma _y}} \right) + \left( {1 - \mathop {\sum }\limits_{i = 1}^3 {f_i}} \right) \cdot N\left( {x,y;\mu _x^{high},\mu _y^{high},{\sigma _x},{\sigma _y}} \right)\\ \end{array}$$

N(x, y, µ x, µ y, σ x, σ y) is a bivariate normal distribution in x and y (PDGFRA and CD24 expression, respectively) with mean (µ x, µ y) and standard deviation (σ x, σ y). This model was fit to a reference data set (typically untreated control cells after 96 h of RA exposure) by maximizing the log-likelihood −log(p). To subsequently calculate the size of the fractions fi for a particular sample we first calculated the probabilities that the expression values (x, y) found in a particular cell were drawn from one the 4 normal distributions N(x, y, µ x, µ y, σ x, σ y). The cell was then ascribed to the distribution from which it was most likely drawn.

### Stochastic simulation of the lineage transition

We simulated the differentiation process using a discretized version of the Langevin equation describing the system (Euler method):

$$\begin{array}{l}\\ dX = \left( {{a_X}\frac{{{X^n}}}{{{\theta ^n} + {X^n}}} + b\frac{{{\theta ^n}}}{{{\theta ^n} + {E^n}}} - kX} \right){\rm{\Delta }} + \sqrt {D{\rm{\Delta }}} \,{\cal N}\left( {0,1} \right)\\ dE = \left( {{a_E}\frac{{{E^n}}}{{{\theta ^n} + {E^n}}} + b\frac{{{\theta ^n}}}{{{\theta ^n} + {X^n}}} - kE} \right){\rm{\Delta }} + \sqrt {D{\rm{\Delta }}} \,{\cal N}\left( {0,1} \right)\\ \end{array}$$

X and E indicate the expression levels of the XEN and ectoderm programs respectively. N(0, 1) indicates a Wiener process with mean 0 and standard deviation 1. D sets the strength of gene expression noise and Δ determines the size of the time step. After initializing X and E randomly between 0 and 0.1 we first equilibrated the system for 100 iterations. Subsequently, we propagated the system for 200 additional iterations. To relate the simulation to experimental time scales, the end point of the simulation was taken to be at 96 h. To model the exit from pluripotency the degradation parameter k was switched from a high value (k = 10) to a low value (k = 1) after 12 h (25 iterations), which allowed X and E to increase. To model timed application of RA the auto-activation parameter for the XEN program aX was switched at various points in time (no RA: a X = 0; RA: a X = 0.5). For each condition we generated 10,000 trajectories and counted the number of trajectories that ended at the XEN or ectoderm attractor (see Supplementary Fig. 8b). The relative frequency of trajectories ending at the XEN attractor is reported in Fig. 5e.

### Used parameters

n = 4

θ = 0.5

Δ = 0.05

pluripotency: k = 10

differentiation: k = 1

a E = 0.5

no RA:, a X = 0

RA: a X = 0.5

In Supplementary Fig. 8b

low noise: D = 0.0001

high noise: D = 0.01

### Code availability

The MATLAB scripts used for data analysis and simulations are freely available on request from the corresponding author.

### Data availability

All raw and processed data is freely available from the GEO repository (https://www.ncbi.nlm.nih.gov/geo/) under accession number GSE79578.

## References

1. 1.

Cohen, D. E. & Melton, D. Turning straw into gold: directing cell fate for regenerative medicine. Nat. Rev. Genet. 12, 243–252 (2011).

2. 2.

Soldner, F. & Jaenisch, R. iPSC Disease Modeling. Science 338, 1155–1156 (2012).

3. 3.

Tabar, V. & Studer, L. Pluripotent stem cells in regenerative medicine: challenges and recent progress. Nat. Rev. Genet. 15, 82–92 (2014).

4. 4.

Semrau, S. & van Oudenaarden, A. Studying lineage decision-making in vitro: emerging concepts and novel tools. Annu. Rev. Cell Dev. Biol. 31, 317–345 (2015).

5. 5.

Balázsi, G., van Oudenaarden, A. & Collins, J. J. Cellular decision making and biological noise: from microbes to mammals. Cell 144, 910–925 (2011).

6. 6.

Pauklin, S. & Vallier, L. The cell-cycle state of stem cells determines cell fate propensity. Cell 155, 135–147 (2013).

7. 7.

Kobayashi, T. et al. The cyclic gene Hes1 contributes to diverse differentiation responses of embryonic stem cells. Gene Dev. 23, 1870–1875 (2009).

8. 8.

Morgani, S. M. & Brickman, J. M. The molecular underpinnings of totipotency. Philos. Trans. R. Soc. Lond. B. 369, 20130549–20130549 (2014).

9. 9.

Gudas, L. J. & Wagner, J. A. Retinoids regulate stem cell differentiation. J. Cell Physiol. 226, 322–330 (2010).

10. 10.

Rhinn, M. & Dolle, P. Retinoic acid signalling during development. Development 139, 843–858 (2012).

11. 11.

Ying, Q.-L. et al. The ground state of embryonic stem cell self-renewal. Nature 453, 519–523 (2008).

12. 12.

Hennig, C. Cluster-wise assessment of cluster stability. Comput. Stat. Data. Anal. 52, 258–271 (2007).

13. 13.

Capo-chichi, C. D. et al. Perception of differentiation cues by GATA factors in primitive endoderm lineage determination of mouse embryonic stem cells. Dev. Biol. 286, 574–586 (2005).

14. 14.

Bain, G., Ray, W. J., Yao, M. & Gottlieb, D. I. Retinoic acid promotes neural and represses mesodermal gene expression in mouse embryonic stem cells in culture. Biochem. Biophys. Res. Commun. 223, 691–694 (1996).

15. 15.

Soumillon, M., Cacchiarelli, D. & Semrau, S. Characterization of directed differentiation by high-throughput single-cell RNA-Seq. bioRxiv 003236 (2014).

16. 16.

Pruszak, J., Ludwig, W., Blak, A., Alavian, K. & Isacson, O. CD15, CD24, and CD29 define a surface biomarker code for neural lineage differentiation of stem cells. Stem Cells 27, 2928–2940 (2009).

17. 17.

Artus, J., Panthier, J.-J. & Hadjantonakis, A.-K. A role for PDGF signaling in expansion of the extra-embryonic endoderm lineage of the mouse blastocyst. Development 137, 3361–3372 (2010).

18. 18.

van der Maaten, L. & Hinton, G. Visualizing Data using t-SNE. J. Mach. Learn. Res. 9, 2579–2605 (2008).

19. 19.

Marks, H. et al. The Transcriptional and Epigenomic Foundations of Ground State Pluripotency. Cell Stem. Cell 149, 590–604 (2012).

20. 20.

Trapnell, C. et al. The dynamics and regulators of cell fate decisions are revealed by pseudotemporal ordering of single cells. Nat. Biotechnol. 32, 381–386 (2014).

21. 21.

Boroviak, T. et al. Lineage-specific profiling delineates the emergence and progression of naive pluripotency in mammalian embryogenesis. Dev. Cell 35, 366–382 (2015).

22. 22.

Boroviak, T., Loos, R., Bertone, P., Smith, A. & Nichols, J. The ability of inner-cell-mass cells to self-renew as embryonic stem cells is acquired following epiblast specification. Nat. Cell Biol. 16, 516–528 (2014).

23. 23.

Roost, M. S. et al. KeyGenes, a tool to probe tissue differentiation using a human fetal transcriptional atlas. Stem Cell Reports 4, 1112–1124 (2015).

24. 24.

Shakiba, N. et al. CD24 tracks divergent pluripotent states in mouse and human cells. Nat. Commun. 6, 1–11 (2015).

25. 25.

Malleshaiah, M. et al. Nac1 coordinates a sub-network of pluripotency factors to regulate embryonic stem cell differentiation. Cell Reports 14, 1181–1194 (2016).

26. 26.

Hamilton, W. B. & Brickman, J. M. Erk signaling suppresses embryonic stem cell self-renewal to specify endoderm. Cell Reports 9, 2056–2070 (2014).

27. 27.

Schröter, C., Rué, P., Mackenzie, J. P. & Arias, A. M. FGF/MAPK signaling sets the switching threshold of a bistable circuit controlling cell fate decisions in embryonic stem cells. Development 142, 4205–4216 (2015).

28. 28.

Kunath, T. et al. FGF stimulation of the Erk1/2 signalling cascade triggers transition of pluripotent embryonic stem cells from self-renewal to lineage commitment. Development 134, 2895–2902 (2007).

29. 29.

Ohnishi, Y. et al. Cell-to-cell expression variability followed by signal reinforcement progressively segregates early mouse lineages. Nat. Cell Bio. 16, 27–37 (2014).

30. 30.

JOHNSON, A. T. et al. Synthesis and characterization of a highly potent and effective antagonist of retinoic acid receptors. J. Med. Chem. 38, 4764–4767 (1995).

31. 31.

Ying, Q.-L., Stavridis, M., Griffiths, D., Li, M. & Smith, A. Conversion of embryonic stem cells into neuroectodermal precursors in adherent monoculture. Nat. Biotechnol. 21, 183–186 (2003).

32. 32.

Huang, S., Guo, Y.-P., May, G. & Enver, T. Bifurcation dynamics in lineage-commitment in bipotent progenitor cells. Dev. Biol. 305, 695–713 (2007).

33. 33.

Picelli, S. et al. Smart-seq2 for sensitive full-length transcriptome profiling in single cells. Nat. Methods. 10, 1096–1098 (2013).

34. 34.

Liu, A. & Joyner, A. L. EN and GBX2 play essential roles downstream of FGF8 in patterning the mouse mid/hindbrain region. Development 128, 181–191 (2001).

35. 35.

Osumi, N., Shinohara, H., Numayama-Tsuruta, K. & Maekawa, M. Concise review: pax6 transcription factor contributes to both embryonic and adult neurogenesis as a multifunctional regulator. Stem Cells 26, 1663–1672 (2008).

36. 36.

Lu, R., Yang, A. & Jin, Y. Dual Functions of T-Box 3 (Tbx3) in the control of self-renewal and extraembryonic endoderm differentiation in mouse embryonic stem cells. J. Biol. Chem. 286, 8425–8436 (2011).

37. 37.

Arnold, S. J. & Robertson, E. J. Making a commitment: cell lineage allocation and axis patterning in the early mouse embryo. Nat. Rev. Mol. Cell. Biol. 10, 91–103 (2009).

38. 38.

Ballim, R. D., Mendelsohn, C., Papaioannou, V. E. & Prince, S. The ulnar-mammary syndrome gene, Tbx3, is a direct target of the retinoic acid signaling pathway, which regulates its expression during mouse limb development. Mol. Biol. Cell. 23, 2362–2372 (2012).

39. 39.

Balmer, J. E. Gene expression regulation by retinoic acid. J. Lipid Res. 43, 1773–1808 (2002).

40. 40.

Raj, A., van den Bogaard, P., Rifkin, S. A., van Oudenaarden, A. & Tyagi, S. Imaging individual mRNA molecules using multiple singly labeled probes. Nat. Methods. 5, 877–879 (2008).

41. 41.

Klein, A. M. et al. Droplet barcoding for single-cell transcriptomics applied to embryonic stem cells. Cell 161, 1187–1201 (2015).

42. 42.

Morgani, S. M. et al. Totipotent embryonic stem cells arise in ground-state culture conditions. Cell Reports 3, 1945–1957 (2013).

43. 43.

Saiz, N., Williams, K. M., Seshan, V. E. & Hadjantonakis, A.-K. Asynchronous fate decisions by single cells collectively ensure consistent lineage composition in the mouse blastocyst. Nat. Commun. 7, 13463 (2016).

44. 44.

Chazaud, C., Yamanaka, Y., Pawson, T. & Rossant, J. Early lineage segregation between epiblast and primitive endoderm in mouse blastocysts through the Grb2-MAPK pathway. Dev. Cell 10, 615–624 (2006).

45. 45.

Yamanaka, Y., Lanner, F. & Rossant, J. FGF signal-dependent segregation of primitive endoderm and epiblast in the mouse blastocyst. Development 137, 715–724 (2010).

46. 46.

Turner, D. A., Trott, J., Hayward, P., Rue, P. & Martinez Arias, A. An interplay between extracellular signalling and the dynamics of the exit from pluripotency drives cell fate decisions in mouse ES cells. Biology Open 3, 614–626 (2014).

47. 47.

Nichols, J., Silva, J., Roode, M. & Smith, A. Suppression of Erk signalling promotes ground state pluripotency in the mouse embryo. Development 136, 3215–3222 (2009).

48. 48.

Thomson, M. et al. Pluripotency factors in embryonic stem cells regulate differentiation into germ layers. Cell 145, 875–889 (2011).

49. 49.

Trott, J. & Martinez Arias, A. Single cell lineage analysis of mouse embryonic stem cells at the exit from pluripotency. Biology Open 2, 1049–1056 (2013).

50. 50.

Marco, E. et al. Bifurcation analysis of single-cell gene expression data reveals epigenetic landscape. Proc. Natl Acad. Sci. 111, E5643–E5650 (2014).

51. 51.

Dunn, S. J., Martello, G., Yordanov, B., Emmott, S. & Smith, A. G. Defining an essential transcription factor program for naive pluripotency. Science 344, 1156–1160 (2014).

52. 52.

Kim, J., Chu, J., Shen, X., Wang, J. & Orkin, S. H. An extended transcriptional network for pluripotency of embryonic stem cells. Cell 132, 1049–1061 (2008).

53. 53.

Ivanova, N. et al. Dissecting self-renewal in stem cells with RNA interference. Nature Publishing Group 442, 533–538 (2006).

54. 54.

Niwa, H., Ogawa, K., Shimosato, D. & Adachi, K. A parallel circuit of LIF signalling pathways maintains pluripotency of mouse ES cells. Nature 460, 118–122 (2009).

55. 55.

Tai, C.-I. & Ying, Q.-L. Gbx2, a LIF/Stat3 target, promotes reprogramming to and retention of the pluripotent ground state. J. Cell Sci. 126, 1093–1098 (2013).

56. 56.

Singer, Z. S. et al. Dynamic Heterogeneity and DNA methylation in embryonic stem cells. Mol. Cell 55, 319–331 (2014).

57. 57.

Mojtahedi, M. et al. Cell fate decision as high-dimensional critical state transition. PLoS Biol. 14, e2000640 (2016).

58. 58.

Richard, A. et al. Single-cell-based analysis highlights a surge in cell-to-cell molecular variability preceding irreversible commitment in a differentiation process. PLoS Biol. 14, e1002585 (2016).

59. 59.

Fox, M. H. A model for the computer analysis of synchronous DNA distributions obtained by flow cytometry. Cytometry 1, 71–77 (1980).

60. 60.

Picelli, S. et al. Full-length RNA-seq from single cells using Smart-seq2. Nat. Protoc. 9, 171–181 (2014).

61. 61.

Trombetta, J. J. et al. Preparation of single-cell RNA-seq libraries for next generation sequencing. Curr. Protoc. Mol. Biol. 107, 4.22.1–4.22.17 (2014).

62. 62.

Trapnell, C. et al. Differential gene and transcript expression analysis of RNA-seq experiments with TopHat and Cufflinks. Nat. Protoc. 7, 562–578 (2012).

63. 63.

Huang, D. W., Sherman, B. T. & Lempicki, R. A. Systematic and integrative analysis of large gene lists using DAVID bioinformatics resources. Nat. Protoc. 4, 44–57 (2008).

64. 64.

Wang, J. et al. Transcriptome analysis of neural progenitor cells by a genetic dual reporter strategy. Stem Cells 29, 1589–1600 (2011).

65. 65.

Hagiwara, K. et al. Molecular and cellular features of murine craniofacial and trunk neural crest cells as stem cell-like cells. PLoS. ONE. 9, e84072 (2014).

66. 66.

Stumpo, D. J. et al. Deficiency of the placenta- and yolk sac-specific tristetraprolin family member ZFP36L3 identifies likely mRNA targets and an unexpected link to placental iron metabolism. Development 143, 1424–1433 (2016).

67. 67.

Julio, M. K.-D. et al. Regulation of extra-embryonic endoderm stem cell differentiation by Nodal and Cripto signaling. Development 138, 3885–3895 (2011).

68. 68.

Wang, L., Feng, Z., Wang, X., Wang, X. & Zhang, X. DEGseq: an R package for identifying differentially expressed genes from RNA-seq data. Bioinformatics 26, 136–138 (2010).

69. 69.

Semrau, S. et al. FuseFISH: robust detection of transcribed gene fusions in single cells. Cell Reports 6, 18–23 (2014).

70. 70.

Buecker, C. et al. Reorganization of enhancer patterns in transition from naive to primed pluripotency. Cell Stem. Cell 14, 838–853 (2014).

## Acknowledgements

S.S. was supported by the Netherlands Organisation for Scientific Research (NWO/OCW), as part of the Frontiers of Nanoscience program and by an NWO Rubicon award. J.G. was supported by the Boehringer Ingelheim Foundation as well as a Jerome and Florence Brill Graduate Student Fellowship. R.J. was supported by NIH grants HD 045022 and RO1-CA084198. A.v.O. was supported by an ERC Advanced grant (ERC-AdG 294325-GeneNoiseControl) and an NWO Vici award.

## Author information

Authors

### Contributions

S.S. designed and performed experiments, analyzed data and wrote the manuscript, J.G. performed experiments and edited the manuscript, M. S. performed experiments, T.S.M. supervised the project, R.J. and A.v.O. supervised the project and edited the manuscript.

### Corresponding author

Correspondence to Stefan Semrau.

## Ethics declarations

### Competing interests

The authors declare no competing financial interests.

Publisher's note: Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

## Rights and permissions

Reprints and Permissions

Semrau, S., Goldmann, J.E., Soumillon, M. et al. Dynamics of lineage commitment revealed by single-cell transcriptomics of differentiating embryonic stem cells. Nat Commun 8, 1096 (2017). https://doi.org/10.1038/s41467-017-01076-4

• Accepted:

• Published:

• ### Investigating higher-order interactions in single-cell data with scHOT

• Shila Ghazanfar
• , Yingxin Lin
• , Xianbin Su
• , David Ming Lin
• , Ellis Patrick
• , Ze-Guang Han
• , John C. Marioni
•  & Jean Yee Hwa Yang

Nature Methods (2020)

• ### Epigenetic regulator function through mouse gastrulation

• Stefanie Grosswendt
• , Helene Kretzmer
• , Zachary D. Smith
• , Abhishek Sampath Kumar
• , Sara Hetzel
• , Lars Wittler
• , Sven Klages
• , Bernd Timmermann
• , Shankar Mukherji
•  & Alexander Meissner

Nature (2020)

• ### Selective USP7 inhibition elicits cancer cell killing through a p53-dependent mechanism

• Nathan J. Schauer
• , Xiaoxi Liu
• , Robert S. Magin
• , Laura M. Doherty
• , Wai Cheung Chan
• , Scott B. Ficarro
• , Wanyi Hu
• , Rebekka M. Roberts
• , Roxana E. Iacob
• , Björn Stolte
• , Andrew O. Giacomelli
• , Sumner Perera
• , Kyle McKay
• , Sarah A. Boswell
• , Ellen L. Weisberg
• , Arghya Ray
• , Dharminder Chauhan
• , Sirano Dhe-Paganon
• , Ken C. Anderson
• , James D. Griffin
• , Jianing Li
• , William C. Hahn
• , Peter K. Sorger
• , John R. Engen
• , Kimberly Stegmaier
• , Jarrod A. Marto
•  & Sara J. Buhrlage

Scientific Reports (2020)

• ### SCRABBLE: single-cell RNA-seq imputation constrained by bulk RNA-seq data

• Tao Peng
• , Qin Zhu
• , Penghang Yin
•  & Kai Tan

Genome Biology (2019)

• ### Computational analysis of single-cell transcriptomics data elucidates the stabilization of Oct4 expression in the E3.25 mouse preimplantation embryo

• Daniela Gerovska
•  & Marcos J. Araúzo-Bravo

Scientific Reports (2019)

format: plain


text:
# How do you evaluate (- 7k + 2s - 3q ) - ( - 2q - 3k ) + ( 4k + q - 2s )?

Jan 30, 2018

0

#### Explanation:

Let’s start by making the equation to addition format (as it is easier to evaluate).

Note: Negative x Negative = Positive

So $\left(- 7 k + 2 s - 3 q\right) - \left(- 2 q - 3 k\right) + \left(4 k + q - 2 s\right)$

Becomes $\left(- 7 k + 2 s - 3 q\right) + \left(2 q + 3 k\right) + \left(4 k + q - 2 s\right)$

Now we can group all the like terms together; move all the same variables together. Ill arrange it in brackets for you to get a better visual.

$\left(- 7 k + 2 s - 3 q\right) + \left(2 q + 3 k\right) + \left(4 k + q - 2 s\right)$

$\left(- 7 k + 3 k + 4 k\right) + \left(2 s - 2 s\right) + \left(- 3 q + 2 q + q\right)$

$\left(- 7 k + 7 k\right) + \left(2 s - 2 s\right) + \left(- 3 q + 3 q\right)$

$0 + 0 + 0 = 0$


text:
# Homework Help: Heat Thermodynamics Questions

1. Jan 25, 2010

### mmmboh

Alright so my friends and I have been trying really hard on our new assignment, and there are some problems that we still can't do. It's like no matter how hard we try the book doesn't explain it well enough for us to be confident in our answers. So Pleeease help!

1. A mixture of hydrogen and oxygen is enclosed in a rigid insulating container and exploded by a spark. The temperature and pressure both increase. Neglect the small amount of energy provided by the spark itself.
A) Has there been a flow of heat into the system.
B)Has any work been done by the system.
C)Has there been any change in internal energy U of the system?

2.The water in a rigid. insulated cylindrical tank is set in rotation and left to itself. It is eventually brought to rest by viscous forces. The tank and the water constitute the system.
A)Is any work done during the process in which the water is brought to rest?
B)Is there a flow of heat?
C)Is there any change in the internal energy U?

3.Compressing the system represented in the figure along the adiabatic path a-c requires 1000J. Compressing the system along b-c requires 1500J but 600J of heat flow out of the system.
A)Calculate the work done, the heat absorbed, and the internal energy change of the system in each process and in the total cycle a-b-c-a.
B)Sketch this cycle on a P-V diagram.
C)What are the limitations on the values that could be specified for process b-c given that 1000J are required to compress the system along a-c

1.A)There is no heat flow because the system is insulated.
B)Yes adiabatic work is done. (I'm not sure why though)
C)Yes the change in internal energy is equal to the adiabatic work done (I think it's an increase in energy.

2.A)Yes there is dissipative work caused by the viscous forces.
B)There is no flow of heat because the system is insulated.
C)Yes the change in internal energy equals the dissipative work (decrease in energy?)

3.A)For a-c I wrote there is 1000J work done, no heat flow, and 1000J change in energy.
For b-c there is 1500J work done, 600J heat flow, 2100J change in energy.
For a-b-c-a there is no heat absorbed, 600J work done, no change in internal energy.
I haven't yet sketched the diagram but I am not sure what to do for the limitations on the values that could be specified for b-c given that 1000J are required to compress the system along a-c/

Diagram for 3:

Thanks for the help :D

2. Jan 25, 2010

### mmmboh

Anyone?

3. Jan 26, 2010

### Jobrag

My old mechanics teacher always used to tell us "draw a box round it" , look at your problems draw a boundry round them and ask "does anything cross the boundry?"

4. Jan 26, 2010

### mathman44

Can anyone help clarify these?

5. Jan 26, 2010

### Jobrag

Take your first question.
Draw a boundry round the the container, does anything pass across the boundry?
If nothing crosses the boundry going in no heat can have entered the system.
If nothing crosses the boundry going out the system has not done any work.
If the system has performed no work and no heat has entered the system can the internal energy have changed?

6. Jan 26, 2010

### mmmboh

Ok got it thanks a lot!. I figured out the third question, but in the case of the "boundary", that would make my second answer correct right?

7. Jan 26, 2010

### fluidistic

I'm confused about this. Because the problem states that there is an increase of temperature. Plus, the internal energy of an ideal gas is directly proportional to the temperature. Hence an increase of temperature should mean an increase of internal energy. But as you implied, according to the first law of thermodynamics, it cannot be. So what's going on?

8. Jan 26, 2010

### Jobrag

fuidistic
Heat isn't the only form of energy.

9. Jan 26, 2010

### fluidistic

Hmm I know (for instance work), but how does it matter here?
I don't want to hijack the thread, so if you want to reply me by PM, it's ok.

10. Jan 26, 2010

### Jobrag

fluidistic
If you accept that the thermal energy of the system has increased (the gas has got hotter) and that there has been no energy imported into the system, what else might have happened.
Work is heat and heat is work is a true statement but NOT excusivly true.

11. Jan 27, 2010

### fluidistic

I accept that the thermal energy (temperature) of the system has increased and that no energy has been imported into the system. It is explicitly said that
I also agree that there is no work done since the volume is constant. Hence one might say that the variation of internal energy of the system is 0, according to the first law.
My problem is looking at the definition of the internal energy of an ideal gas: it depends on temperature. So a hotter gas means a gas with a greater internal energy.
But how can that be? It seems to be in contradiction with the first law of thermodynamics here.
Do you understand where I'm misunderstanding?

12. Jan 27, 2010

### Jobrag

fluiodistic
What happens when hydrogen and oxygen are mixed and a spark applied?

13. Jan 27, 2010

### fluidistic

It forms water if I remember well. (from high school since I've no chemistry in my physics degree).
Hmm ok... But to be rigorous, don't you think we'd have to show that it is under the liquid form? Because if it would be a gas, it seems the first law of thermodynamics couldn't be applied. Is it sufficient to affirm that it is under the liquid form otherwise the first law wouldn't be right?

14. Jan 27, 2010

### Jobrag

It forms water but what is given off during the reaction?

15. Jan 27, 2010

### fluidistic

I don't remember. Heat I guess.

16. Jan 27, 2010

### mmmboh

Hydrogen and oxygen form water and release energy. Having said that I hope my second answer is right...

17. Jan 27, 2010

### Jobrag

OK lets wrap this up. The hydrogen and oxygen turn into water in its gaseous phase (steam). The energy released by the chemical reaction provides the heat to raise the temperature and pressure, nothing needs to enter the system as energy is changed from one state to another.

18. Jan 28, 2010

### fluidistic

Ok thank you very much. So as in this (https://www.physicsforums.com/showthread.php?p=2552092#post2552092) thread says, there is the chemical energy that Resnick-Halliday never mentioned.

19. Jan 28, 2010

### Andrew Mason

A) There is heat flow into the system. The post-explosion system is definitely at a higher temperature than the pre-explosion system (higher internal energy) and there is no work done (volume remains constant). Apply the first law to see that Q cannot be 0: $\Delta Q = \Delta U + W$

B) If there is no change in volume, how can there be any work done?\\

C) The change in energy is equal to the heat flow.

AM

20. Jan 28, 2010

### fluidistic

Dear Andrew Mason,
I think you are wrong by saying that there was a heat flow. I don't know if you have read the previous posts, in case of yes, please correct Jobrag and I. And I also suggest you to read https://www.physicsforums.com/showthread.php?p=2552092#post2552092 where Matterrave suggest to use a more generalized form of the firs law to explain the increase of temperature of the gas when both the work done by the gas and the heat flow are null.


text:
Calculate the energy change = energy in – energy out. Since the maximum value of the bond energy occurs when the particles are widely separated, and because of the way the pair-wise potential is defined, the bond energy of liquids and solids must be less than zero; that is, the bond energy is negative. Let us consider a structure of four atoms keeping them as close as possible to each other, as shown below. You can calculate the energy change in a reaction using average bond energies. Therefore: $E_{bond}=-3\times PE(r_o)=3\times(-0.6\times 10^{-21})=-1.8\times 10^{-21}J$.

Guarantee the serial number you enter is valid. The only task is to calculate the total number of nearest neighbors. What matters is that that the change in bond energy when bonds are broken is positive, which is the case here since we start with a negative bond energy and end up with zero bond energy. Of course, here we can just count the pairs directly (there are three! An example of a common one, known a the face-centered cubic (fcc) structure has 12 nearest neighbors as shown in the figure below. However, it is no longer possible to arrange four atoms in two-dimensions so all of them are at the $$r_o$$ separation. Unless otherwise noted, LibreTexts content is licensed by CC BY-NC-SA 3.0.

[ "article:topic", "authorname:ucd7", "license:ccby", "showtoc:no" ], $PE_{LJ}(r=1.94\sigma)=4\varepsilon\Big [\Big(\frac{\sigma}{1.94\sigma}\Big)^{12}-\Big(\frac{\sigma}{1.94\sigma}\Big)^6 \Big]=-0.074\varepsilon$, By convention, all pair-wise potentials are defined to be zero when the particles are separated sufficiently so that the force acting between the particles is zero and negative when the particles are bound. Therefore, the bond energy of any condensed substance (solid of liquid) is always negative. Include all pair-wise interactions.

However, if you plan to save an inventory of bonds, you may want to enter serial numbers.) We can think of this quantity as the change of bond energy of the two particle system initially at equilibrium. To calculate bond energy Add together the bond energies for … Another method is used to calculate the approximate bond length. We discussed in the previous section that as long as $$E_{tot}$$ is less than zero the particles remain bound. Hydrogen and chlorine react to form hydrogen chloride gas: The energy change is negative, due to the fact that the energy released by the bonds formed is greater than the energy absorbed by the bonds broken. It is possible to read $$PE(r=2.24\times 10^{-10}m)$$ from the plot, but it is more accurate to plug this values into the Lennard-Jones equation: $PE_{LJ}(r=2.24\sigma)=4\varepsilon\Big [\Big(\frac{\sigma}{2.24\sigma}\Big)^{12}-\Big(\frac{\sigma}{2.24\sigma}\Big)^6 \Big]=-0.031\varepsilon$, $E_{bond}=-2\times (1\times 10^{-21}J)-0.031\times 10^{-21} J=-2.031\times 10^{-21}J$. Legal. We call atoms that are separated by the shortest possible distance, nearest neighbors, sometimes abbreviated as nn. For Cy, $$PE(r_o)=-0.6\times 10^{-21} J$$. We saw when calculating the energy required to break just a few atoms, we set the initial energy of the structure when all atoms are motionless at equilibrium.

We will see our Particle Model of Thermal Energy that there are changes in the thermal energy at a phase change as well as in the bond energy. Energy required to break the bonds or the change in bond energy is simply the magnitude of the bond energy and is always a positive number, even though the bond energy is negative. This is sometimes hard to get our minds around. To calculate the energy required to break a bond, let us assume that the two atoms are motionless at equilibrium separation of $$r_o$$. The energy required to break a specific covalent bond in one mole of gaseous molecules is called the bond energy or the bond dissociation energy. Thus, we will neglect edge effects. Bond dissociation enthalpy and mean bond enthalpy. b) Calculate E bond for another configuration shown below.


text:
# Algorithm for solving system of integer polynomial equations

I'm looking for an algorithm to solve a system of $n$ positive integer polynomial equations, assuming that there exists one and only one solution to the system.

Could anyone point me in the right direction or literature?

Example: Assume that the number of variables is $n$, usually we expect $10 \le n \le 50$. Assume that $n=5$ for the example. We need to solve: $$x_1 = 5$$ $$x_1^2 x_2 x_3 = 100$$ $$x_2^2 x_3 + x_1 + x_4^2 = 62$$ $$x_5 x_4 + 2 x_5 = 9$$ $$x_1 + x_5 x_2 + x_3 x_4 = 21$$ $$x_1 + x_2 + x_3 + x_4 + x_5 = 17$$ Solution: From $x_1 = 5$ we get that $x_2 x_3 = 4$, thus $(x_2,x_3)=(1,4),(2,2),(4,1)$. From the third equation $x_5(x_4 + 2) = 9$, thus $(x_4, x_5) = (1, 3), (7, 1)$. From the last equation we have that $x_2 + x_3 + x_4 + x_5 = 12$, since the sum of $x_2 + x_3 = 4, 5$ and $x_4 + x_5 = 4, 8$, it means only option is $x_4 = 7, x_5 = 1, x_2 = 2, x_3 = 2$.

Note that I need an algorithm to do this , such that it is automated.

• If you do not show the system of equations - how to solve it? – individ Aug 26 '16 at 16:05
• I'm looking for a general algorithm to solve any such system. I don't have any concrete system, I can make one up, but the goal is to have a general algorithm for solving it (e.g. like for instance the Greatest Common Divider algorithm, or Gaussian Elimination) – Alex Botev Aug 26 '16 at 16:18
• For systems of linear equations - standard method of Gauss. For quadratic and slightly higher extent have to wait until about it will tell. While about it to speak. Some systems can be found there. artofproblemsolving.com/community/… It is better to write the system - could work to solve it. – individ Aug 26 '16 at 16:32
• What degrees? How many variables? You should consider posting a minimal example for a small $n$ at least. – dxiv Aug 26 '16 at 16:47
• In general, deciding whether a solution to a system of polynomial equations exists is undecidable (this was Hilbert's tenth problem). This implies that any algorithm giving a solution to a system of equations guaranteed to have a root has time complexity growing faster than any computable function. I'd suspect the situation is no better for systems with a unique solution, but I can't tell. – Milo Brandt Aug 26 '16 at 21:12

If (and it's a big if) some variable has degree $1$ in one of the polynomials, you can substitute for that variable, obtaining a smaller set of equations in fewer variables, but in general with higher degrees.
Somewhat more generally, you can take a Groebner basis of your system of polynomials. In the example, I took a grevlex basis (using Maple's Groebner package) and obtained $[1]$, indicating that the system is inconsistent. With $20$ replaced by $21$ in the fifth equation, the grevlex basis was $[-1+x_5, x_4-7, -2+x_3, -2+x_2, x_1-5]$, from which you can read off the solution.