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

Abstract

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.

Keywords: cell fate; clonal data; computational biology; developmental biology; differentiation; none; stem cell self-renewal; stochastic modelling; systems biology; universality.

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 $A$ 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, $S_k$, 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 $\mu =0$ is at the apex, while other SCCs have $\mu <0$. (Right): We distinguish two compartments, the Renewing compartment ${\displaystyle \mathcal{ℛ}}$, consisting of the apex SCC, with $\mu =0$, and the Committed compartment ${\displaystyle \mathcal{𝒞}}$ consisting of the remainder, with $\mu <0$.

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

In (a), $\tau =20/{\alpha }_{\min }$, in (b), $\tau$ 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.

Figure 3.. Rescaled clone size distributions (expected relative frequency P of clone sizes) for random GIA models (a), and GPA models (b), in terms of the rescaled clone size x=n/n¯s, at final time t=τ (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.

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

Sensitivity to parameter ${\widehat{\lambda }}_R$ is shown for one illustrative case in panel (a), and all GIA models for ${\widehat{\lambda }}_R=30$ in panel (b). The distributions are shown in terms of the rescaled variables $x=n/{\overline{n}}_s$ for panel (a) and $\tilde{x}=(n-{\overline{n}}_s)/{\sigma }_n$, where ${\sigma }_n$ 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 $t=\tau$ is not reached (due to computational limitations) are not included, resulting in 922 model instances.

Appendix 1—figure 1.. Invariant Asymmetry (IA) test cases clone size distribution P⁢(n), that is the distribution of the total number of cells n forming the progeny of a single initial cell in ℛ.

For each case, the distribution is shown at τ (defined in Figure 2, main text), which is well representative of the steady state condition. Tested parameters for cases IA#1-3 are provided in Appendix 1—table 1; the numerical simulation results are compared to the expected Poisson distribution. The detailed discussion is reported in 'Invariant Asymmetry and Population Asymmetry models'.

Appendix 1—figure 2.. Population Asymmetry (PA) test cases clone size distribution P⁢(n), that is the distribution of the total number of cells n forming the progeny of a single initial stem cell.

For each case, the distribution is shown at the final time τ, at which the total extinction of the process is not yet achieved. Tested parameters for cases PA#-3 are provided in Appendix 1—table 1; the numerical simulation results are compared to the solution of the numerical integration of the master Equation 8 and, for test case PA#1, also to the reference analytic solution from Antal and Krapivsky, 2010. The detailed discussion is reported in 'Invariant Asymmetry and Population Asymmetry models'.

Appendix 1—figure 3.. Metastate (MS) test cases simulation results in terms of mean number of cells in the surviving clones n¯s and extinction probability P(n=0) as function of time (scaled by the final simulation time τ).

As well as for the PA test cases, at τ the total extinction of the process is not yet achieved. Profiles from the numerical simulation for cases MS#,3 are compared to the corresponding PA#1,3 test cases which are based on parameters provided in Appendix 1—table 1. The detailed discussion is reported in 'Population Asymmetry model using metastates'.

Appendix 1—figure 4.. Metastate (MS) test cases simulation results in terms clone size distribution P⁢(n), that is the distribution of the total number of cells n forming the progeny of a single initial stem cell.

As well as for the PA test cases, the distribution is shown at the final time, τ, at which the total extinction of the process is not yet achieved. Profiles from the numerical simulation for cases MS#,3 are compared to the corresponding PA#1,3 test cases which are based on parameters provided in Appendix 1—table 1. The detailed discussion is reported in 'Population Asymmetry model using metastates'.

Appendix 1—figure 5.. GIA⁰ test case parameters λ^1 and λ^2 over the contour map of the expected steady state mean number of cells in state X2, n¯2*.

The tested conditions are divided in three groups representing the limiting behaviors discussed in in 'GIA⁰ test case: steady state distribution and limiting behavior', and for which the steady state distribution is shown respectively in Appendix 1—figure 6, Appendix 1—figure 7 and Appendix 1—figure 8.

Appendix 1—figure 6.. GIA⁰ test case (see GIA⁰ test case: steady state distribution and limiting behavior') results in terms of steady state distribution P*⁢(n2) of the the number of cells in state X2, n2.

The tested parameters correspond to the condition ${\widehat{\lambda }}_2=0.01$, as representative of the limiting case ${\widehat{\lambda }}_2\to 0$, and to different values of ${\widehat{\lambda }}_1$. 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.. GIA⁰ test case (see 'GIA⁰ test case: steady state distribution and limiting behavior') results in terms of steady state rescaled distribution P*⁢(x2) of the the number of cells in state X2, where x2=n2/n¯2*.

The tested parameters correspond to the condition ${\widehat{\lambda }}_2=0.99$, as representative of the limiting case ${\widehat{\lambda }}_2\to 1$, and to different values of ${\widehat{\lambda }}_1$. 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.. GIA⁰ test case (see 'GIA⁰ test case: steady state distribution and limiting behavior') results in terms of steady state rescaled distribution P*⁢(x2) of the the number of cells in state X2, where x2=n2/n¯2*.

The tested parameters correspond to the condition ${\widehat{\lambda }}_1=60$, as representative of the limiting case ${\widehat{\lambda }}_1\to \mathrm{\infty }$, and to different values of ${\widehat{\lambda }}_2$. The results from the numerical simulations are compared to the analytic solution (Equation 23), and its approximation that is the Normal distribution (Equation 40).

Appendix 1—figure 9.. GIA random models (generated as described in 'Generation of random models' and analyzed in the main text) equivalent parameters λ^1=λ^R and λ^2=λ^C (see section Approximation of generic GIA models) over the contour map of the expected steady state mean number of cells in the committed compartment, n¯C*.

Some illustrative cases, for which the steady state distribution is shown in Appendix 1—figure 11, Appendix 1—figure 12 and Appendix 1—figure 13, are highlighted.

Appendix 1—figure 10.. Relative error of the the equivalent model approximation, ϵ, (see definition in 'Approximation of generic GIA models') as function of λ^2=λ^C for the GIA random models (generated as described in 'Generation of random models' and analyzed in the main text).

The selected cases correspond to some illustrative cases for which the steady state distribution is shown in Appendix 1—figure 11, Appendix 1—figure 12 and Appendix 1—figure 13.

Appendix 1—figure 11.. GIA random models selected cases (see Appendix 1—figure 9 and Appendix 1—figure 10) where λ^2<0.2: the steady state distribution P*⁢(nC) of the number of cells in the committed compartment, nC, is compared to that of the equivalent model (Equation model in the legend) analytic solution and its approximation for low λ^2 (i.e. the Poisson distribution, Poisson (λ^1)).

Results discussion is reported in 'Approximation of generic GIA models'.

Appendix 1—figure 12.. GIA random models selected cases (see Appendix 1—figure 9 and Appendix 1—figure 10) where λ^2>0.8: the steady state rescaled distribution P*⁢(xC) of the number of cells in the committed compartment, where xC=nC/n¯C*, is compared to that of the equivalent model (Equation model in the legend) analytic solution and its approximation for high λ^2 (i.e. the Gamma distribution, Gamma (λ^1,1/λ^1)).

Results discussion is reported in 'Approximation of generic GIA models'.

Appendix 1—figure 13.. GIA random models selected cases (see Appendix 1—figure 9 and Appendix 1—figure 10) where 0.2<λ^2<0.8: the steady state distribution P*⁢(nC) (or the rescaled distribution P*⁢(xC)) of the number of cells in the committed compartment, nC (or in the rescaled case xC=nC/n¯C*), is compared to that of the equivalent model (Equation model in the legend) analytic solution.

Results discussion is reported in 'Approximation of generic GIA models'.

Appendix 1—figure 14.. Rescaled clone size distribution for the random GIA models when λ^R=30 at the final simulation time, which corresponds to 20/αmin (αmin is the minimum process rate).

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. A reference curve corresponding to a Normal distribution of unitary mean and variance equal to $1/{\widehat{\lambda }}_R=1/30$ is shown in green. Distributions of the total number of cells $n$ are scaled by the mean number of cells $\overline{n}$, being $x=n/\overline{n}$. Simulations for which the final condition (20 times the inverse of the minimum process rate) is not achieved (due to computational limitations) are not included, resulting in 922 models. Results discussion in provided in 'GIA model for large'.

Appendix 1—figure 15.. Rescaled distribution of the cells number in the committed compartment in the G⁢I⁢AB test cases at time τ, which is 20/αmin (αmin is the minimum process rate).

The distributions $P({\tilde{x}}_C)$ of the number of cells in the committed compartment $n_C$ is rescaled considering that ${\tilde{x}}_c=(n_C-{\overline{n}}_C)/{\sigma }_{n_c}$, where ${\sigma }_{n_c}$ is the variance of $n_c$. In addition to the stochastic simulation results for different settings (see Appendix 1—table 2), the reference Normal and bimodal distributions are also shown. Results discussion is provided in 'GIA^B test case: bimodal distribution'.

Appendix 1—figure 16.. Clonal size distribution (corresponding to the 50%ile curve) in the GPA random models at different extinction fraction (i.e. different time).

The curves are compared to the expected Exponential distribution (see 'Analysis of the generalized Population Asymmetry model').