Quantifying Lgr5-positive stem cell behaviour in the pyloric epithelium

Using in-vivo lineage tracing data we quantified clonal expansion as well as proliferation and differentiation of the Lgr5-positive stem cell population in pyloric gastric glands. Fitting clone expansion models, we estimated that there are five effective Lgr5-positive cells able to give rise to monoclonal glands by replacing each other following a pattern of neutral drift dynamics. This analysis is instrumental to assess stem cell performance; however, stem cell proliferation is not quantified by clone expansion analysis. We identified a suitable mathematical model to quantify proliferation and differentiation of the Lgr5-positive population. As expected for populations in steady-state, the proliferation rate of the Lgr5-positive population was equal to its rate of differentiation. This rate was significantly faster than the rate at which effective cells are replaced, estimated by modelling clone expansion/contraction. This suggests that the majority of Lgr5-positive cell divisions serve to renew epithelial cells and only few result in the effective replacement of a neighbour to effect expansion to the entire gland. The application of the model under altered situations with uncoupled differentiation and proliferation was demonstrated. This methodology represents a valuable tool for quantifying stem cell performance in homeostasis and importantly for deciphering altered stem cell behaviour in disease.


Generalizing clonal expansion models for any initial clone size
Lopez-Garcia (2010) 1 , described clonal expansion in the small intestinal crypt as a stochastic process where the probabilistic dynamics of the size of the clone is that of a one-dimensional random walk with closed boundaries. This model hypothesises that the size of the clone moves at a constant rate, λ, with steps of size 1 between 0 and N cells, both of which are absorbent states 1-3 The forward stochastic differential equations for a one-dimensional random walk with absorbent boundary states are as follows:  Where P ij (t) = P{S (t) = j | S (0) = i} is the probability that the clone size at time t, S (t), is equal to j given that its initial size at time 0 was equal to i, 0 < i < N. N is the number of cells in the ring and λ is the rate at which clones gain or lose cells.
Therefore clonal study involves only those N cells in the gland base able to give rise to clonal expansion. By hypothesis, the clones are equipotent and follow a neutral drift dynamics so that the probability of gaining a cell is equal for all of them and equal to the probability of losing a cell. Monoclonality occurs when a given clone reaches size N, while clone takes place when the clone size reaches the other absorbent state, i.e. 0 cells.
Equations (S1) have been explicitly solved for the special case in which the initial size of the clone is 1 cell 1, 3 . We have extended the solution for other initial states so that the model can be fitted to data in which the initial clone size is greater than 1 cell after tamoxifen activation.
To find the solution for any initial clone size, we started from the solution of Equations (S1) for 1 ≤ j ≤ N-1 3 : Where A k is the constant of integration.
Considering that the size of a clone at the initial time, S (0), is equal to i, which is between 1 and N-1, the initial value for the stochastic differential Equations (S1) are defined as follows: With this initial value, A k in Equation (S2) has the form: By substituting the value of A k in Equation (S2), the explicit solution for P ij (t), from Equations (S1) can be found. Replacing in Equations (S1) with the solution for P ij (t) in (S4), P i0 (t) and P iN (t) are found by integration. Adding the transition probabilities of the absorbent states, 0 and N, the whole set of solution is then:  In the stomach, we have observed that after tamoxifen activation, the initial clone size may not be a constant value but vary from gland to gland. Thus, the variation in the initial number of labelled Lgr5-positive cells between clones has been included in the model by assuming that the initial clone size, S(0), has a given The initial value for Equations (S1) is now And the value of A k is equal to: With this value for the integrating constant, the solution of Equations (S1) is The same result can be obtained using the law of total probability as follows In the newly generated dataset in this work, the distribution of the initial size of the clones immediately after tamoxifen action was in good agreement with a Poisson The maximum likelihood estimator for the parameter α of the Poisson distribution is Total clones of Lgr5+ cellŝ ln Unlabelled clones of Lgr5+ cells When no unlabelled clones are detected, α can be estimated as follows: Assuming in this situation a Poisson distribution as initial clone size distribution will underestimate the probability of clones with size from 1 to N, while it will give high values to the probability of clones with a non-possible initial size, i.e. greater than N.
In that case a binomial distribution with parameters N, number of cells able to give rise to monoclonality and p, probability of having a cell labelled after Cre activation, could be used. This distributions are closely related so that when N/p is a large number, the binomial distribution converges to the Poison distribution with parameter α = N . p; Similarly, a superposed Poison distribution can be replaced by a multinomial distribution for mouse reporter models with more than one colour label induced with relatively high doses of Tamoxifen.
Other distribution functions as well as the empirical distribution derived from the histogram of observations could be used to define the distribution of the initial clone size. The empirical distribution of the initial clone size, S (0), was defined as: The parameters λ and N were estimated by fitting equations (S5) and (S8) to datasets, simulated or experimentally observed, respectively, by MCMC methods.
The chosen distribution function for the likelihood of the data was a multinomial with parameters 0 ( ).. ( ) i iN P t P t which are described in Equations (S5) and The selected non-informative prior distributions for N and λ were: where E [λ] = g 1 / g 2 and Var [λ] = g 1 / g 2 2 Prior to fitting the clonal model to the dataset generated in this work , is the probability of reaching a clone size equal to k when t → ∞ given that the initial clone size was equal to i. The limiting probability of reaching a given clone size is equal to 0 for all cases with the exception of 0 and N which means that in the long term clone extinction and monoclonality are the only possible events and therefore their limiting probabilities add up to 1. These limiting probability functions depend on the number of Lgr5-positive cells at the base of the gland, N, and on the initial clone size, i. They are not affected by the clone replacement rate, λ.

Parameter estimation
The limiting distribution assuming that the initial clone size varies in between glands according to a distribution, F, as described above is as follows: As above any other suitable distribution function as well as the empirical observed distribution could be used.
The probability of extinction,   Figure 1A shows that the true value of the parameters N and λ was accurately recovered only when the initial clone size assumed in the model was coincidental with the initial size used to simulate the data.
The model assuming 1 initial cell per clone resulted in accurate estimations of the parameters N and λ only when fitted to the dataset generated with an initial clone size equal to 1 cell (Supplementary Figures 1A and 1B). The greater the initial clone size used for data generation, the larger the bias of the estimates of N and λ obtained by fitting the model developed with 1 initial cell per clone (Supplementary Figures 1A   and 1B). Supplementary Figure 1C shows that the dataset generated with 4 initial cells per clone was accurately described by the model that assumed that initial value, while the performance of the model developed with 1 single initial cell was very poor (Supplementary Figure 1D).
Moreover, the importance of using an accurate measurement of the initial clone size to fit the probability of extinction was also demonstrated by fitting Equation (2) to the dataset in Figure 3A but assuming different values for the initial clone size.
Data in Supplementary Figure 3E shows that the initial clone size affects remarkably the value of the estimation of N. For instance, when the initial clone size was assumed

Birth-Death stochastic model to quantify population kinetics of Lgr5positive cells
Proliferation and differentiation of Lgr5-positive cells was modelled as a birth death stochastic process. The forward differential stochastic equations are 5 The parameters δ and µ were estimated by MCMC methods. The chosen distribution function to construct the likelihood of the data was a multinomial with parameters 0 ( ).. ( ) n P t P t which are described in Equations (S17) and (S18) Fitting diagnosis and posterior distribution were analysed for δ, µ and µδ.

Parameters estimation
In order to estimate η, the model was fitted to data by the MCMC approach Lgr5-positive cells located at the base of the gland were considered the effective cells giving rise to monoclonal crypts and therefore they could replace each other and move across the gland base. After division of one of these cells, the cell to be replaced was either the cell located in the ring immediately above with probability 1-P or one of the other effective Lgr5-positve cells with probability P ( Figure 4A). Two scenarios were implemented for the spatial replacement of these cells ( Figure 4A). In the first scenario i) cell replacement affected either the cell on the left or on the right with the same probability, i.e cells move along a one-dimensional ring; in scenario ii) cell replacement takes place in two dimensions, so that any of the other effective Lgr5 cells have the same chance to be replaced. In both scenarios, the just replaced cell could either move upwards with probability 1-P or could replace one of the other effective Lgr5-positive cells. A necessary restriction is that one cell cannot be affected twice in the same division event.
Cells located out of the 5 basal effective positions move always upwards to the ring immediately above. Vertical cell displacement spreads upwards until the last ring modelled in the gland.