Activation of naive CD8 T-cells can lead to the generation of multiple effector and memory subsets. Multiple parameters associated with activation conditions are involved in generating this diversity that is associated with heterogeneous molecular contents of activated cells. Although naive cell polarisation upon antigenic stimulation and the resulting asymmetric division are known to be a major source of heterogeneity and cell fate regulation, the consequences of stochastic uneven partitioning of molecular content upon subsequent divisions remain unclear yet. Here we aim at studying the impact of uneven partitioning on molecular-content heterogeneity and then on the immune response dynamics at the cellular level. To do so, we introduce a multiscale mathematical model of the CD8 T-cell immune response in the lymph node. In the model, cells are described as agents evolving and interacting in a 2D environment while a set of differential equations, embedded in each cell, models the regulation of intra and extracellular proteins involved in cell differentiation. Based on the analysis of in silico data at the single cell level, we 1 show that immune response dynamics can be explained by the molecular-content heterogeneity generated by uneven partitioning at cell division. In particular, uneven partitioning acts as a regulator of cell differentiation and induces the emergence of two coexisting sub-populations of cells exhibiting antagonistic fates. We show that the degree of unevenness of molecular partitioning, along all cell divisions, affects the outcome of the immune response and can promote the generation of memory cells.

ABSTRACT We employed our previously-described single-cell gene expression analysis CALISTA (Clustering And Lineage Inference in Single-Cell Transcriptional Analysis) to evaluate transcriptional uncertainty at the single-cell level using a stochastic mechanistic model of gene expression. We reconstructed a transcriptional uncertainty landscape during cell differentiation by visualizing single-cell transcriptional uncertainty surface over a two dimensional representation of the single-cell gene expression data. The reconstruction of transcriptional uncertainty landscapes for ten publicly available single-cell gene expression datasets from cell differentiation processes with linear, single or multi-branching cell lineage, reveals universal features in the cell differentiation trajectory that include: (i) a peak in single-cell uncertainty during transition states, and in systems with bifurcating differentiation trajectories, each branching point represents a state of high transcriptional uncertainty; (ii) a positive correlation of transcriptional uncertainty with transcriptional burst size and frequency; (iii) an increase in RNA velocity preceeding the increase in the cell transcriptional uncertainty. Finally, we provided biological interpretations of the universal rise-then-fall profile of the transcriptional uncertainty landscape, including a link with the Waddington’s epigenetic landscape, that is generalizable to every cell differentiation system.

Abstract Differentiation is the process whereby a cell acquires a specific phenotype, by differential gene expression as a function of time. This is thought to result from the dynamical functioning of an underlying Gene Regulatory Network (GRN). The precise path from the stochastic GRN behavior to the resulting cell state is still an open question. In this work we propose to reduce a stochastic model of gene expression, where a cell is represented by a vector in a continuous space of gene expression, to a discrete coarse-grained model on a limited number of cell types. We develop analytical results and numerical tools to perform this reduction for a specific model characterizing the evolution of a cell by a system of piecewise deterministic Markov processes (PDMP). Solving a spectral problem, we find the explicit variational form of the rate function associated to a large deviations principle, for any number of genes. The resulting Lagrangian dynamics allows us to define a deterministic limit of which the basins of attraction can be identified to cellular types. In this context the quasipotential, describing the transitions between these basins in the weak noise limit, can be defined as the unique solution of an Hamilton-Jacobi equation under a particular constraint. We develop a numerical method for approximating the coarse-grained model parameters, and show its accuracy for a symmetric toggle-switch network. We deduce from the reduced model an approximation of the stationary distribution of the PDMP system, which appears as a Beta mixture. Altogether those results establish a rigorous frame for connecting GRN behavior to the resulting cellular behavior, including the calculation of the probability of jumps between cell types.

Abstract Cell differentiation requires the integration of two opposite processes, a stabilizing cellular memory, especially at the transcriptional scale, and a burst of gene expression variability which follows the differentiation induction. Therefore, the actual capacity of a cell to undergo phenotypic change during a differentiation process relies upon a modification in this balance which favors change-inducing gene expression variability. However, there are no experimental data providing insight on how fast the transcriptomes of identical cells would diverge on the scale of the very first two cell divisions during the differentiation process. In order to quantitatively address this question, we developed different experimental methods to recover the transcriptomes of related cells, after one and two divisions, while preserving the information about their lineage at the scale of a single cell division. We analyzed the transcriptomes of related cells from two differentiation biological systems (human CD34+ cells and T2EC chicken primary erythrocytic progenitors) using two different single-cell transcriptomics technologies (sc-RT-qPCR and scRNA-seq). We identified that the gene transcription profiles of differentiating sister-cells are more similar to each-other than to those of non related cells of the same type, sharing the same environment and undergoing similar biological processes. More importantly, we observed greater discrepancies between differentiating sister-cells than between self-renewing sister-cells. Furthermore, a continuous increase in this divergence from first generation to second generation was observed when comparing differentiating cousin-cells to self renewing cousin-cells. Our results are in favor of a continuous and gradual erasure of transcriptional memory during the differentiation process.

Abstract The rise of single-cell data highlights the need for a nondeterministic view of gene expression, while offering new opportunities regarding gene regulatory network inference. We recently introduced two strategies that specifically exploit time-course data, where single-cell profiling is performed after a stimulus: HARISSA, a mechanistic network model with a highly efficient simulation procedure, and CARDAMOM, a scalable inference method seen as model calibration. Here, we combine the two approaches and show that the same model driven by transcriptional bursting can be used simultaneously as an inference tool, to reconstruct biologically relevant networks, and as a simulation tool, to generate realistic transcriptional profiles emerging from gene interactions. We verify that CARDAMOM quantitatively reconstructs causal links when the data is simulated from HARISSA, and demonstrate its performance on experimental data collected on in vitro differentiating mouse embryonic stem cells. Overall, this integrated strategy largely overcomes the limitations of disconnected inference and simulation. Author summary Gene regulatory network (GRN) inference is an old problem, to which single-cell data has recently offered new challenges and breakthrough potential. Many GRN inference methods based on single-cell transcriptomic data have been developed over the last few years, while GRN simulation tools have also been proposed for generating synthetic datasets with realistic features. However, except for benchmarking purposes, these two fields remain largely disconnected. In this work, building on a combination of two methods we recently described, we show that a particular GRN model can be used simultaneously as an inference tool, to reconstruct a biologically relevant network from time-course single-cell gene expression data, and as a simulation tool, to generate realistic transcriptional profiles in a non-trivial way through gene interactions. This integrated strategy demonstrates the benefits of using the same executable model for both simulation and inference.

Abstract According to Waddington’s epigenetic landscape concept, the differentiation process can be illustrated by a cell akin to a ball rolling down from the top of a hill (proliferation state) and crossing furrows before stopping in basins or “attractor states” to reach its stable differentiated state. However, it is now clear that some committed cells can retain a certain degree of plasticity and reacquire phenotypical characteristics of a more pluripotent cell state. In line with this dynamic model, we have previously shown that differentiating cells (chicken erythrocytic progenitors (T2EC)) retain for 24 hours the ability to self-renew when transferred back in self-renewal conditions. Despite those intriguing and promising results, the underlying molecular state of those “reverting” cells remains unexplored. The aim of the present study was therefore to molecularly characterize the T2EC reversion process by combining advanced statistical tools to make the most of single cell transcriptomic data. For this purpose, T2EC, initially maintained in a self-renewal medium (0H), were induced to differentiate for 24h (24H differentiating cells); then a part of these cells was transferred back to the self-renewal medium (48H reverting cells) and the other part was maintained in the differentiation medium for another 24h (48H differentiating cells). For each time point, cell transcriptomes were generated using scRT-qPCR and scRNAseq. Our results showed a strong overlap between 0H and 48H reverting cells when applying dimensional reduction. Moreover, the statistical comparison of cell distributions and differential expression analysis indicated no significant differences between these two cell groups. Interestingly, gene pattern distributions highlighted that, while 48H reverting cells have gene expression pattern more similar to 0H cells, they retained traces of their engagement in the differentiation process. Finally, Sparse PLS analysis showed that only the expression of 3 genes discriminates 48H reverting and 0H cells. Altogether, we show that reverting cells return to an earlier molecular state almost identical to undifferentiated cells and demonstrate a previously undocumented physiological and molecular plasticity during the differentiation process, which most likely results from the dynamic behavior of the underlying molecular network.

The in vivo erythropoiesis, which is the generation of mature red blood cells in the bone marrow of whole organisms, has been described by a variety of mathematical models in the past decades. However, the in vitro erythropoiesis, which produces red blood cells in cultures, has received much less attention from the modelling community. In this paper, we propose the first mathematical model of in vitro erythropoiesis. We start by formulating different models and select the best one at fitting experimental data of in vitro erythropoietic differentiation. It is based on a set of linear ODE, describing 3 hypothetical populations of cells at different stages of differentiation. We then compute confidence intervals for all of its parameters estimates, and conclude that our model is fully identifiable. Finally, we use this model to compute the effect of a chemical drug called Rapamycin, which affects all states of differentiation in the culture, and relate these effects to specific parameter variations. We provide the first model for the kinetics of in vitro cellular differentiation which is proven to be identifiable. It will serve as a basis for a model which will better account for the variability which is inherent to the experimental protocol used for the model calibration.

Inference of gene regulatory networks from gene expression data has been a long-standing and notoriously difficult task in systems biology. Recently, single-cell transcriptomic data have been massively used for gene regulatory network inference, with both successes and limitations. In the present work we propose an iterative algorithm called WASABI, dedicated to inferring a causal dynamical network from time-stamped single-cell data, which tackles some of the limitations associated with current approaches. We first introduce the concept of waves, which posits that the information provided by an external stimulus will affect genes one-by-one through a cascade, like waves spreading through a network. This concept allows us to infer the network one gene at a time, after genes have been ordered regarding their time of regulation. We then demonstrate the ability of WASABI to correctly infer small networks, which have been simulated in-silico using a mechanistic model consisting of coupled piecewise-deterministic Markov processes for the proper description of gene expression at the single-cell level. We finally apply WASABI on in-vitro generated data on an avian model of erythroid differentiation. The structure of the resulting gene regulatory network sheds a fascinating new light on the molecular mechanisms controlling this process. In particular, we find no evidence for hub genes and a much more distributed network structure than expected. Interestingly, we find that a majority of genes are under the direct control of the differentiation-inducing stimulus. In conclusion, WASABI is a versatile algorithm which should help biologists to fully exploit the power of time-stamped single-cell data.

Recent rise of single-cell studies revealed the importance of understanding the role of cell-to-cell variability, especially at the transcriptomic level. One of the numerous sources of cell-to-cell variation in gene expression is the heterogeneity in cell proliferation state. How cell cycle and cell size influences gene expression variability at single-cell level is not yet clearly understood. To deconvolute such influences, most of the single-cell studies used dedicated methods that could include some bias. Here, we provide a universal and automatic toxic-free label method, compatible with single-cell high-throughput RT-qPCR. This led to an unbiased gene expression analysis and could be also used for improving single-cell tracking and imaging when combined with cell isolation. As an application for this technique, we showed that cell-to-cell variability in chicken erythroid progenitors was negligibly influenced by cell size nor cell cycle.

It has been suggested that a switch from glycolysis, with lactate production, toward mitochondrial oxidative phosphorylation (OXPHOS) could be a driving force during stem cell differentiation. Based upon initial results, from our previous work, showing a drop in LDHA mRNA level during the differentiation of chicken erythroid progenitors, we studied metabolism behavior to question whether such switch might also be operating in those cells. We first analyzed the level of 9 enzymes, including LDHA, involved either in glycolysis or OXPHOS, in self-renewing and differentiating cells. Our results suggest that erythroid differentiation might be accompanied by an enhancement of the respiratory chains and glycolysis activities at 12h, followed by a strong decline of the glycolytic pathway and a stabilization of OXPHOS. To confirm that OXPHOS might be increased and glycolysis decreased during erythroid differentiation, we measured lactate concentration and mitochondrial membrane potential (MMP) of self-renewing and differentiating cells. Our findings show that at 12h-24h of differentiation, a surge of energy is needed, which could be fueled jointly by glycolysis and OXPHOS. Then the energy demand comes back to normal and might be supplied by OXPHOS instead of lactate production through glycolysis. These results support the hypothesis that erythroid differentiation is associated with a metabolic switch from glycolysis to OXPHOS. We also assessed LDHA role in erythroid progenitors self-renewal and the metabolic status changes. Inhibition experiments showed that LDHA activity could be involved in the maintenance of erythroid progenitors self-renewal, and its decline could influence their metabolic status. Finally, we investigated whether these metabolic rearrangements were necessary for erythroid differentiation. The addition of an inhibitor of the respiratory chains affected progenitors ability to differentiate, suggesting that the metabolic switch from glycolysis toward OXPHOS might act as a driving force for erythroid differentiation.