Universality of clonal dynamics poses fundamental limits to identify stem cell self-renewal strategies.

How adult stem cells maintain self-renewing tissues is commonly assessed by analysing clonal data from in vivo cell lineage-tracing assays. To identify strategies of stem cell self-renewal requires that different models of stem cell fate choice predict sufficiently different clonal statistics. Here, we show that models of cell fate choice can, in homeostatic tissues, be categorized by exactly two 'universality classes', whereby models of the same class predict, under asymptotic conditions, the same clonal statistics. Those classes relate to generalizations of the canonical asymmetric vs. symmetric stem cell self-renewal strategies and are distinguished by a conservation law. This poses both challenges and opportunities to identify stem cell self-renewal strategies: while under asymptotic conditions, self-renewal models of the same universality class cannot be distinguished by clonal data only, models of different classes can be distinguished by simple means.

Introduction

Adult stem cells are the key players for maintaining and renewing biological tissue, due to their ability to persistently produce tissue cells through cell division and differentiation (National Institute of Health, 2009). For maintaining tissues in a homeostatic state, it is crucial that stem cells adopt suitable self-renewal strategies, a pattern of stem cell fate choices that balances proliferation and differentiation; otherwise, imbalanced proliferation may lead to hyperplasia and cancer. Therefore, the understanding and identification of stem cell self-renewal strategies has been a major goal of stem cell biology ever since the discovery of adult stem cells.

Classically, two stem cell self-renewal strategies have been proposed (Potten and Loeffler, 1990; Simons and Clevers, 2011a): following the Invariant Asymmetric division (IA) strategy, stem cells undertake only asymmetric divisions, whose outcome is one differentiating cell and one stem cell as daughter cells. The other proposed strategy, Population Asymmetry (PA) (Potten and Loeffler, 1990; Simons and Clevers, 2011a; Watt and Hogan, 2000; Klein and Simons, 2011), features additionally symmetric divisions, which produce either two stem cells or two differentiating cells as daughters, yet in balanced proportions. Both patterns of cell fate choice leave the number of cells on average unchanged and thus can maintain homeostasis. Assessing stem cell self-renewal strategies experimentally is difficult in vivo, since direct observation of cell divisions is rarely possible. Yet, through genetic cell lineage-tracing assays, the statistics of clones – the progeny of individual cells – can be obtained, and via mathematical modeling assessing cell fate dynamics became possible. With such an approach several studies suggested that population asymmetry prevails in many mouse tissues (e.g. Clayton et al., 2007; Lopez-Garcia et al., 2010; Simons and Clevers, 2011b; Doupé et al., 2012; Klein et al., 2010).

However, the interpretation of those studies has been challenged by a suggested alternative self-renewal strategy, called Dynamic Heterogeneity (DH), featuring some degree of cell fate plasticity (Greulich and Simons, 2016). In this model, all stem cell divisions are asymmetric, yet it is in agreement with the experimental clonal data that had previously been shown to agree also with the population asymmetry strategy. Thus, those two strategies are not distinguishable in view of the clonal data.

This raises the question to what extent different stem cell self-renewal strategies can be distinguished at all via clonal data (Klein and Simons, 2011; Greulich, 2019). Here, we address this question by studying models for stem cell fate choice, which define the self-renewal strategies, in their most generic form. We show that many cell fate models predict, under asymptotic conditions, the same clonal statistics and thus cannot be distinguished via clonal data from cell lineage-tracing experiments. In particular, we find that there exist two particular classes of stem cell self-renewal strategies: one class of models which all generate an Exponential distribution of clone sizes (the number of cells in a clone) after sufficiently large time, and one which generates a Normal distribution under sufficiently fast stem cell proliferation. Crucially, these two classes are not differentiated via the classical definitions of symmetric and asymmetric stem cell divisions, but by whether or not a subset of cells is conserved. These classes thus bear resemblance to 'universality classes' known from statistical physics, as suggested in Klein and Simons, 2011. This leads us to a more generic, and in this context more useful, definition of the terms ‘symmetric’ and ‘asymmetric’ divisions. Notably, however, we find that the conditions for the emergence of universality are not always fulfilled in real tissues, which provides chances, but also further challenges, for the identification of stem cell fate choices in homeostatic tissues.

Strategies for stem cell self-renewal

The two classical stem cell self-renewal strategies, Invariant Asymmetry (IA) and Population Asymmetry (PA) (Potten and Loeffler, 1990; Simons and Clevers, 2011a; Watt and Hogan, 2000; Klein and Simons, 2011), are commonly described in terms of two cell types: stem cells (S) which can self-renew (i.e. divide without reducing their potential to divide in the future); and differentiating cells (D). Both strategies can be expressed in terms of a single parametrized stochastic model, a multi-type branching process (Haccou et al., 2005), defined by the outcomes of cell divisions (the cell fate choices),where cells of type S divide with rate λ. Here, a daughter cell configuration corresponds to symmetric self-renewal division and to symmetric differentiation, while daughter cells of different type, , marks an asymmetric division. In the basic model version, a cell of type D is eventually lost with rate γ, (corresponding to death, shedding, or emigration of D-cells), while other versions may include the possibility of limited proliferation as committed progenitor cells. The two self-renewal strategies, IA and PA, are distinguished by the value of the symmetric division fraction r: the PA model corresponds to any ; the IA model is defined by , that is, only asymmetric divisions occur.

To maintain homeostasis, the number of cells must stay, on average, constant. Thus cells following the PA strategy must regulate the probabilities of symmetric self-renewal and differentiation to be exactly equal, whereas for the IA model this is trivially assured. However, only for the IA model is the number of stem cells strictly conserved, that is, no gain or loss of stem cells is possible.

A way to assess self-renewal strategies experimentally is via genetic cell-lineage tracing (Kretzschmar and Watt, 2012; Blanpain and Simons, 2013): By marking single cells with an inheritable genetic marker (through a Cre-Lox system [Soriano, 1999; Sauer, 1998]) each cell’s progeny, called a clone, which retain that marker, can be traced. The number of cells per clone, that is the clone size, is measured and the statistical frequency distribution of clone sizes (clone size distribution) determined. To test the cell fate choice models on that data, one evaluates the models with a single cell as initial condition and samples the outcome in terms of the final cell numbers – the size of a virtual clone. In the basic version of the model (i.e. when ), the IA and PA models predict, respectively, a Poisson and an Exponential clone size distribution for large times (Klein and Simons, 2011; Antal and Krapivsky, 2010) (see also the Appendix, 'Invariant Asymmetry and Population Asymmetry models'). Thus, they are fundamentally different and can easily be distinguished when compared with clonal data. By a series of lineage-tracing experiments it was confirmed that Exponential clone size distributions prevail for most mouse tissues, which thus exclude the IA model and support the PA strategy (Clayton et al., 2007; Lopez-Garcia et al., 2010; Simons and Clevers, 2011b; Doupé et al., 2012; Klein et al., 2010).

While this seemed to settle the case in favour of the PA strategy, at least for most adult mouse tissues, this was challenged by a third type of strategy, the DH model (Greulich and Simons, 2016). Motivated by the emerging view of prevailing cell plasticity (Blanpain and Fuchs, 2014; Tetteh et al., 2015; Tetteh et al., 2016; Donati and Watt, 2015), the DH model considers the possibility of reversible switching between two cell types:where symbols at arrows denote the process rates (frequency of events). This strategy is also capable of maintaining a homeostatic population if . Notably, the DH model only features asymmetric divisions (in that daughter cells are of different type), like the IA model, yet the DH model predicts clonal statistics that are indistinguishable from the PA model (Greulich and Simons, 2016). This means that in view of the existing clonal data for mouse tissues, the DH model, may as well describe the real cell fate dynamics. More fundamentally, this implies that the PA and DH model cannot be distinguished via plain clonal data, which poses fundamental limitations to the common approach to use lineage tracing for determining cell fate choices.

This demonstrates that the classical definition of asymmetric and symmetric divisions is not always suitable to distinguish cell fate strategies in view of clonal data alone. In general, cell fate dynamics may be much more complex than the simplified models described above, as there may exists a plethora of cell (sub-)types in a tissue. However, to what extent would it be possible to distinguish details of potentially rather complex cell fate dynamics models through comparison with clonal data at all? This is only the case if the clonal statistics are sufficiently different. In the following, we study cell fate models in their most generic form, and analyze what clonal statistics would be expected.

Results

Model generalization

Let us consider the dynamics of a generic system of cells, characterized by a number m of possible cell states X, . We define a cell state here as a group of cells showing common properties (e.g. any cell sub-type classification). Most generally, cells in a state X may be able to divide, producing daughter cells of any cell states X and X (where , that is, simple cell duplication, is possible). Furthermore, any cell state X may turn into another state X or may be lost (through emigration, shedding, or death). Hence, we can write a generic cell fate model as,where . In this model, is the rate of division of cells in state X and the parameter corresponds to the proportion of division outcomes producing daughter cells of state X and state X; ω is the transition rate from state X to state X and the loss rate from state X.

The dynamics of each cell in Equations 3-5 could depend on the cell environment through spatial, cell-extrinsic regulation of cell fate. However, the clonal statistics of spatial models that include cell-extrinsic regulation of cell fate (models of the voter type [Clifford and Sudbury, 1973]) are, in the long term, the same as for the corresponding branching process models (Haccou et al., 2005), as Equations 3-5 are, except for one-dimensional arrangements of cells (as shown in Klein and Simons, 2011; Bramson and Griffeath, 1980). Here, we are focussing on the long-term clonal statistics of self-renewal strategies, and since this is not affected by cell-extrinsic regulation, for tissues with two-dimensional or three-dimensional arrangements of dividing cells (like epithelial sheets, and volumnar tissue), we wish to keep the analysis simple and therefore choose dynamics (and thus the parameters ) to be independent of the cell environment.

In the following, we study the dynamics of cell numbers in each state X, . To gain initial insight into those dynamics, let us first consider the time evolution of the mean cell numbers, , given by,in which is the probability of having a daughter cell in state X produced upon division of a cell in state X. This linear system of differential equations can be written more compactly in terms of the mean cell number vector ,with being the matrixwhere we defined the total transition rate , combining all transitions from X to X by cell divisions and direct transitions, and the local loss rate .

Models of the form Equations 3–5 are not generally in homeostasis, which in this context is defined by the existence of a stationary state , with , that is (Lyapunov) stable and non-trivial (for a discussion, see the Appendix 'Conditions for homeostasis'). This can in principle be assessed through the spectral properties of A (Åström and Murray, 2008), but applying spectral conditions explicitly is unwieldy and difficult to interpret biologically. For a more intuitive view, we interpret the system, Equation 7, as a network (graph): the matrix A can be interpreted as the adjacency matrix of the cell state network. This is a weighted directed graph in which cell states correspond to the graph’s nodes and a link from state X to X exists where a transition is possible, that is, when . The value of also denotes the link weights (diagonal elements of A can be considered as self-links). Now, we note that Equation 7 is linear and cooperative, that is, the off-diagonal elements of matrix A are non-negative, and for such systems more simple and intuitive conditions for homeostasis exist (Greulich et al., 2019), based on a decomposition into the network’s Strongly Connected Component (SCC). An SCC is a sub-graph that groups nodes which are strongly connected, that is, which are mutually connected by paths (more accurately: two nodes, X and X are strongly connected if there exists a path from X to X and from X to X on the network). An example of such a decomposition, which yields an acyclic condensed network that contains SCCs as nodes and directed links between them, is shown in Figure 1.

Figure 1.

Illustration of the decomposition of a homeostatic cell state network into SCCs and the compartment representation, Equation 9.

(Left): An example cell state network representing the matrix in Equation 8 (self-links not displayed). The dashed circles denote the network’s Strongly Connected Components (SCCs ) (see definition in text). (Middle): The Condensed network is the corresponding network of SCCs, , wherein SCCs are the nodes and a link between two SCCs exists if any of their states are connected. For homeostatic networks, an SCC with dominant eigenvalue is at the apex, while other SCCs have . (Right): We distinguish two compartments, the Renewing compartment , consisting of the apex SCC, with , and the Committed compartment consisting of the remainder, with .

The stability of systems like Equation 7 is then determined by the dominant eigenvalues of each strongly connected component , for where is the number of SCCs, and their topological arrangement (the Perron-Frobenius theorem assures that for adjacency matrices of SCCs of cooperative systems, a unique, real, maximal eigenvalue exists, which is the dominant eigenvalue [Arrow, 1989; Greulich et al., 2019]). In brief, according to Greulich et al., 2019, the conditions for existence of a homeostatic state are that, at the apex of each lineage (the condensed cell state network), there must be an SCC with dominant eigenvalue , while all SCCs downstream of the former must have (see detailed discussion in the Appendix, 'Conditions for homeostasis'). Given this structure of homeostatic models, we can define two compartments in the cell state transition network: (1) the (self-) Renewing compartment (), which is the SCC at the apex of the lineage tree; and (2) the Committed compartment (), which consists of all SCCs with , that is, those downstream of the apex SCC. Importantly, cells in states forming have the potential to return to any state within the same compartment and this population maintains itself. Instead, the cell population in would vanish without external input, since the combined dominant eigenvalue of all those SCCs is negative (it is the maximum of all SCCs’ ), thus the progeny of each cell in the committed compartment will eventually be lost. We can thereby classify cells as being of a (self-)Renewing type (R) if their state is within , and of a Committed type (C) if their state is in . With this coarse-grained classification, a generic homeostatic model can be represented in terms of compartments and as,where the symbols above arrows are the effective rates of those events, denoting the average frequency at which they occur (loss events are not explicitly included, since they can be approximated by a short lived state in , as ). To be compatible with a homeostatic condition, it is further required that (i) the R-population remains on average constant (), that is, , and (ii) the loss rate of C must exceed its proliferation rate (), that is, . Figure 1 shows how a generic homeostatic cell state network can be condensed into an effective model of renewing and committed cell states, according to Equation 9. It has to be noted, however, that the events depicted in Equation 9 are not Markovian, that is, the timing of events is not independent from each other and depends on their history. Thus, the ‘rates’ , , , and γ are not constant rates in the Markovian sense, yet we can define them by the mean frequency of events occurring (see Appendix 'Approximation of generic GIA models' and 'Asymptotic clone size distributions: mathematical analysis').

The formulation in terms of renewing and committed states can help us to gain insights into potential behaviors of generic homeostatic cell fate models. In particular, we define generalized asymmetric divisions as events of the type , and generalized symmetric divisions as events of the type (symmetric renewal) and (symmetric commitment). With these definitions, we can categorize homeostatic cell fate models into two classes: Generalized Invariant Asymmetry (GIA) models are those which only exhibit divisions in the renewing compartment, while Generalized Population Asymmetry (GPA) are models for which such restriction does not hold. We note that the two classes are equivalently characterized by a conservation law: For GIA models, the number of cells in is strictly conserved, while for GPA models, no such conservation law holds. Since is necessary for conservation, the only possible conserved cell states in homeostasis are those in . Naturally, the previously discussed IA model is a GIA model and the PA model is a GPA model. Notably, the DH model (Equation 2) is of the GPA category, since in that model S and D cells form a single SCC at the apex of the lineage hierarchy, and thus they are both part of . Therefore, a division in the DH model, which is asymmetric in the conventional sense, corresponds to in terms of compartments (Equation 9) and thus it is a generalised symmetric division. According to this classification, PA and DH models are both in the same category (GPA), and indeed, both predict the same type of clone size distribution, an Exponential one (Greulich and Simons, 2016).

Numerical simulation of random cell fate models

To check whether the correspondence between model class, GIA vs. GPA, and predicted clonal statistics holds in general, we analyze the clonal dynamics numerically, by generating and testing a large number of random stochastic models, implemented via random generation of the parameters λ, ω, γ and . To simulate clones, we perform stochastic simulations based on the Gillespie algorithm (Gillespie, 1977), assuming a Markov process following the rules of Equation 3-5. We run, for each model, a large number of simulations with initially one cell in the compartment , thus the cell population of each simulation run represents one clone. Then we sample their outcomes, the total cell numbers per clone (the clone size) , to obtain predictions for clonal statistics, namely the frequency distribution of clone sizes (clone size distribution) and mean clone sizes (see Materials and methods).

We first study the mean clone size of surviving clones (with ), , shown in Figure 2, respectively, for the GIA and GPA models, as a function of time (the final time where is the minimal process rate, ). We note that indeed a common behavior is seen in each case. While for every simulated GIA model, saturates at a plateau value, it steadily increases for every GPA model. This is expected, and can be understood given that clones in a GPA model can go extinct while those in a GIA model not. Assume that there are initially a large number of clones, such that the total number of cells is . Since the system is homeostatic, it will reach a constant steady state after a sufficient amount of time, meaning that the mean clone size is . If no clones go extinct, as in GIA models, is constant and thus approaches a constant. However, in non-conserved multi-type branching processes, as GPA models are, the clone number decreases through progressive extinction of clones (Haccou et al., 2005), and therefore increases, despite the cell population as a whole staying stationary.

Figure 2.

Mean size of surviving clones, , as a function of time for random GIA models (a), and GPA models (b).

In (a), , in (b), is the time at 98% clone extinction. The grey shade represents the percentile of all the simulations (black lines limit the 5-95%ile range); the blue curves correspond to some illustrative selected simulations. Simulations for which the final mean is below two and where the final condition is not achieved (due to computational limitations) are not included: this results in 238 and 571 models, respectively for the GIA and GPA cases.

The resulting clone size distributions for the two model classes are shown in Figure 3. Here, clones sizes are rescaled by the mean value and compared to an Exponential distribution of unitary mean (green curve). As conjectured, all simulated GPA models shown in panel (b) predict asymptotically the same rescaled clone size distribution, namely a standard Exponential distribution. Deviations exist for small times and small clone sizes, but these deviations vanish in the large time limit (details on the convergence are shown in the Appendix, 'Analysis of the generalized Population Asymmetry model'). This means that different models within the GPA class cannot be distinguished in the long-term limit, since they differ only by the mean clone size, which is a free fit parameter. In analogy to statistical physics, we can categorize them as a universality class (Klein and Simons, 2011), meaning that the details of the model do not affect the (scaled) outcomes for assymptotic conditions, which is a form of weak convergence of random variables (Billingsley, 1968). However, the same cannot be said about the GIA models. In fact, we see all kind of shapes in the clone size distributions, both peaked distributions and non-peaked ones, and in fact, some distributions are even close to an Exponential form, and can thus not be distinguished from GPA models. The question is whether we can yet find other parameters for which, when large, also GIA models exhibit universality, that is, yield the same rescaled clone size distribution. For this purpose, we will in the following sections develop a deeper theoretical understanding of the model classes.

Figure 3.

Rescaled clone size distributions (expected relative frequency of clone sizes) for random GIA models (a), and GPA models (b), in terms of the rescaled clone size , at final time (see Figure 2 for definition).

The grey shade represents the percentile of all the simulations (black lines limit the 5-95%ile range); the blue curves correspond to some selected simulations. A reference curve corresponding to an Exponential distribution of unitary mean (’Exp(1)’) is shown in green.

Mathematical analysis: Markovian approximation of compartment model

To obtain a deeper understanding of the numerical results, we study the cell fate models in terms of the compartment representation, Equation 9. In this representation models are not Markovian, yet we can study their Markovian counterpart, as an approximation. While this is not expected to yield accurate clone size distributions in general, the limiting distributions of non-Markovian processes are commonly well estimated by their Markovian counterparts.

For GIA models, which only feature transitions between the renewing compartment, , and the committed compartment, , a corresponding Markovian model reads,in which represents a single state in and in , and symbols at arrows are the process rates. The number of cells in , , is conserved, that is, given an single -cell initially, it always remains at . Thus, we only need to consider the dynamics of cells in , . This Markov process can be solved analytically, and for sufficiently large steady state mean number of -cells, (see Appendix, 'GIA0 test case: steady state distribution and limiting behavior'), the rescaled distribution of cells in is,in which , and  and is the Gamma function (Abramowitz and Stegun, 1972). We note that this distribution exhibits a large variety of shapes: for large the distribution is peaked, while for small is loses its peak. Notably, for and , the distribution becomes Exponential and in this case it cannot be distinguished from the GPA case. On the other hand, for , that is, when the ratio of asymmetric divisions over the loss rate is high, this distribution tends to a Normal distribution with unitary mean and variance equal to . These different behaviors are graphically shown in the Appendix (see Appendix 1—figure 6, 7 and 8).

Appendix 1—figure 6.

GIA0 test case (see GIA0 test case: steady state distribution and limiting behavior') results in terms of steady state distribution of the the number of cells in state , .

The tested parameters correspond to the condition , as representative of the limiting case , and to different values of . The results from the numerical simulations are compared to the analytic solution (Equation 20), and its approximation, that is, the Poisson distribution (Equation 27).

Appendix 1—figure 7.

GIA0 test case (see 'GIA0 test case: steady state distribution and limiting behavior') results in terms of steady state rescaled distribution of the the number of cells in state , where .

The tested parameters correspond to the condition , as representative of the limiting case , and to different values of . The results from the numerical simulations are compared to the analytic solution (Equation 23), and its approximation that is the Gamma distribution (Equation 34).

Appendix 1—figure 8.

GIA0 test case (see 'GIA0 test case: steady state distribution and limiting behavior') results in terms of steady state rescaled distribution of the the number of cells in state , where .

The tested parameters correspond to the condition , as representative of the limiting case , and to different values of . The results from the numerical simulations are compared to the analytic solution (Equation 23), and its approximation that is the Normal distribution (Equation 40).

For the GPA models, a Markovian approximation reads, accordingly,whereby for homeostasis to prevail, and must hold. We note that the dynamics of are independent of and thus for the number of cells in in homeostasis holdswhich corresponds to a simple continuous-time branching process with two offspring, for which it is known that the resulting distribution of cell numbers is Exponential, that is, , where is the mean number of cells in the surviving clones (Haccou et al., 2005).

cells produced according to 12 follow the same fate as in the two-state GIA model above. While it is not assured that the distribution of cells is identical to that of Equation 11 (due to simultaneous production events of type ), we show in the Appendix, 'Asymptotic clone size distributions: mathematical analysis', that for large rates of production of C-cells, the distribution of C-cells – here: cells in state – attains a Normal distribution with mean equal to its variance . As each cell contributes independently to the production of -cells, we have that . Crucially, this means that in terms of the rescaled variable the standard deviation vanishes for large times, since . Hence, given fixed , can be approximated by a constant random number . Therefore, the rescaled distribution of the total number of cells is , where . Thus, the rescaled distribution of the total clone size, , is as well an Exponential.

Universality of generic cell fate models

For generic GIA or GPA models, the compartment representation, Equation 9, is not Markovian and one would not expect exactly the distributions we found in the previous section. Fortunately, the limiting distributions of non-Markovian processes and their Markovian counterparts are often, under certain conditions on the parameters, the same. While we reserve the technical arguments for the Appendix ('Asymptotic clone size distributions: mathematical analysis'), we note that this independence of the limiting distribution on the Markov property related to the central limit theorem, which does not rely on the Markov property.

To identify the correct limiting parameters for more complex cell fate models, we need to express the effective non-Markovian rates (i.e. the mean frequency of events) of representation nine in terms of the original model, 3–5. As discussed in the Appendix ('Approximation of generic GIA models' and 'Asymptotic clone size distributions: Mathematical analysis'), we identify those effective rates by the total rates of cell divisions, , , and where, for each compartment, is the probability of a single cell being in state X of , respectively ( are the solutions to Equation 6). In the Appendix, 'Asymptotic clone size distributions: mathematical analysis', we reason that all GPA models are expected to generate Exponential clone size distributions for large times t. This is indeed what is observed in Figure 3(b). Correspondingly, for GIA models we expect that for large the clone size distribution of GIA models would tend to a Normal distribution. To test this prediction, we simulated the same GIA models as for Figure 3 before, but we tuned parameters in such that the effective parameter becomes large (see details in the Appendix, 'GIA model for large '). The result is shown in Figure 4: for an illustrative case shown in panel (a), increasing changes the distribution from an exponential form to a peaked form akin to a Normal distribution, and for all simulated random GIA models, shown in panel (b), a Normal distribution is approached when becomes large.

Figure 4.

Rescaled clone size distributions (expected relative frequency P of clone sizes) for random GIA models as in Figure 3, at time (see definition in Figure 2).

Sensitivity to parameter is shown for one illustrative case in panel (a), and all GIA models for in panel (b). The distributions are shown in terms of the rescaled variables for panel (a) and , where is the distributions variance, in panel (b). In (b), the grey shade represents the percentile of all simulations (black lines limit the 5-95%ile range); the blue curves correspond to some selected simulations. A reference curve corresponding to a Normal distribution of zero mean and unitary variance is shown in green. Simulations for which is not reached (due to computational limitations) are not included, resulting in 922 model instances.

We note that when taking the limit of large , as shown in Figure 4, also all other process rates ω with i,j within increased as well. What if instead some process rates in do not scale to become large with ? To assess this situation, we studied a simple test case similar to model 10 but containing two states in , connected via direct state transition (see Appendix, 'GIA test case: bimodal distribution'). As discussed there, if all rates within are large compared to the rates in then indeed we observe a Normal clone size distribution, as expected. However, if the direct transition rates between the states of are smaller or of equal magnitude as γ, and in addition, one of the two division rates is higher then the other, then we observe a bimodal clone size distribution. The reason is that if the transitions between the two states in are rare compared to the life time of cells, , they become essentially separated and each of those states generate separate Normal distributions with different mean (due to different cell division rates in those two states) which, when overlaid, generate a bimodal clone size distribution (see detailed arguments in the Appendix, 'Asymptotic clone size distributions: mathematical analysis').

Finally, from those considerations follows:

GPA models attain an Exponential clone size distribution for time .

GIA models attain a Normal clone size distribution if all process rates within are much larger than the inverse lifetime of C-cells, γ.

Hence, the GIA and GPA model classes, each represent a universality class, that is, a scaling limit exists in which all models of the same class yield the same rescaled clonal statistics.

Discussion

Our analysis shows that intrinsic limitations exist for identifying strategies of stem cell self-renewal through clonal data from cell lineage-tracing experiments. This is due to different models of cell fate choice generating the same type of clonal statistics (clone size distributions), so that model inference based on clonal statistics – currently still the most prevalent method to determine stem cell self-renewal strategies – fails to distinguish them. The feature that different models asymptotically generate the same statistics is a form of weak convergence of random variables (Billingsley, 1968) and corresponds to universality, as known from statistical physics.

Cell fate models can in principle be very complex, with a plethora of cell (sub-)types in a tissue. We introduced a new categorization of cell types, distinguishing between cell states that are committed (C-cells), whose progeny is inevitably lost eventually, and non-committed or (self-)renewing cell states (R-cells), which retain the potential to remain or return to the apex of the lineage hierarchy. According to this categorization we classified generic models of cell fate choice as Generalized Invariant Asymmetry (GIA), if only generalized asymmetric divisions of the form occur for R-cells, and Generalized Population Asymmetry (GPA), when all kind of divisions can occur, as long as gain and loss of R-cells are balanced. Models of the GIA category are also characterized by a conservation law, since the number of R-cells is strictly conserved, while GPA models do not exhibit such a conservation law.

We found that the classification in GIA and GPA models mirrors the clonal statistics generated by them: models of the GPA class all generate clonal statistics which with time converge to an Exponential clone size distribution. Thus, two GPA models can therefore not be distinguished through clonal data, once some time has passed after induction of clones. For GIA models, distributions can generally vary, but if the rates of divisions and transitions in the compartment are much larger that the rate of cell loss, the clone size distribution of all those models becomes a Normal distribution. In that case, two GIA models can not be distinguished by the clonal data. While here we do not explicitly consider cell-extrinsic regulation of cell fate, this kind of regulation does not affect long-term clone size distributions, except when cells are arranged one-dimensionally (Klein and Simons, 2011; Bramson and Griffeath, 1980). Thus, our results cover cell dynamics in most renewing tissues, such as epithelial sheets or volumnar organs, but not (quasi-)one-dimensional arrangements of stem cells, as found in the seminiferous tubule, or in intestinal crypts, where clonal statistics may differ. Hence, our analysis shows that models of cell fate choice cannot in general be distinguished with further resolution beyond the R vs. C categorization of cell types. The universality of the model dynamics also shows that effective, simplistic models are often equally accurate to model experimental data, yet with a higher statistical power due to less free parameters.

While at first glance, this analysis seems to discourage efforts to unravel details of cell fate dynamics, room remains in regimes where the limiting conditions for asymptotic distributions are not fulfilled. In particular, if fast cycling committed progenitor cells are present, while stem cells are slow cycling, then the condition that the division rate of R-cells is much larger than the cell loss rate is not fulfilled. In that case, details of the model dynamics may affect the shape of the clone size distribution and thus allow distinction between models. However, caution should be given when an Exponential clone size distribution is observed, since this could indicate either a GIA model with high activity of committed progenitor cells, or a GPA model. In that case, the mean clone size needs to be consulted to distinguish models (see Figure 2). Differentiating between models within the GPA category is more difficult, since the predicted statistics from different models always become more similar over time. Short-term measurements would in principle allow such a distinction, but since in reality the underlying processes are not truly Markovian (as assumed for the modeling purpose) they are not necessarily a good representation of the real cell dynamics at short times. At long times, however, Markovian approximations are increasingly accurate, precisely because of the feature of universality.

How could the resolution of cell fate modeling be improved? The state-of-the-art approach to determine cell fate trajectories is via analysis and modeling of single-cell RNA-sequencing (scRNA-seq) data. However, many limitations to this method exist, discussed in Weinreb et al., 2018, and neither reversible trajectories nor the modes of cell division, such as asymmetric vs symmetric divisions, can be inferred. Intravital live imaging, on the other hand, allows to trace individual clones over time (Ritsma et al., 2012; Pittet and Weissleder, 2011; Hara et al., 2014; Rompolas et al., 2016), and thus can obtain details of cell fate trajectories, yet this technique is limited to few tissue types which are accessible for invasive long-term imaging. Nonetheless, while each of those experimental assays alone is prone to limitations in defining self-renewal strategies, advanced model inference schemes, that integrate data from different experimental sources, might be the way forward in the future to finally reveal the details of stem cell self-renewal strategies.

Materials and methods

The numerical analysis of the random cell fate model was implemented in Matlab. The description of the stochastic models definition, the random model generation and the simulation campaign is detailed in the Appendix, 'Stochastic process modelling'. Additionally, as a validation of the implemented simulator, based on the Gillespie algorithm (Gillespie, 1977), the IA and PA models were simulated and the results analyzed in the Appendix, 'Invariant Asymmetry and Population Asymmetry models'.

Analytical solutions were partially obtained using Mathematica.

In the interests of transparency, eLife publishes the most substantive revision requests and the accompanying author responses.

Acceptance summary:

This work builds on several recent high-profile applications of statistical physics to infer the mode of stem cell fate choice in mammalian tissues based on clone size distributions obtained via cell lineage tracing. The work makes a solid contribution especially in clarifying how a broad class of models conceptually and quantitatively fits with the population asymmetry or the invariant asymmetry model. Such universality argument may have been already taken for granted by physicists/mathematicians but may have not been obvious to biologists who have specific systems in mind that have much richer cell state structure than the original simplified models.

Through a detailed mathematical analysis, the authors demonstrate that this large set of models can, in homeostatic tissues, be decomposed into two 'universality classes' in terms of their predicted asymptotic clone size statistics. This limits our ability to qualitatively distinguish between different models of stem cell fate choice based on clone size data only.

Decision letter after peer review:

Thank you for submitting your article "Universality of clonal dynamics poses fundamental limits to identify stem cell self-renewal strategies" for consideration by eLife. Your article has been reviewed by two peer reviewers, and the evaluation has been overseen by a Reviewing Editor and Naama Barkai as the Senior Editor The reviewers have opted to remain anonymous.

The reviewers have discussed the reviews with one another and the Reviewing Editor has drafted this decision to help you prepare a revised submission.

We would like to draw your attention to changes in our revision policy that we have made in response to COVID-19 (https://elifesciences.org/articles/57162). Specifically, we are asking editors to accept without delay manuscripts, like yours, that they judge can stand as eLife papers without additional data, even if they feel that they would make the manuscript stronger. Thus the revisions requested below only address clarity and presentation.

Summary:

This manuscript provides an important contribution to the field of stem cell biology regarding their 'division modes' (symmetric vs. asymmetric divisions), as detailed in individual comments from reviewers.

Essential revisions:

The reviewers felt that the concerns can be mainly addressed textual changes as detailed in individual comments.

Reviewer #1:

This manuscript describes a general framework of understanding fate decision models, which have been previously proposed in explaining the dynamics of tissue stem cells in homeostasis. The work makes a solid contribution especially in clarifying how a broad class of models conceptually and quantitatively fits with the population asymmetry or the invariant asymmetry model. Such universality argument may have been already taken for granted by physicists/mathematicians but may have not been obvious to biologists who have specific systems in mind that have much richer cell state structure than the original simplified models.

The main criticism I have about the work is that the manuscript may misguide the readers to think that the GIA and GPA models are the only models that have been discussed in explaining the clonal statistics. As one of the authors showed in their previous paper (Greulich and Simons, 2016), the basic statistics such as the clone size distribution will change when adding spatial regulation to the voter model class, which is distinct from what is obtained by GIA or GPA. In fact, for all the tissues that have been studied so far, including the gut (10.1016/j.cell.2010.09.016, Lopez-Garcia et al., 2010), seminiferous tubule (Klein et al., 2010, Hara et al., 214, and 10.1016/j.stem.2018.11.013), and the skin (10.1016/j.stem.2018.09.005), experiments suggest that the spatial correlation between the fates are non-negligible and therefore the asymptotic statistics should follow the voter model. I would assume that there needs to be at least a justification in only studying the cases without spatial interactions here, especially if the authors do not know any specific example that nicely corresponds to their model.

Reviewer #2:

This work builds on several recent high-profile applications of statistical physics to infer the mode of stem cell fate choice in mammalian tissues based on clone size distributions obtained via cell lineage tracing.

The authors note that two conceptually distinct 'models' of stem cell fate choice, the Population Asymmetry (PA) and Dynamic Heterogeneity (DH) models, exhibit the same asymptotic clone size distributions. They then consider a more general set of possible models for the proliferation, differentiation and loss of a multi-type population of cells in a tissue.

Through a detailed mathematical analysis, the authors demonstrate that this large set of models can, in homeostatic tissues, be decomposed into two 'universality classes' in terms of their predicted asymptotic clone size statistics. This limits our ability to qualitatively distinguish between different models of stem cell fate choice based on clone size data only.

Overall I found this work, though technically complex, to be well motivated and well explained. The authors provide a significant advance to our understanding of how different descriptions of stem cell fate choice may be qualitatively distinguished experimentally (though, as the authors note, we may nevertheless apply parameter inference techniques and model selection to quantitatively distinguish competing models based on given datasets).

Reviewer #1:

[…] The main criticism I have about the work is that the manuscript may misguide the readers to think that the GIA and GPA models are the only models that have been discussed in explaining the clonal statistics. As one of the authors showed in their previous paper (Greulich and Simons, 2016), the basic statistics such as the clone size distribution will change when adding spatial regulation to the voter model class, which is distinct from what is obtained by GIA or GPA. In fact, for all the tissues that have been studied so far, including the gut (10.1016/j.cell.2010.09.016, Lopez-Garcia et al., 2010), seminiferous tubule (Klein et al., 2010, Hara et al., 2014, and 10.1016/j.stem.2018.11.013), and the skin (10.1016/j.stem.2018.09.005), experiments suggest that the spatial correlation between the fates are non-negligible and therefore the asymptotic statistics should follow the voter model. I would assume that there needs to be at least a justification in only studying the cases without spatial interactions here, especially if the authors do not know any specific example that nicely corresponds to their model.

We thank the reviewer for this remark. The reviewer notes correctly that in most tissues cell-extrinsic regelation via cell-cell signalling determines fate choices. This can indeed be represented by a (multi-type) voter model – a lattice model where lost cells are compensated by dividing neighbours. However, it has been shown that a voter model in any dimension larger than one has the same asymptotic rescaled clone size distribution as the corresponding branching process without interaction, which is the type of model we use (shown rigorously by Bramson and Griffeath, 1980, and discussed in the context of cell fate dynamics by Klein and Simons, 2011). For mean clone sizes as function of time t, there is no difference for dimensions larger than two (in the two dimensional scenario merely a logarithmic prefactor distinguishes the mean clone size in both models; since this factor is nearly constant for large t, it often cannot be distinguished by experimental data (Klein and Simons, 2011)).

Therefore, while a voter model may be a more realistic depiction of a tissue with spatial cell-cell interaction, the asymptotic rescaled clone size distribution of such a model is the same as for a corresponding model without spatial regulation, if considering tissues with cell arrangements of dimensions larger than one. Hence, our results would not be different for a corresponding voter model where spatial regulation is included, for two- or three-dimensional tissues such as epithelial sheets and volumnar tissues (e.g. liver or stroma). Yet, we agree that for quasi one-dimensional arrangements of stem cells, as found in seminiferous tubules and intestinal crypts, the clone size distribution differs, and these are indeed not covered by our model.

In the original version of the manuscript we expressed those limitations after the model definition: "This model is a general multi-type branching process, which is suitable to describe cell population dynamics in any dimension larger than one, even under cell-extrinsic regulation (Klein and Simons, 2011; Bramson and Griffeath, 1980)", but we acknowledge that this point warrants a clearer explanation. We therefore now include at this point a paragraph that clarifies the validity range of results to be expected from our model. Furthermore, we now discuss this in the Discussion section.

Appendix 1—table 1.

IA and PA test cases simulation parameters (see 'Invariant Asymmetry and Population Asymmetry models').

Case	λ	γ	r
IA#1	1.0	1.0	-
IA#2	2.0	1.0	-
IA#3	5.0	1.0	-
PA#1	1.0	1.0	1/4
PA#2	2.0	1.0	1/4
PA#3	2.0	1.0	1/6

Appendix 1—table 2.

GIA test case simulation parameters (see 'GIA test case: bimodal distribution').

Case	ω^	λ_1/λ_2
GIAB#1	3 101	1
GIAB#2	3 10-2	1
GIAB#3	3 102	10
GIAB#4	3 101	10
GIAB#5	3 100	10
GIAB#6	3 10-1	10
GIAB#7	3 10-2	10