Quantification and statistical analysis—alternative models and parameter-fitting

To fit the model to the experimental data, we constructed distributions of clone size and composition numerically using stochastic simulations of at least 10,000 clones, each derived from a single bipotent progenitor cell labelled at the initial time point. We made use of standard bootstrapping methods on the lineage tracing data to build standard deviations, as well as 95% confidence intervals for our predictions. To constrain and test the model, we focussed on the E12.5–P14 R26-CreERT2/R26R-Confetti tracings, for which the largest data set was available, and where islet clones were resolved with high statistical confidence. However, to further challenge the model, we used the fits to predict clone size and compositional dependences using the E15.5–P14 and E18.5–P14 R26-CreERT2/R26R-Confetti tracings, as well as the population-level Ngn3-CreER^TM tracings. To further constrain the model, we made use of several robust features of the data:

First, we noted that the overall α- and β-cell distributions were well-described by a simple exponential size dependence (Fig. 4h, j, k and Supplementary Fig. ^5m) (as was the overall clone size distribution), and that α- and β-cell unipotent clones also had similar sizes (Supplementary Fig. ^5k). This is not what one would expect in models where the division rate of β-cells would be much larger than α-cells, suggesting that λ_A(t) ≈ λ_B(t). Simulating very different proliferation rates for α- and β-cells produced a marked departure from exponentiality (see Supplementary Fig. ^3k, where n/N = 0.333 and λ_A = 4λ_b). Indeed, the higher abundance of β-cells over α-cells overall could be explained by the production of a larger number of unipotent β-cell clones as compared to unipotent α-cell clones, and larger numbers of β-cells in bipotent clones (Fig. 4f and Supplementary Fig. ^{3a, b}). This situation prevailed in all tracings, confirming that differences in α- and β-cell production was not due to temporal differences in differentiation rates. This argues that that the rate of transition into the β-cell sublineage is larger than into the α-cell sublineage, k_B > k_A. This observation was further confirmed by the E12.5–P28 tracing (Supplementary Fig. ^5j–o), where the phenomenology was reversed: α-cell clones had not increased in size since P14, while β-cell clones had more than doubled in size, resulting in exponential distributions with different slopes. This was indicative of a difference in the division rates, λ_A(t) < λ_B(t), appearing between P14 and P28, i.e. that the expansion post-P14 of the β-cell compartment becomes driven by preferential proliferation, and not preferential differentiation of bipotent progenitors P into β-cell precursors B, as earlier on. Here, we note that we are able only to refer to effective proliferation rates: If, for example, there were an undisclosed residual base-line rate of cell apoptosis, numerical simulations shows that the join clone size distribution would be largely unaffected (Supplementary Fig. ^3l), while the true proliferation rate would be proportionately higher than the inferred effective rate to accommodate the apoptosis rate.

Second, as time is counted here in terms of numbers of divisions, the relative division rate of the bipotent progenitor population is irrelevant and, for simplicity, can be set equal to that of unipotent precursors, and taken as constant: λ ≈ λ_A ≈ λ_B (with the total number of divisions and timing of commitment/differentiation being the only relevant quantity).

Third, as mentioned above, if bipotent progenitors and unipotent precursors coexisted in the long term and had highly different dynamics (for instance in their division rate), one would expect to see more complex clone size distributions, which is ruled out by the simple exponential scaling form for all time points. To show this, given that roughly one half of the clones traced from E12.5 displayed bipotent outcomes, we simulated a scenario in which 50% of induced clones were unipotent precursors (either α- or β-cell lineage in a 1:2 ratio), and 50% were bipotent until the end of the tracing. However, if both populations had the same division rate, a fraction of clones derived from bipotent progenitors would still end up with unipotent output due to chance and small number statistics (for instance, a two-cell clone make derive from the stochastic specification of β-cells), and the overall bipotency fractions would be too low to infer bipotency with confidence. The only solution to remedy this is to assume a much larger proliferation rate of bipotent progenitors (so that their clone size is so large as to make chance unipotent outcomes unlikely), while maintaining the overall average clone size of 5.6 cells. We thus simulated stochastic unipotent precursors undergoing N = 0.75 divisions and stochastic bipotent progenitors undergoing N = 3.25 divisions. However, this resulted in non-exponential clone size distributions, which despite the higher number of fitting parameters compared to our previous simpler model, worsened the fit in particular for the total clone size (Supplementary Fig. ^3m).

Thus, the simplest model, consistent with the qualitative features described above (from E12.5 to P14) is one where (i) division rates are constant, and (ii) a transition takes place from bipotency to unipotency (viz. rates k_A(t) = k_B(t) = 0 for t < n/N, and k_A(t) = K_A and k_B(t) = K_B = 1 − K_A for t > n/N, with constant probabilities K_A (and therefore K_B) chosen to match the proportions of α-cells and β-cells found in the organ overall.

From a more quantitative standpoint, the simplest model thus relies on three key parameters: the relative probability to differentiate into α- vs. β-cell fate, k_A, the total number of cell divisions between E12.5 and P14,N, and the timing of the transition n/N (expressed as a fraction of N). The first and second parameters can be estimated with high precision directly from the data simply by (i) calculating the total ratio of α to β-cells, which leads to the estimate KA=0.37=1−KB, and (ii) using the average clone size of 6.1 ± 0.6 cells to estimate N = 2.6 ± 0.1 (mean ± SEM). The third parameter needs to be fit via numerical simulations of the model. To constrain the model with the smallest fraction of data (in order to test for consistency on the remainder of the clonal data), we thus resorted to a simple fit of the fraction of bipotent clones (which increases monotonically with the timing of commitment, n/N). We performed 1000 iterations of the numerical simulation, each time boot-strapping the experimental dataset (to calculate the confidence intervals on bipotency fraction), and found the best fit for n/N (each time simulating 10,000 clones). This allowed us to find a best fit and construct a 95% confidence interval with n/N = 65% ± 15%.

With this best-fit value, several aspects of the data, not used in the fitting procedure, can be predicted quantitatively by the model, including the joint clone size distributions of α and β-cell compartments (Fig. 4i). Features of this two-dimensional distribution include: (i) the exponential-like clone size distributions of α- and β-cell compartments (Fig. 4h); (ii) the relative fraction of unipotent α- and β-cell clones (Fig. 4f); (iii) the fact that α- to β-cell ratio is identical when looking either at bipotent clones only or α/β-cell unipotent clones only (Supplementary Fig. ^3a); and iv) the correlation and variance of α vs. β-cell number in a clone (Supplementary Fig. ^{3d, f}), showing that the variance increases with the number of cells in a clone, as expected for a binomial distribution arising from random fate assignment between α- and β-cell fate. However, we note that the E12.5–P14 tracing data contained 6 clones (out of 138), which were manifest outliers compared to the model and the rest of the data. These 6 clones all had an average clone size between 3 and 5 standard deviations away from the mean. To improve visual comparison between joint distributions of α- and β-cells, we thus “zoomed in” (in Fig. 4i, compared to Supplementary Fig. ^3d) on smaller clones sizes to concentrate on the bulk of the data. Interestingly, despite their size, these outlier clones were “normal” in their ratio of α- and β-cells (i.e. an approximate ratio of 1:2, Supplementary Fig. ^3d). Given that we did not find these clones to correlate with areas of larger clonal induction, we reasoned that they were unlikely to arise from fusion events, and might thus arise due to rare hotspots of proliferation during early embryogenesis. This would be consistent with the fact that such outliers were not present in the E15.5 or E18.5 tracings.

To perform a sensitivity analysis, the predictions were contrasted with a model involving earlier sublineage commitment (i.e. mostly unipotent α- and β-cell clones being present, and thus with highly independent α- vs. β-cell compartment sizes, see Supplementary Fig. ^3g), or later sublineage commitment (where α- and β-cell compartment size in a clone are maximally correlated, as both derive from the same bipotent progenitor pool, see Supplementary Fig. ^3h). Notably, we found an intermediate level of positive correlation between α- and β-cell numbers in E12.5 clones (slope 0.3 ± 0.04, R² = 0.27), in quantitative agreement once again with the predictions of the model (Fig. 4i and Supplementary Fig. ^{3d, i}). Alternative models, such as one where a bipotent islet progenitor self-renews while continuously producing unipotent α- and β-cells with stochastic fate allocation, provided poor fits to the data (Supplementary Fig. ^3j), including the average and variance of the number of α-cells in clones with a given number of β-cells (Supplementary Fig. ^{3d, f})

Although our data ruled out distinct phases of α- and β-cell allocation (a feature confirmed by the Ngn3-CreER^TM tracing), one can, of course, introduce additional refinements of the minimal model incorporating, for instance, a defined time for the fate specification of α- and β-cells, as well as small time delays between the allocation of both. Although these additional parameters only improve marginally the quality of the fits, the data can accommodate, for instance, a phase of α-cell allocation (initiating at around T_A ≈ 50%) and β-cell allocation (starting later at around T_B ≈ 65%). Note that these delays are within the confidence intervals of n/N from the minimal model, and that the model still requires the late specification of α-cells (50% corresponds to around E15.5–E18.5 in embryonic time from the clone size of the different tracings) compared to past proposals in the literature. Finally, the same model as for the E12.5–P14 tracing could explain faithfully the E12.5–P28 tracing (Supplementary Fig. ^5j–o) by simply assuming that α-cells do not proliferate within the later time interval (P14–P28), while β-cells undergo an additional round of division. This is consistent with clone sizes of each compartment (Supplementary Fig. ^5k), as well as overall clonal potency (Supplementary Fig. ^5l), the α- and β-cell clone size distributions (Supplementary Fig. ^5m), and the joint probability distributions (Supplementary Fig. ^{5n, o}).