Skip to main content

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.

  • Article
  • Published:

Inferring population dynamics from single-cell RNA-sequencing time series data

Nature Biotechnology volume 37, pages 461–468 (2019)Cite this article

Subjects

Abstract

Recent single-cell RNA-sequencing studies have suggested that cells follow continuous transcriptomic trajectories in an asynchronous fashion during development. However, observations of cell flux along trajectories are confounded with population size effects in snapshot experiments and are therefore hard to interpret. In particular, changes in proliferation and death rates can be mistaken for cell flux. Here we present pseudodynamics, a mathematical framework that reconciles population dynamics with the concepts underlying developmental trajectories inferred from time-series single-cell data. Pseudodynamics models population distribution shifts across trajectories to quantify selection pressure, population expansion, and developmental potentials. Applying this model to time-resolved single-cell RNA-sequencing of T-cell and pancreatic beta cell maturation, we characterize proliferation and apoptosis rates and identify key developmental checkpoints, data inaccessible to existing approaches.

This is a preview of subscription content, access via your institution

Access options

Buy this article

  • Purchase on Springer Link
  • Instant access to full article PDF
Buy now

Prices may be subject to local taxes which are calculated during checkout

Fig. 1: A population-based view of single-cell RNA-seq time-series experiments: concept of pseudodynamics and example fits on a mouse embryonic stem cell differentiation data set.
Fig. 2: A trajectory model for T-cell maturation yields a description with higher resolution than a discretized description of the cell state space.
Fig. 3: Pseudodynamics density and parameter fits extend the stationary description of T-cell maturation.
Fig. 4: Pseudodynamics annotates a trajectory model with the position of a developmental checkpoint based on knock-out data, with the developmental potential and with developmental phases.
Fig. 5: Cell state space discretization and time-dependence of models of pancreatic beta cell maturation and proliferation.
Fig. 6: Likelihood-based model selection favors a state-dependent birth-death model over a state- and time-dependent model for beta cell maturation.

Similar content being viewed by others

Data availability

The Rag2 knockout Drop-seq data have been deposited in GEO under accession code GSE126579. The remaining thymic hematopoietic cell data are available from GEO under accession code GSE107910.

Code availability

The pseudodynamics model code and the presented examples are available through GitHub (https://github.com/theislab/pseudodynamics) and in the Supplementary Code.

References

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

    Article  CAS  Google Scholar 

  2. Taniguchi, K., Kajiyama, T. & Kambara, H. Quantitative analysis of gene expression in a single cell by qPCR. Nat. Methods 6, 503–506 (2009).

    Article  CAS  Google Scholar 

  3. Bandura, D. R. et al. Mass cytometry: technique for real time single cell multitarget immunoassay based on inductively coupled plasma time-of-flight mass spectrometry. Anal. Chem. 81, 6813–6822 (2009).

    Article  CAS  Google Scholar 

  4. Haghverdi, L., Büttner, M., Wolf, F. A., Buettner, F. & Theis, F. J. Diffusion pseudotime robustly reconstructs lineage branching. Nat. Methods 13, 845–848 (2016).

    Article  CAS  Google Scholar 

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

    Article  CAS  Google Scholar 

  6. Qiu, X. et al. Reversed graph embedding resolves complex single-cell trajectories. Nat. Methods 14, 979–982 (2017).

    Article  CAS  Google Scholar 

  7. Wolf, F. A. et al. PAGA: Graph abstraction reconciles clustering with trajectory inference through a topology preserving map of single cells. Genome Biol. 20, 59 (2019).

    Article  Google Scholar 

  8. Saelens, W., Cannoodt, R., Todorov, H. & Saeys, Y. A comparison of single-cell trajectory inference methods: towards more accurate and robust tools. Preprint at bioRxiv https://doi.org/10.1101/276907 (2018).

  9. Kafri, R. et al. Dynamics extracted from fixed cells reveal feedback linking cell growth to cell cycle. Nature 494, 480–483 (2013).

    Article  CAS  Google Scholar 

  10. Kuritz, K., Stöhr, D., Pollak, N. & Allgöwer, F. On the relationship between cell cycle analysis with ergodic principles and age-structured cell population models. J. Theor. Biol. 414, 91–102 (2017).

    Article  CAS  Google Scholar 

  11. Weinreb, C., Wolock, S., Tusi, B. K., Socolovsky, M. & Klein, A. M. Fundamental limits on dynamic inference from single-cell snapshots. Proc. Natl Acad. Sci. USA 115, E2467–E2476 (2018).

    Article  CAS  Google Scholar 

  12. Schiebinger, G. et al. Optimal-transport analysis of single-cell gene expression identifies developmental trajectories in teprogramming. Cell 176, 928–943.e22 (2019).

    Article  CAS  Google Scholar 

  13. Hashimoto, T., Gifford, D. & Jaakkola, T. Learning population-level diffusions with generative RNNs. in International Conference on Machine Learning 48, 2417–2426 (2016).

  14. Buchholz, V. R. et al. Disparate individual fates compose robust CD8+ T cell immunity. Science 340, 630–635 (2013).

    Article  CAS  Google Scholar 

  15. Cho, H. et al. Modelling acute myeloid leukaemia in a continuum of differentiation states. Lett. Biomath. 5, S69–S98 (2018).

    Article  CAS  Google Scholar 

  16. Segerstolpe, Å. et al. Single-cell transcriptome profiling of human pancreatic islets in health and type 2 diabetes. Cell Metab. 24, 593–607 (2016).

    Article  CAS  Google Scholar 

  17. Haber, A. L. et al. A single-cell survey of the small intestinal epithelium. Nature 551, 333–339 (2017).

    Article  CAS  Google Scholar 

  18. Yui, M. A. & Rothenberg, E. V. Developmental gene networks: a triathlon on the course to T cell identity. Nat. Rev. Immunol. 14, 529–545 (2014).

    Article  CAS  Google Scholar 

  19. Kernfeld, E. M. et al. A single-cell transcriptomic atlas of thymus organogenesis resolves cell types and developmental maturation. Immunity https://doi.org/10.1016/j.immuni.2018.04.015 (2018).

  20. Haghverdi, L., Buettner, F. & Theis, F. J. Diffusion maps for high-dimensional single-cell analysis of differentiation data. Bioinformatics 31, 2989–2998 (2015).

    Article  CAS  Google Scholar 

  21. Vosshenrich, C. A. J. et al. A thymic pathway of mouse natural killer cell development characterized by expression of GATA-3 and CD127. Nat. Immunol. 7, 1217–1224 (2006).

    Article  CAS  Google Scholar 

  22. Ribeiro, V. S. G. et al. Cutting edge: Thymic NK cells develop independently from T cell precursors. J. Immunol. 185, 4993–4997 (2010).

    Article  CAS  Google Scholar 

  23. Tang, Y. et al. Emergence of NK-cell progenitors and functionally competent NK-cell lineage subsets in the early mouse embryo. Blood 120, 63–75 (2012).

    Article  CAS  Google Scholar 

  24. Cook, A. M. Proliferation and lineage potential in fetal thymic epithelial progenitor cells. Edinburgh Research Archive, 1–195 (2010).

  25. Germain, R. N. T-cell development and the CD4-CD8 lineage decision. Nat. Rev. Immunol. 2, 309–322 (2002).

    Article  CAS  Google Scholar 

  26. Qiu, W.-L. et al. Deciphering Pancreatic Islet β Cell and α Cell Maturation Pathways and Characteristic Features at the Single-Cell Level. Cell Metab. 25, 1194–1205.e4 (2017).

    Article  CAS  Google Scholar 

  27. Bader, E. et al. Identification of proliferative and mature β-cells in the islets of Langerhans. Nature 535, 430–434 (2016).

    Article  CAS  Google Scholar 

  28. Herbach, N., Bergmayr, M., Göke, B., Wolf, E. & Wanke, R. Postnatal development of numbers and mean sizes of pancreatic islets and beta-cells in healthy mice and GIPRdn transgenic diabetic mice. PLoS One 6, e22814 (2011).

    Article  CAS  Google Scholar 

  29. Hija, A. et al. G0-G1 Transition and the restriction point in pancreatic β-cells in vivo. Diabetes 63, 578 (2014).

    Article  CAS  Google Scholar 

  30. Scaglia, L., Cahill, C. J., Finegood, D. T. & Bonner-Weir, S. Apoptosis participates in the remodeling of the endocrine pancreas in the neonatal rat. Endocrinology 138, 1736–1741 (1997).

    Article  CAS  Google Scholar 

  31. Waddington, C. H. Organisers and Genes (Cambridge Univ. Press, 1940).

  32. McKenna, A. et al. Whole-organism lineage tracing by combinatorial and cumulative genome editing. Science 353, aaf7907 (2016).

    Article  Google Scholar 

  33. Spanjaard, B. et al. Simultaneous lineage tracing and cell-type identification using CRISPR–Cas9-induced genetic scars. Nat. Biotechnol. 36, 469 (2018).

    Article  CAS  Google Scholar 

  34. La Manno, G. et al. RNA velocity of single cells. Nature 560, 494–498 (2018).

    Article  Google Scholar 

  35. Cohen, S. D., Hindmarsh, A. C. & Dubois, P. F. CVODE, a stiff/nonstiff ODE solver in C. Computers in Physics 10/2, 138–143 (1996).

  36. Fröhlich, F., Theis, F. J., Rädler, J. O. & Hasenauer, J. Parameter estimation for dynamical systems with discrete events and logical operations. Bioinformatics 33, 1049–1056 (2017).

    PubMed  Google Scholar 

  37. Stapor, P. et al. PESTO: Parameter EStimation TOolbox. Bioinformatics 34, 705–707 (2018).

    Article  CAS  Google Scholar 

  38. Raue, A. et al. Lessons learned from quantitative dynamical modeling in systems biology. PLoS One 8, e74335 (2013).

    Article  CAS  Google Scholar 

  39. Stubbington, M. J. T. et al. T cell fate and clonality inference from single-cell transcriptomes. Nat. Methods 13, 329–332 (2016).

    Article  Google Scholar 

  40. Satija, R., Farrell, J. A., Gennert, D., Schier, A. F. & Regev, A. Spatial reconstruction of single-cell gene expression data. Nat. Biotechnol. 33, 495–502 (2015).

    Article  CAS  Google Scholar 

  41. Macosko, E. Z. et al. Highly parallel genome-wide expression profiling of individual cells using nanoliter droplets. Cell 161, 1202–1214 (2015).

    Article  CAS  Google Scholar 

  42. Waltman, L. & van Eck, N. J. A smart local moving algorithm for large-scale modularity-based community detection. Eur. Phys. J. B 86, 471 (2013).

    Article  Google Scholar 

  43. Wolf, F. A., Angerer, P. & Theis, F. J. SCANPY: large-scale single-cell gene expression data analysis. Genome Biol. 19, 15 (2018).

    Article  Google Scholar 

Download references

Acknowledgements

We would like to express our gratitude towards Ping Xu and Kashfia Neherin for mouse husbandry and tissue preparation. F.J.T. acknowledges financial support by the Graduate School QBM, the German Research Foundation (DFG) within the Collaborative Research Centre 1243, Subproject A17, by the Helmholtz Association (Incubator grant sparse2big, ZT-I-0007) and by the Chan Zuckerberg Initiative DAF (advised fund of Silicon Valley Community Foundation, 182835). R.M. acknowledges financial support by the Leona M. and Harry B. Helmsley Charitable Trust (2015PG-T1D035), a Charles H. Hood Foundation Child Health Research Award, the Glass Charitable Foundation and the National Institutes of Health (1DP3DK111898, R01 AI132963, UC4 DK104218). J.H. acknowledges financial support by the German Research Foundation (DFG) (HA 7376/1–1) and the German Federal Ministry of Education and Research (BMBF) within the SYS-Stomach project (01ZX1310B). H.L. acknowledges financial support by the Helmholtz Society and German Center for Diabetes Research (DZD e.V.). D.S.F. acknowledges financial support by a German research foundation (DFG) fellowship through the Graduate School of Quantitative Biosciences Munich (QBM) (GSC 1006) and by the Joachim Herz Stiftung.

Author information

Author notes

  1. These authors contributed equally: David S. Fischer and Anna K. Fiedler.

Authors and Affiliations

  1. Institute of Computational Biology, Helmholz Zentrum München, Neuherberg, Germany

    David S. Fischer, Anna K. Fiedler, Jan Hasenauer & Fabian J. Theis

  2. TUM School of Life Sciences Weihenstephan, Technical University of Munich, Freising, Germany

    David S. Fischer & Fabian J. Theis

  3. Department of Mathematics, Technical University of Munich, Garching bei München, Germany

    Anna K. Fiedler, Jan Hasenauer & Fabian J. Theis

  4. Program in Molecular Medicine, Diabetes Center of Excellence, University of Massachusetts Medical School, Worcester, MA, USA

    Eric M. Kernfeld, Ryan M. J. Genga & Rene Maehr

  5. Institute of Diabetes and Regeneration Research, Helmholtz Zentrum München, Neuherberg, Germany

    Aimée Bastidas-Ponce, Mostafa Bakhti & Heiko Lickert

  6. Institute of Stem Cell Research, Helmholtz Zentrum München, Neuherberg, Germany

    Aimée Bastidas-Ponce, Mostafa Bakhti & Heiko Lickert

  7. Medical Faculty, Technical University of Munich, Munich, Germany

    Aimée Bastidas-Ponce & Heiko Lickert

  8. German Center for Diabetes Research (DZD), Neuherberg, Germany

    Aimée Bastidas-Ponce, Mostafa Bakhti & Heiko Lickert

Authors
  1. David S. Fischer
    View author publications

    You can also search for this author in PubMed Google Scholar

  2. Anna K. Fiedler
    View author publications

    You can also search for this author in PubMed Google Scholar

  3. Eric M. Kernfeld
    View author publications

    You can also search for this author in PubMed Google Scholar

  4. Ryan M. J. Genga
    View author publications

    You can also search for this author in PubMed Google Scholar

  5. Aimée Bastidas-Ponce
    View author publications

    You can also search for this author in PubMed Google Scholar

  6. Mostafa Bakhti
    View author publications

    You can also search for this author in PubMed Google Scholar

  7. Heiko Lickert
    View author publications

    You can also search for this author in PubMed Google Scholar

  8. Jan Hasenauer
    View author publications

    You can also search for this author in PubMed Google Scholar

  9. Rene Maehr
    View author publications

    You can also search for this author in PubMed Google Scholar

  10. Fabian J. Theis
    View author publications

    You can also search for this author in PubMed Google Scholar

Contributions

F.J.T. conceived the study; F.J.T. and J.H. conceived the model and supervised analyses; R.M., H.L. and D.S.F. conceived the experiments; A.K.F. implemented the models and performed the parameter estimation in all examples; D.S.F. and E.M.K. performed the initial computational analysis of the thymus data; R.M.J.G. performed the experiments for the thymus study; D.S.F. performed the initial computational analysis of the pancreas data; A.B.P. and M.B. performed the experiments for the pancreas study; and D.S.F., A.K.F, J.H. and F.J.T. wrote the manuscript with assistance from all authors.

Corresponding author

Correspondence to Fabian J. Theis.

Ethics declarations

Competing interests

The authors declare no competing interests.

Additional information

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

Integrated supplementary information

Supplementary Figure 1 Overview of all thymic hematopoietic cells (n = 10,895) observed with single-cell RNA-seq.

(a-f) Time (a-c) or louvain group (d-f) label superimposed on cells in t-SNE based on data processed for louvain clustering (a,d), in in t-SNE based on data processed for diffusion pseudotime analysis used throughout the paper (b,e) and in diffusion map based on data before filtering of putative myeloid or dendritic cells (c,f). See also Supplementary data 1.2 for the further details on this analysis. (g) Cell type labels for louvain groups from marker genes. stem: hematopoietic stem cells, DN2a/b (double-negative T-cells), phase 2 (Nat. Rev. Immunol. 14, Yui, M. A. & Rothenberg, E. V., 2014), phase 3 (Nat. Rev. Immunol. 14, Yui, M. A. & Rothenberg, E. V., 2014) and DP (double-positive T-cells) are T-cell maturation stages, NCL_ILC: non-conventional lymphoid cells (Immunity, Kernfeld, E. M. et al., 2018) with innate lymphoid cell character, NCL_ydTC: non-conventional lymphoid cells with γδ-T-cell cell character, myeloid: putative myeloid and dendritic cells. (h) T-cell development phase across merged louvain groups from (g). (i) Abstracted graph from paga (bioRxiv, Wolf F. A. et al., 2017) at a connectivity threshold of 0.02 with cluster identity superimposed. See also Supplementary data 1.2 for the further details on this analysis.

Supplementary Figure 2 Myeloid and dendritic cell subclusters of the thymic hematopoietic cells observed with single-cell RNA-seq.

Here, we sub-clustered cluster 14_myeloid from Supplementary Fig. 2a into three clusters, again with louvain clustering. We show the expression of a dendritic and a myeloid cell marker across the t-SNE based on all thymic hematopoietic cells (n = 10895) observed with single-cell RNA-seq (a, b) to to show the overall restriction of these markers to the cluster 14_myeloid and across the three sub-clusters (c, d) to show the heterogeneity in this cluster. The violin plots in (c, d) are based on kernel density estimates and on n = 130 cells in group 0, n = 118 in group 1 and n = 102 in group 2. (a, c) Dendritic cell marker Cd74 log expression superimposed on t-SNE (a) and by louvain sub-cluster (c). (b, d) Myeloid cell marker Lyz2 log expression superimposed on t-SNE (b) and by louvain sub-cluster (d). See also Supplementary data 1.2 for the further details on this analysis.

Supplementary Figure 3 Myeloid, dendritic cell, leukocyte and T-cell marker gene expression by Louvain group.

The cluster labels are explained in Supplementary Fig. 1. The violin plots are based on kernel density estimates of log gene expression and are based on the following number of cells per louvain group: 0_stem: n = 1314, 1_DN2a_1: n = 906, 2_DN2a_2: n = 588, 3_Phase2_1: n = 568, 4_Phase2_2: n = 625, 5_Phase2_3: n = 687, 6_Phase2_4: n = 1050, 7_Phase3_1: n = 1234, 8_Phase3_2: n = 784, 9_DP_1: n = 463, 10_DP_2: n = 759, 11_DP_3: n = 881, 12_NCL_ydTC: n = 402, 14_myeloid: n = 350,13_NCL_ILC: n = 284. (a-f) Myeloid and dendritic cell marker genes: Itgax (a), Itgam (b), Adgre1 (c), Lyz2 (d), Cd74 (e), Anpep (f). (g-j) Leukocyte and T-cell marker genes: Ptprc (g), Cd34 (h), Notch1 (i), Bcl11b (j).

Supplementary Figure 4 CD3 subunit and T-cell development marker gene expression by Louvain group.

The cluster labels are explained in Supplementary Fig. 1. The violin plots are based on kernel density estimates and are based on the following number of cells per louvain group: 0_stem: n = 1314, 1_DN2a_1: n = 906, 2_DN2a_2: n = 588, 3_Phase2_1: n = 568, 4_Phase2_2: n = 625, 5_Phase2_3: n = 687, 6_Phase2_4: n = 1050, 7_Phase3_1: n = 1234, 8_Phase3_2: n = 784, 9_DP_1: n = 463, 10_DP_2: n = 759, 11_DP_3: n = 881, 12_NCL_ydTC: n = 402, 14_myeloid: n = 350,13_NCL_ILC: n = 284. (a-c) Cd3 subunits: Cd3d (a), Cd3e(b), Cd3g (c). (d-i) T-cell development marker genes: Flt3 (d), Cd44 (e), Il2ra (f), Ptcra (g), Cd8a (h), Cd8b1 (i).

Supplementary Figure 5 αβ T-cell, γδ T-cell and natural killer cell marker gene expression by Louvain group.

The cluster labels are explained in Supplementary Fig. 1. The violin plots are based on kernel density estimates and are based on the following number of cells per louvain group: 0_stem: n = 1314, 1_DN2a_1: n = 906, 2_DN2a_2: n = 588, 3_Phase2_1: n = 568, 4_Phase2_2: n = 625, 5_Phase2_3: n = 687, 6_Phase2_4: n = 1050, 7_Phase3_1: n = 1234, 8_Phase3_2: n = 784, 9_DP_1: n = 463, 10_DP_2: n = 759, 11_DP_3: n = 881, 12_NCL_ydTC: n = 402, 14_myeloid: n = 350,13_NCL_ILC: n = 284. (a,b) αβ-T-cells marker genes: Tcra (TCRα) (a), Tcrb (TCRβ) (b). (c-k) γδ-T-cells and natural killer cell marker genes: Tcrg (TCRγ) (c), Tcrd (TCRδ) (d), Ifng (e), Il17a (f), Il2rb (g), Ncr1 (h), Klrd1 (i), Klrb1b (j), Klrb1f (k).

Supplementary Figure 6 Cell type partitioning in the diffusion map of the thymic hematopoietic cells observed with single-cell RNA-seq.

DC: diffusion component. (a) Diffusion map on all thymic hematopoietic cells (n = 10895). Colour: branch allocation according to diffusion pseudotime. (b,c,d) Diffusion map on gated thymic hematopoietic cells: T-cells and non-conventional lymphocytes (n = 10872). (b) Branch allocation according to diffusion pseudotime. (c) Adapted diffusion pseudotime-based branch allocation for pseudodynamics based on a linear cut in the diffusion component 1 versus 2 plane. (d) Branching and non-branching region allocation for pseudodynamics: Cells are in the branching region if they have a diffusion pseudotime coordinate smaller than 0.25 and are not early thymic progenitors.

Supplementary Figure 7 Surface marker and transcription factor expression across cell state in thymic hematopoietic cells observed with single-cell RNA-seq.

(a-c) Surface marker expression across cell state. (a) Natural cubic spline fits to log expression data with cell state covariate (df=10). Shown are markers that are traditionally used to distinguish αβ-T-cell stages. Vertical lines indicate the right boundaries of cell stages in cell state. (b,c) Heatmap of all surface marker genes that are differentially expressed in cell state (log counts with limma (Nucleic Acid Res. 43, Ritchie, M. E. et al., 2015), natural cubic spline model with df = 4, q-value threshold of 1e-5) either in αβ-T-cell (b, n = 10079) or the non-conventional lymphocyte lineage (c, n = 1070). Genes are ordered separately for both lineages by the cell state at which their expression is maximal. Lists of all genes in these two heatmaps are supplied in Supplementary data 2.2.1 (αβ-T-cell lineage) and 2.2.2 (non-conventional lymphocyte lineage). The gene annotation used is supplied in Supplementary data 2.3. The underlying differential expression analysis results are supplied Supplementary data 2.4. (d,e) Transcription factor expression across cell state. (d) Natural cubic spline fits to log expression data with cell state covariate (df = 10). (e) Heatmap of all transcription factor genes that are differentially expressed in cell state (log counts with limma (Nucleic Acid Res. 43, Ritchie, M. E. et al., 2015), natural cubic spline model with df = 4, q-value threshold of 1e-5) in the αβ-T-cell lineage. Genes are ordered by peak time. A list of all genes in this heatmap is supplied in Supplementary data 2.1. The gene annotation used is supplied in Supplementary data 2.3. The underlying differential expression analysis results are supplied Supplementary data 2.4.

Supplementary Figure 8 Mean gene expression of marker gene groups in αβ T-cell branch in cell state bins.

(a) Fraction Cd4/Cd8 positive cells in cell state bins. (b) Surface markers expression. (c) Bcl2-family expression. (d) Nfat-family expression.

Supplementary Figure 9 Pseudodynamics model fits to T-cell maturation data with diffusion pseudotime as cell state.

(a-b) Cross validation results (leave-one-time-point-out) of pseudodynamics fits on T-cell maturation with diffusion pseudotime as cell state. (a) Overall prediction error on withheld data by regularization parameter value omega. (b) Regularized log-likelihood value of prediction at held-out time point (prediction error) by regularization parameter value omega and time point. (c-i) Observed density, model fit to full data (simulation) and imputed density (simulation_cv) on T-cell lineage at a given time point. The cell state shown is the cell state used in the main text linearly scaled into the interval [0,0.9] and extended to 1. Accordingly, there are no observations in (0.9,1]. The imputed density is the model fit of a model trained on all remaining time points with a regularization parameter of 10 (leave-one-time-point-out cross validation). (c) E13.5, (d) E14.5, (e) E15.5, (f) E16.5, (g) E17.5, (h) E18.5, (i) P0. (j) Population size estimates. Observed (points) and simulated (line) total size of thymic hematopoietic compartment in a thymic lobe with 95% confidence interval on simulated data (shaded) and observed data with one standard deviation around the mean as error bars. The population size observations are based on 5 replicates for t = 12.5 to t = 17.5 and on two replicates for t = 18.5 and t = 19.5. Replicates are independent measurements based on separate thymus samples for t = 12.5 to 18.5 and are independent measurements based on the two lobes of a single thymus for t = 19.5.

Supplementary Figure 10 Pseudodynamics parameter fits across regularization hyperparameters for T-cell maturation data with diffusion pseudotime cell state.

(a-c) Maximum likelihood estimator spline fit of birth-death (a) and drift parameter (b) and diffusion parameter (c) by regularization parameter (rho).

Supplementary Figure 11 Monocle2 embedding of thymic hematopoietic cells observed with single-cell RNA-seq.

(a,b) Monocle2 pseudotime (a) and time (b) superimposed on monocle2 embedding based on n = 10705 cells. (c) Gene expression as function of pseudotime with spline interpolation (Nat. Methods, Qiu X. et al., 2017) of genes related to beta-selection: Rorc, Bcl2l1 (Bcl-xL) and Bcl2. (d-k) Cell counts in hexagonal cell state bins in monocle2 embedding by clusters of thymocytes (defined in Supplementary Fig. 2b) with n = 1314 cells in 0_stem (d), 1494 cells in 1_DN2a (e), n = 2930 cells in 2_Phase2 (f), n = 2018 cells in 3_Phase3 (g), n = 2103 cells in 4_DP (h), n = 117 cells in 5_NCL_ILC (i), n = 402 cells in 6_NCL_ydTC (j), n = 327 cells in 7_myeloid (k).

Supplementary Figure 12 Monocle2 pseudotime assignment (cell state) of thymic hematopoietic cells observed with single-cell RNA-seq and pseudodynamics model fits to this Monocle2 cell state.

(a,b) Distribution of cells by sample across cell state (monocle2 pseudotime) on T-cell lineage. Colour: time point in days after fertilisation. (a) Kernel density estimate of union of all samples per time point (n = 442 at t = 12.5, n = 1795 at t = 13.5, n = 1052 at t = 14.5, n = 1013 at t = 15.5, n = 2616 at t = 16.5, n = 1966 at t = 17.5, n = 882 at t = 18.5, n = 939 at t = 19.5). (b) Box plot of each sample per time point (n = {152, 58, 232} at t = 12.5; n = {628, 476, 691} at t = 13.5; n = {508, 544} at t = 14.5; n = {437, 576} at t = 15.5; n = {870, 918, 828} at t = 16.5; n = {975, 991} at t = 17.5; n = {455, 427} at t = 18.5; n = {422, 517} at t = 19.5). Here, pseudotime coordinates are computed based on all replicates. Replicates are independent Drop-seq samples which are based on separate thymus samples, the two replicates at P0 are based on the two lobes of a single thymus. The center of each boxplots is the sample median, the whiskers extend from the upper (lower) hinge to the largest (smallest) data point no further than 1.5 times the interquantile range from the upper (lower) hinge. (c) Scatter plot of diffusion pseudotime cell state versus moncle2 cell state (pseudotime) on T-cell lineage cells. (d-f) Pseudodynamics parameter fits to monocle2 cell state: (d) Birth-death, (e) drift and (f) diffusion parameter estimates for three different regularisation parameters (rho).

Supplementary Figure 13 Rag1 and Rag2 knockout cells fall into the wild-type developmental manifold of T-cell maturation.

(a-c) Diffusion map projections of the combined wild-type and knock-out data sets show that the knockout cells lie within the wild-type manifold in each projection. KO: knockout, WT: wild-type. (d) The overall structure of the manifold is conserved between a diffusion map computed just on wild-type and one computed on wild-type and mutant samples: Scatter plot of diffusion pseudotime coordinates of wild-type cells in wild-type only (WT) and in wild-type with mutants (WT+RagKO) data sets.

Supplementary Figure 14 Developmental arrest of T-cell maturation in Rag1 and Rag2 knockout mice at β-selection.

(a,b) Rag1 and Rag2 knockout mice have T-cell populations which are delayed in developmental progress along T-cell maturation measured with diffusion pseudotime compared to age-matched wild-type mice. The p-value is the result of the one-sided Kolmogorov-Smirnov (KS) test which was used to test whether the knock-out empirical cumulative density function of the knock-out lies below that of the wild type cells. (a) Rag2 knockout (Rag2KO) (n = 423 cells) and wild-type (n = 971 cells) samples at E14.5. (b) Rag1 knockout (Rag1KO) (n = 1565 cells) and wild-type (n = 2484 cells) samples at E16.5. (c,d) Mean marker gene expression by sample of age matched wild-type (WT) and Rag1/Rag2 knockout (Rag1KO, Rag2KO) embryos at E14.5 (c) and E16.5 (d). (e) Beta-selection point estimators by regularization parameter (reg). (f) Scatter plot of diffusion pseudotime coordinates obtained based on wild-type and on merged wild-type and mutant mouse samples. The dependency can be approximated with a smooth function (blue line: natural cubic spline fit of degree five). Yellow line: identity function.

Supplementary Figure 15 Pseudodynamics model fit to the continuous pancreatic beta cell maturation data set.

(a-f) Kernel density estimator of observed density (shaded) and best model fit (lines) of β-cell population density across cell state by time point for models with state- (s) or state- and time-dependent (st) birth-death rates. (g) Observed number of β-cells in black with one standard deviation around the mean as error bars and population size fits by model with the following number of replicates per time point: n = 3 at t = 0, n = 3 at t = 5.5, n = 3 at t = 10.5, n = 8 at t = 11.5, n = 3 at t = 15.5, n = 5 at t = 46.5. Replicates were separately measured in one unique animal per replicate.

Supplementary information

Supplementary Figures and Text

Supplementary Figures 1–15 and Supplementary Notes 1–3

Reporting Summary

Supplementary Video 1

Pseudodynamics fit to in vitro differentiation of embryonic stem cells

Supplementary Video 2

Pseudodynamics fit to T-cell maturation

Supplementary Video 3

Cell flux captured with pseudodynamics during T-cell maturation

Supplementary Data 1

Single-cell RNA-seqm analysis workflows

Supplementary Data 2

Gene-wise results tables

Supplementary Data 3

Pseudodynamics input

Supplementary Data 4

Pancreatic beta cell data

Supplementary Software

Pseudodynamics code used in this study.

Rights and permissions

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Fischer, D.S., Fiedler, A.K., Kernfeld, E.M. et al. Inferring population dynamics from single-cell RNA-sequencing time series data. Nat Biotechnol 37, 461–468 (2019). https://doi.org/10.1038/s41587-019-0088-0

Download citation

  • Received:

  • Accepted:

  • Published:

  • Issue Date:

  • DOI: https://doi.org/10.1038/s41587-019-0088-0

This article is cited by

Search

Advanced search

Quick links

Nature Briefing

Sign up for the Nature Briefing newsletter — what matters in science, free to your inbox daily.

Get the most important science stories of the day, free in your inbox. Sign up for Nature Briefing