Critical growth of cerebral tissue in organoids: theory and experiments

We develop a Fokker-Planck theory of tissue growth with three types of cells (symmetrically dividing, asymmetrically dividing and non-dividing) as main agents to study the growth dynamics of human cerebral organoids. Fitting the theory to lineage tracing data obtained in next generation sequencing experiments, we show that the growth of cerebral organoids is a critical process. We derive analytical expressions describing the time evolution of clonal lineage sizes and show how power-law distributions arise in the limit of long times due to the vanishing of a characteristic growth scale. We discuss that the independence of critical growth on initial conditions could be biologically advantageous.

Introduction The mechanisms of tissue growth and renewal are a core topic of stem-cell research [1,2].In particular, the role of stochasticity in cell differentiation is discussed [3][4][5][6][7][8].Single cell sequencing [9], combined with the labeling of cells with inheritable DNA sequences, enables large scale, quantitative studies of cell populations in biological tissues, where offspring populations can be traced back to their individual ancestral cells [10,11].Such lineage tracing experiments have revealed that offspring numbers in mammalian cerebral tissue can vary by several orders of magnitude [6,12], which supports the hypothesis that stochasticity is an important property of cell proliferation and differentiation in the developing cerebral cortex.
This work presents a study of lineage tracing data obtained by sequencing 15 cerebral organoids at different stages of their development [12].Cerebral organoids are highly controllable, self organized in vitro models of the human cerebral cortex grown from stem cells.Organoids are unique because they model human tissues which cannot be studied in vivo.They have therefore become important biological tools to study neural development and brain diseases [13][14][15].We take a physics point of view on the population dynamics of cell lineages in cerebral organoids, show that organoid growth is a critical process, and discuss the biological implications.
Model Our study begins with the observation that the numbers of descendants of an individual stem cell in the organoid (lineage sizes) are roughly distributed according to a 3/2-power-law.This behavior becomes more and more pronounced at late stages of the organoid development.To model the growth process mathematically, we build upon the theory of continuous state branching processes [16,17], pioneered by Feller [18].These processes are known to lead to power-law distributed poluation sizes [19][20][21].Stem cells (green, S) either divide at rate gS (S→2S), differentiate at rate gA (S→A) or die at a rate g0 (S→Ø).Differentiated cells (orange, A) that committed to a developmental trajectory either divide asymmetrically (A→A+N) at a rate gA, or differentiate directly into non-dividing N cells (rate gF ).On average, each A cell produces N N-cells until it looses the ability to divide.At criticality, the rates gS (division) and gA (differentiation) are equal.
We thus introduce the SAN model.It consists of three agents: symmetrically dividing S-cells that represent stem cells, asymmetrically dividing A-cells and nondividing N-cells (fully developed cells, e.g.neurons).Scells undergo symmetric division (S→2S) at rate g S , differentiation (S→A) at rate g A and death (S→0) at rate g 0 .A-cells have committed to a developmental trajectory and produce N-cells through asymmetric divisions (A→A+N) at a rate g N until the process is terminated by direct differentiation (A→N) (rate g F ).The branching process of symmetric division and differentiation of S-cells with rates g S and g A is at the heart of the SANmodel.Criticality is reached, when the two rates are equal: g S = g A .The model is illustrated in Fig. 1.We solve the SAN-model analytically in the continuum limit and show how, at long times, 3/2-power-law distributions of cell populations asymptotically arise near criticality.Fitting the model predictions to the empirical data of Ref. [12], we show that cerebral organoid growth is indeed critical (see Fig. 2 and Table I).In experiments, organoids are grown for forty days.As we will demonstrate, our model can describe tissue growth as a dynamical process with a limited number of parameters -three rates and one initial condition.
Fokker-Planck description of lineage dynamics We start with a master equation for the SAN process.Individual populations of S, A and N-cells are not accessible in experiments, since all descendants of a stem cell inherit the same lineage identifier.We therefore need to calculate the probability distribution of the total lineage size x tot = s + a + n.This calls for a slight simplification.We replace the processes A→A + N and A→N by the assumption that each A-cell produces N = 1 + g N /g F N-cells over the course of its existence (where the +1 stems from the final A→N conversion) (Fig. 1).Since the N-cell output per A-cell varies by multiple orders of magnitude less than the total offspring of an S-cell [22], this simplification does not affect the model's main predictions, making it analytically tractable.N is a fitting parameter in our theory.The total lineage size becomes x tot = s + n.The master equation for the probability distribution f of S-and N-cell numbers s, n at time t then reads The right hand side terms of Eq. ( 1) correspond to the different processes that the cells undergo: S-cell death at rate g 0 , S-cell division at rate g S and differentiation into N N-cells at rate g A .Focusing on large cell counts, we translate the discrete process of Eq. ( 1) into a continuous version given by the Fokker-Planck equation [23]: with the differential operator Here, x = (s, n) is a vector with continuous cell numbers as components, α = g S − g A − g 0 and β = g S + g A + g 0 .Regimes of growth: power-laws and avalanches Next, we want to examine the implications of the Fokker-Planck Eq. ( 2) for the lineage sizes within a tissue sample.We solve Eq. ( 2) with the initial condition f (x, t = 0) = δ (s − s 0 ) δ (n).s 0 corresponds to the initial number of stem cells of a lineage -not necessarily unity, since stem cells proliferate at initial stages of organoid preparation, which are not considered here otherwise.Using the Fourier transform of Eqs. ( 2), (3) with respect to x and the method of characteristics to solve the resulting first order partial differential equation (see supplementary material), we find the characteristic func- The distribution of total lineage sizes x tot is given by an integral of f (x, t) over all states with equal x tot : The presence of a critical point can be clearly seen from the expectation value of the stem cell population: for finite α, ⟨s⟩ = s 0 e αt , whereas for the critical α = 0, ⟨s⟩ = s 0 .For t → ∞, f tot (x tot , t) approaches a limiting distribution f t→∞ tot (x tot ) -a truncated 3/2-power-law Lévy distribution [24][25][26]: ( This formula holds for α ≪ β.For a more general expression see supplemetary Eq. (S 28).At criticality the distribution becomes a true 3/2-power-law.
The way in which f tot (x tot , t) approaches the limit of Eq. ( 5) is very different for α > 0 and α < 0. For positive α and large enough lineage size x tot > x * tot ∼ e αt , there is a region where f tot (x tot , t) is not approximated by Eq. ( 5).This is the avalanche region illustrated in Fig. 2 c).Here, lineages have a high percentage of proliferating S-cells which are driving the system's growth.The lineage sizes are very large but the probability of their occurrence is very small.For x tot < x * tot most lineages are fully differentiated and stopped growing.In the avalanche regime, f tot (x tot , t) can be approximated by where I 1 (z) is the modified Bessel function of the first kind and q * tot (t) and a (t) are defined in supplementary Eqs.(S 37) and (S 42).
For α < 0, we also find an avalanche region with an active S-cell population approximately described by Eq. ( 6) (even though the avalanche will stop eventually).As pointed out in the supplement, the approximation breaks down for αt > 1, but is still reasonable for the data at hand.This region is located at large x tot , for which f tot (x tot , t) > f t→∞ tot (x tot ), followed by a rapid truncation at even larger x tot (Fig. 2 b)).In contrast to the behavior at α > 0, f tot (x tot , t) convergences to f t→∞ tot (x tot ) uniformly: the avalanches becomes less and less pronounced as t → ∞, because the S-cells of all lineages eventually fully differentiate if α ≤ 0 [18].Our analysis of the data shows that α ≲ 0 holds in experiments (see Table I).
Dynamics and criticality in experiments The experimental data consists of lineage identifier counts for 15 organoids that were sequenced on days 16, 21, 25, 32 and 40 -three copies for each day [12].Accounting for statistical and readout errors, they can be related to lineage sizes [22].Each organoid consists of ∼ 10 4 lineages, while total cell counts evolve from ∼ 10 5 at day 16 to ∼ 10 6 at day 40.The organoids grow undisturbed from day 11 onward, when they are not subjected to intrusive procedures anymore.Four parameters are determined from experimental data: α, β, s 0 and N .The data at day 40 is very roughly distributed according to an x −3/2 power-law (see Fig. 2), indicating near critical growth with α ≪ β.
To determine the parameters precisely, we use the empirical characteristic function of the data Here tot are the experimentally determined lineage sizes.The characteristic function of our model f (q tot , t) (Eq.(S 33)) is then least squares fitted to fexp (q tot ).Moving Table I.Parameters of the SAN model estimated from experimental data with one standard deviation errors.α is the net growth rate of the stem cell population, and β characterizes the stochasticity of the system.s0 and N are the average initial number of stem cells and the average number of N-cells produced by each stem cell, respectively.
the fitting procedure to Fourier space has several advantages: fexp (q tot ) is less noisy than the empirical probability distribution and no smoothing procedures such as kernel density estimates need to be used.Approximations needed when transforming f (q tot , t) to real space are avoided.Model fitting via the characteristic function has been considered e.g. in Refs.[27,28].Our estimates for the model parameters are shown in Table I.In Fig. 2 a) we show a comparison between the probability densities of the SAN model f tot (x tot , t) and binned histograms of the data at days 16 and 40.We used 500 bins for the data at day 40, and 80 bins for day 16.In the supplement, we show that the agreement is independent on bin-size.Small lineages (< 20 cells) were excluded since these lineages, with a high probability, died out at preparatory stages.f tot (x tot , t) is found using a Fast Fourier Transform of f (q tot , t).
The estimates of Table I show that organoid growth is indeed a critical process with |α| ≪ β.It also shows that even at day 40 the process is far from its t → ∞ limit (αt ≈ 0.8) and the organoids maintain an avalanche-like population of S-cells at large x tot .
The SAN model is only a minimal model of biological reality and as such cannot account for the full spread of the experimental data.It does, for example, not account for the data at small x tot (see Fig. 2).Many small lineages have died out in the early stages of the organoid development and are obviously not described by SAN dynamics.We assumed that all lineages consist of s 0 stem cells at day 11.This is only an average.While for large lineages the initial number of stem cells is unimportant since the growth is stochastic, it matters for lineages that differentiated quickly.Other neglected aspects are the time dependence of rates, and the non-markovian nature of division and differentiation.Despite these caveats, the SAN model proves to be surprisingly robust and fits the experimental data well.In supplementary section D, we perform Kolmogorov-Smirnov (K-S) tests on different intervals for the data at day 40.K-S tests are a sensitive tool to test the power-law behavior of empirical data [29][30][31].We find that the test produces high p-values (up to p ≈ 0.8) in the region between x tot ≈ 200 and x tot ≈ 2000, where the power-law is most pronounced.
Extinction trajectories We now turn to the influence of criticality on the dynamics of single lineages.We focus  9)).For large extinction times τ , the trajectories for s0 = 1 and s0 = 9 become equal for large times: In the critical state the population dynamics is dominated by stochasticity, and hence obtains a universal form independent of the initial conditions.on α = 0 and restrict ourselves to S-cells as the only dynamical component.Neglecting the dynamics of the n variable, we drop all but the first two terms in L in Eq. ( 3).The reduced Fokker-Planck equation is equivalent to the stochastic differential equation w (t) is the standard Brownian motion.It is known that any s (t) described by Eq. ( 8), at some time τ , reaches s (τ ) = 0 [18], i.e. the S-cell population of the lineage goes extinct due to differentiation.Using the Onsager-Machlup formalism [32], we find the extinction trajectory It describes the most probable path between s (0) = s 0 and s (τ ) = 0 (see Fig. 3).Here, τ = 4s 0 / √ 3β .Details of the calculation are given in the supplementary material.Eq. ( 9) and Fig 3 show that at criticality, the cells loose memory of the initial condition s 0 , since for This is not the case for α ≷ 0. Most importantly, since the lineage sizes are power-law distributed, an overwhelming majority of organoid cells will belong to a few very large lineages.For these lineages, τ ≫ τ will hold, meaning that their development will be largely independent of the initial condition s 0 (Eq.( 10)).This is an essential feature of critical growth: it is independent of the initial conditions.
Discussion Finally, we discuss possible implications of critical tissue growth.Tuning itself to the critical point, the organoid maximizes the stochasticity of the growth process.We are not dealing with stochastic growth at a given rate.Instead, the characteristic time scale of growth vanishes (α = 0) and growth is fully determined by the stochasticity of the process.This is in contrast to many other examples of organ growth and regeneration [33][34][35] where the growth process is arrested in a coordinated manner.In critical growth, the long term dynamics of large lineages, which make up most of the organoid, does not depend on initial conditions.Stochastic fluctuations have erased all memory of the lineage's initial size.This might hint at a biological advantage of the critical regime: the outcome of the growth process is less influenced by perturbations in its initial stages, when the tissue is most susceptible to disturbances.
Second, the critical nature of the process implies a mechanism balancing the rates of division and differentiation of stem cells.Such balancing mechanisms are known from homeostatic stem cell renewal [36].One distinguishes between mechanisms based on asymmetric division (only one daughter cell is a stem cell while the other differentiates) and population asymmetry (the rates of stem cell loss and division are balanced on population level).The first strategy cannot produce critical lineage dynamics because it is strictly deterministic.On the population level, stem cell niches are a well known regulatory mechanism for homeostatic tissues [7,37,38] found in intestinal crypts [37] and the adult human brain [39].The available niche space limits the number of possible divisions.For growing tissues, the competition for a scarce resource, e.g.space or nutrients, can provide a feedback loop that limits the stem cells' ability for division at the population level.Such a feedback loop is often encountered in models of self organized criticality (SOC) [19], where the event probabilities are balanced and tuned to the critical point e.g. by energy conservation.For cerebral tissue, the currently available data gives no clues to the nature of the balancing mechanism.The search, however, could inspire future experimental research.Recent advances in lineage tracing allow to study the spatial distribution of lineages in cerebral organoids using light-sheet microscopy and spatial transcriptome sequencing [40].Inferring the spatial dynamics of the growth process could help to identify the the mechanism behind organoid SOC.Quantitative lineage tracing experiments with other organoid types such as intestinal [41,42], retinal [43,44] or cardiac organoids [45,46] could reveal whether critical growth is specific to cerebral tissue, or whether it is a more general organizing principle.

SUPPLEMENTARY MATERIAL A. Solving the Fokker-Planck equation
All supplementary equations are marked with an "S", while equation numbers without "S" refer to the main text.The idea that a laplace transform of the S-cell (branching) part of the Fokker-Planck equation ( 2) yields a solvable first order partial differential equation was used by Feller in Ref. [16].Here we use the Fourier transform to obtain the distributions of total population sizes of the branching based SAN-model near criticality.
Although cell counts are discrete, we chose a continuous model to describe the data.For small s the results of the continuous and discrete approaches will differ, yet both are reasonable approximations of biological reality.While using the discrete version may seem preferential at first, cell division is only the end result of a series of changes to a cell's internal state as it transitions through the G1, S, G2 and M stages of the cell cycle.Because the continuous process accounts for these changes by allowing cell counts to change gradually, it may be in fact the model that is closer to reality.

The characteristic function
The Fokker-Planck equation ( 2) can be solved in the Fourier domain.After a Fourier transform, Eq. ( 2) becomes where q = (s, n).In the following we will drop the factors of N for notational simplicity and restore them later on with the substitution Let us choose the initial condition which corresponds to s 0 S-cells at t = 0.In Fourier space this initial condition translates to Eq. (S 1) contains only first order derivatives and can be solved with the method of characteristics.The characteristic equations are The first equation is trivial.It is solved by t = τ .The second equation is solved by The solution for the Fourier transformed Fokker-Planck equation (S 1) is found by solving for the integration constant C 1 : qa) .(S 10) The next step is to look for a function of C 1 , G (C 1 (q s , q n , t)), which satisfies so that a function f (q, t) that satisfies the initial condition of Eq. (S 5) can be written as Such a function is given by G (q s , q n , t) = q s ρ (q n , t) − 2 β λ (q n , t) D (q a ) 1 + iλ (q n , t) q s (S 13) The function f (q, t) of Eq. (S 12) is the characteristic function of f (x, t).For later purposes we want to separate κ (q n ) into real and imaginary parts.We find ) ) 2. Large t limit.Distribution of N cells.
Let us gain some intuition into the dynamics of the growth process.If the S-cells are dividing at a near critical near zero rate α, while they are differentiating at a much larger rate g A , we can expect that, as t grows, most lineages will consist of differentiated cells.Lineages that escaped differentiation will have to be very lucky and sufficiently large.As a first step, we want to consider the distribution differentiated cells.
The inverse Fourier transform of Eq. (S 12) with respect to q s reads f (s, q a , t) = dq s 2π e iqss−is0G(qs,qn,t) .( This expression can be rewritten as We have thus separated the function f (x s, q n , t) into a part which depends on x s and describes lineages that are still evolving and a second part with x s = 0.This latter part contains information on the distribution of N-cells in fully differentiated lineages.Its inverse Fourier transform with respect to q n is given by Supplementary Figure S 1.The branch-cut and contour integration in Eq. (S 23).
Notice that f (s = 0, n) does not depend on time.We will later see that f (s = 0, n) is the t → ∞ limit of the distribution of total lineage sizes (see Eq. ( 4) in the main text).Note that the integrand of the first right hand side term in the second line of Eq. (S 17) approaches zero as q s → ∞.This is a consequence of separating out the delta function and is necessary to regularize the integral.We will attend to this issue below.To do the integral in Eq. (S 18), it is useful to take a look at the analytic structure of the integrand.Nonanalyticities arise from the square root structure of κ (q n ).Two branch cuts start at the two imaginary roots of κ (q n ) and run to +i∞ and −i∞ along the imaginary q n axis (see Fig. S 1).The positive imaginary root is while for the negative imaginary root we find For small |α| ≪ β, we have We thus make a crucial observation: The branch-cut in the upper complex half-plane of q n descends to the origin for α = 0, while the branch cut in the lower half plane always starts below −2iβ/ (β − g A ).For n > 0, the integral over the real line in the inverse Fourier transform in Eq. (S 18) can be deformed to an integral around the branch cut running from +i∞ to q * a to the left hand side of the branch cut (at a small distance ε, say), and then running back to infinity on the right hand side (see Fig. S 1).The half-circle around q * a is of order ε and can be neglected.We can write On the new contour, we need It is clear from Eq. (S 23) that the most important contribution to the integrals will come from the vicinity of q * n , because the integrand decays exponentially for larger q n the faster, the larger n.Expanding κ (iq n ) around q * n , we obtain (approaching the branch cut from different sides changes the sign of the square root).Thus, for large n, we approximate Eq. (S 23) as For the remaining integral we obtain This is a truncated, one-sided Lévy distribution [24][25][26].The appearance of the Lévy distribution could have been anticipated, since the integrand of Eq. (S 18) for α = 0 can be mapped to the characteristic function of the Lévy distribution after the argument of the square root in κ (q n ) has been expanded for small q n .For small α and g 0 we find g A = (−α + β − g 0 ) /2 ≈ β/2, as well as p ≈ β 2 /2, and for n ≫ 1 the truncated power law of Eq. (S 28) becomes Eq. ( 5) of the main text with N = 1.In the main text, we stated that the t → ∞ limit of f tot (x tot , t) -the distribution of total lineage sizes x tot -is governed by the same expression as f (s = 0, n) (see Eq. ( 5) of the main text).We will justify this statement below.On an intuitive level this can be understood as a consequence of the fact that, as t increases, most cells will already have differentiated and f tot will be dominated by the N-cell distribution.Since q * n descends to the origin as α → 0, the truncation of the 3/2-power-law of (S 28) vanishes and (S 28) becomes a true Lévy stable power-law distribution.Returning to the observation that the integrand of (S 18) is nonanalytic in the lower complex half plane, we conclude that f (s = 0, n) is finite for n < 0, meaning that there is a finite probability to find a negative number of N-cells.This is a consequence of approximating the discrete master equation by a continuous Fokker-Planck equation which is accurate for large s and n.However, one can easily convince oneself that f (s = 0, n) decays very fast as n decreases below zero: Following the reasoning of the above branch-cut integration, we see that f (s = 0, n) will be suppressed by an exponential factor e −q ′ n n .On the other hand, we see from Eq. (S 22) that for small α q holds.Thus we conclude that the probability for a negative n is very small, and can be neglected as an artifact of the Fokker-planck approximation.Notice that the probability for negative s is zero, because Eq. (S 12) is analytic in q s , except for a single pole which is always in the upper complex half plane (since Re (λ (q n ) > 0) as Eqs.(S 14), (S 13) indicate).Indeed, while the diffusion coefficient associated with the second derivative in s in Eq. (3) vanishes at s = 0, the diffusion coefficient associated with n-diffusion does not vanish at the origin and allows for a small leakage of probability towards negative n.So far we have been investigating the analytic structure of the s = 0 contribution to f (q, t).However similar arguments carry over to the general case.Multiplying the numerator and denominator of the function G (q s , q n , t) in Eq. (S 13) by λ (q n , t) −1 and keeping q s real for the moment, we see, that the characteristic function f (q, t) (Eq.(S 12)) is analytic in q n except for the branch-cuts of κ (q n ) and the pole at λ (q n , t) = i/q s .This pole occurs when λ (q n , t) is purely imaginary.f (s, n, t) will decay for n < 0 the faster, the lower the nonanalyticities with Im (q n ) < 0 lie in the complex plane.Eq. (S 29) shows that the contribution of the lower half branch-cut is indeed negligible since it starts much further away from the origin than the upper branch-cut.A numerical inspection of λ (q n , t) shows that regions where λ (q n , t) is purely imaginary in the lower complex half plane of q n are indeed sufficiently far from the origin for all reasonable parameter choices.We conclude that the finite values of f (s, n, t) for negative n are simply an artifact of the Fokker-Planck approximation to the original master equation ( 1) and are small.

The lineage size distribution ftot (xtot, t)
As pointed out in the main text, we are ultimately interested in the distribution of the lineage size This distribution is obtained by summing f (s, n, t) over all states which satisfy the condition (S 30): (see also Eq. ( 4) in the main text).Having argued in Sec.Supplementary Sec.A 2, that the distribution f (s, n, t) corresponding to the characteristic function (S 12) is confined to positive s and -to a good approximation -to positive n, we can extend the integration in Eq. (S 31) to the complete real axis: A Fourier transform then gives f (q tot , t) ≈ f (q s = q tot , q n = q tot , t) , ( where the characteristic function of Eq. (S 12) is on the right hand side.This yields Analytic expressions for the power-law and avalanche regimes The characteristic function given in Eq. (S 33) has a peculiar behavior at q tot = 0 which determines the asymptotics for large x tot .To see this, let us first investigate the t → ∞ and the q tot → ∞ limits of f (q tot , t).Since Reκ (q tot ) > 0 holds for all real q tot (see Eq. (S 14)) and Reκ (q tot ) ∼ |q tot | for large arguments, the approximation holds in both, the t → ∞ and the q tot → ∞ limits.Using (S 34), Eq. (S 33) can be approximated by For α > 0, this expression approaches exp (−2αs 0 /β) as q tot → 0, whereas for α ≤ 0, it approaches unity.For the full characteristic function, however, f (q tot , t) = 1 is true in all cases.This is depicted in Fig. S 2. We conclude that Eq. (S 35) is not a good approximation for small |q tot | if α > 0 holds.To find an approximation for small |q tot | Supplementary Figure S 2. The full characteristic function f (qtot, t) of Eq. (S 33) and the approximations for t → ∞ or large qtot, f (qtot, t → ∞), (Eq.(S 35)), as well as for small qtot, f (qtot → 0, t), (Eq.S 36).We used a positive growth rate α = 0.2/day and β = 10/day.For α > 0 and large but finite t, f (qtot, t → ∞) (green dashed curve) is a good approximation to the characteristic function f (qtot, t) (blue curve), except for a small region near qtot = 0 which is well approximated by f (qtot → 0, t) (orange curve).This region governs the behavior at large lineage sizes xtot.
and positive growth rates, we expand the numerator and denominator inside the exponential function in Eq. (S 33) around q tot ≈ 0 and find where In the last line we assumed that |α| ≪ β.Therefore, the appropriate approximation for α > 0 and small |q tot | reads It remains to determine the distribution function f (x tot , t).At small and intermediate x tot , we expect f (x tot , t) to be governed by the behavior of the characteristic function at larger q tot , whereas the behavior of f (x tot , t) for large x tot will be determined by the region around q tot ≈ 0. This has to do with the analytic structure of Eq. (S 33).The question is, which nonanalyticity dominates the Fourier transform as the integration contour is deformed according to Fig. S 1.Besides the branch-cuts of κ (q n ), f (q tot , t) exhibits a pole on the imaginary axis (see Fig. S 3).Let this pole be located at q tot = iq * tot .If it is closer to the origin than the upper starting point of the branch-cut, i.e. q * tot < q * n , where q * n is defined in Eq. (S 19), the pole becomes dominant at large x tot .This is because both nonanalyticities, the branch cut starting at q * n and the pole at q * tot , are located on the positive imaginary axis.Hence, when the integration along the deformed contour is performed, the contributions of the pole and the branch-cut will roughly behave as e −q * tot xtot , or e −q * n xtot , respectively.We can use Eq.(S 36) to track the behavior of q * tot for different parameter values: Eq. (S 37) shows that the analytic structure of the characteristic function is very different for positive and negative α.For positive α, it is q * tot < q * n , and therefore f (x tot , t) will be dominated by the single pole that is captured in Eq. (S 37).For α < 0, the opposite scenario is true: we find q * n < q * tot .
Supplementary Figure S 3. Branch-cut integration (see Fig. S 1 and Eq.(S 23)) and the pole at q * tot (Eq.(S 37)) characterizing the avalanche part of the distribution.Both, the pole and the branch-cut must be captured by the deformed contour (blue dashed line).For α > 0, the pole at q * tot is always nearer to the origin than the starting point of the branch-cut q * n and thus determines the large xtot behavior.For α < 0, q * tot is within the integration contour, moving towards larger imaginary values as t is increased.For sufficiently small t, the avanlanche part of the distribution is still well approximated by the integration around the q * tot pole (see Fig. 2 b)).
We now calculate the inverse Fourier transforms Eqs.(S 35) and (S 36) which will give us approximate expressions for f (x tot , t) at intermediate and large x tot .The inverse Fourier transform of Eq. (S 35) reads We encountered this integral when we calculated the distribution of N-cells in Eq. (S 18).Similarly to Eq. (S 28) we find In order to account for the N N-cells that each S-cell is producing we have to make the substitution q n → N q n in Eq. (S 32).This substitution carries over to Eq. (S 35) where we have to replace q tot by N q tot .The integration of Eq. (S 38) is evaluated according to For |α| ≪ β, we obtain the result of Eq. ( 5), Sec. of the main text: (S 41) According to the above discussion, Eq. (S 41) is valid for intermediate x tot if α > 0. Forα < 0, the contribution of the pole becomes stronger and stronger sub-leading to the branch-cut contribution as t is increased (see Eq. (S 37)).However, if t is sufficiently small, the pole is in the vicinity of the branch-cut starting point q * n and has a strong influence on the distribution function.For large enough t however, Eq. (S 41) is a good approximation for f tot (x tot ) everywhere.
To determine the behavior at large x tot for α > 0 we make use of Eq. (S 36) which gives a good approximation for the characteristic function at small q tot .It is useful to rewrite Eq. (S 41) as The inverse Fourier transform is strictly speaking divergent (as is the one in Eq. (S 17)).However it can be regularized by adding a small exponential factor e −qtot|ε| to the integrand and letting ε → 0 after the integral is done.The large x tot behavior that we are interested in is dominated by the pole at q * tot for α > 0. We will see later, that for α < 0, the pole still dominates the distribution for reasonably large x tot (avalanche regime) if t is small, although it does not govern the asymptotics for x tot → ∞ .The pole's contribution can be found by integrating along a circle around q * tot (see Fig. S 3).Writing q tot = iq tot + re iφ the integral around the pole at q tot becomes Since the radius r is arbitrary, we are free to chose r = a q 2 tot x tot .
After some algebra, we arrive at Using the integral representation of Bessel functions where C is a contour enclosing the origin, with the substitution t = e iφ Eq. (S 44) becomes For α > 0, this is the large x tot approximation for the distribution function f (x tot , t) given in Eq. ( 6), where we have restored the variable N .As is demonstrated in Fig. (2), Eq. (S 45) provides a good approximation for the large x tot behavior of f tot (x tot , t) for a positive S-cell growth rate α.However, even for α < 0, at reasonably large x tot (namely in the avalanche regime), the behavior of the lineage size probability density is well described by Eq. (S 45), if αt ≲ 1 holds.This is due to the fact that the pole at q * tot , while located above q * n on the imaginary q tot line (see Fig. S 3), is still sufficiently near q * n to dominate the avalanche part of the distribution function (see Eq. (S 37)).However, the asymptotics of f tot (x tot , t) for x tot → ∞ will not be given by Eq. (S 45).Since in both cases, for positive and negative α, Eq. (S 45) is an approximation for the avalanche part of the lineage size probability density, we call it f av.tot (x tot , t).

B. Most-likely paths
We now derive the most likely path s τ (t) that a critical S-cell population (i.e.α = 0) takes from s 0 cells at t = 0 to extinction at t = τ .Since we consider only S-cells and only critical populations, equation (3), accorting to the rules of Ito's calculus, reduces a process described by the stochastic differential equation (SDE) where w(t) is the standard Brownian motion.If the set of possible paths was finite-dimensional, we could proceed by finding the density W defined by equation (S 46) on the set of possible paths starting at s 0 and maximizing W (s τ ) subject to s τ (τ ) = 0 to find s τ .It turns out, however, that this approach only cleanly generalizes to infinite-dimensional path spaces in the case of a constant diffusion term [32].To avoid these technical difficulties, we perform a change of variable to transform equation (S 46) into a process with a constant diffusion term.By setting y(t) = 2 s (t) /β and applying Itô's lemma we get dy = A(y)dt + dw(t) where A(y)=− 1 2y (S 47) (were once y reaches zero, we define it to remain there, despite A becoming undefined).For this process, the density functional W expressed in terms of its Lagrangian L (also called Onsager-Machlup function) is [32] W [y(t)] ∼ e −S(y) (S 48) where The functional W defines a probability density on the set of paths in the following sense: The probability for a random path to lie within a small tube of diameter ϵ around a given differentiable path y is asymptotically proportional to W [y]H(ϵ) for some function H independent of y (this factorization fails if the diffusion term is not constant [32]).We can therefore proceed as we would in the finite-dimensional case and maximize W to find y τ .Since maximizing W means minimizing the action S(y) = 1 (where we chose the solution that ensures λt 2 + µt+ν >= 0 on [0, τ ]).Inserting these parameters into the general solution (S 49) and transforming back from y to s produces after some rearrangements the extinction trajectory stated in equation (9).
C. Full Fokker-Planck equation and the simplification of the A→A+N process The full Fokker-Planck equation for the SAN process of Fig where x = (s, a, n) is the vector of S, A, and N-cell counts.The simplification of the A→A+N process used in our calculations is that each A-cell, on average, produces N N-cells.This, essentially, amounts to cutting the last line of the above equation and rescaling the the a variable by a factor of N , giving Eq. (2) of the main text.Our justification for this simplification is that, while fluctuations in the number of S-cells s scale as ∼ √ s (see also Eq. (S 46)), the fluctuations of the N-cell count are -for a given a -Gaussian.This can be estimated from the diffusion term in the last line.While the large ∼ √ s fluctuations of the stem-cell number ultimately lead to the power-law at criticality, the Gaussian fluctuations of the number of neurons lead only to a broadening of the distributions.Llorca et.al. -eLife, 2019 (see main text for the full reference) considered the neuronal output of precurser cells (radial glia cells, which correspond to our A cells), and found that its uncertainty is only about an order of magnitude in size.As this is uncertainty is small compared to the uncertainty of the power-law distribution of total cell counts, we choose the simplified model, which is analytically tractable.

Figure 1 .
Figure 1.The SAN model of dividing and differentiating cells.Stem cells (green, S) either divide at rate gS (S→2S), differentiate at rate gA (S→A) or die at a rate g0 (S→Ø).Differentiated cells (orange, A) that committed to a developmental trajectory either divide asymmetrically (A→A+N) at a rate gA, or differentiate directly into non-dividing N cells (rate gF ).On average, each A cell produces N N-cells until it looses the ability to divide.At criticality, the rates gS (division) and gA (differentiation) are equal.

Figure 3 .
Figure 3. Single lineage stem cell populations for initial sizes s0 = 1 (solid lines) and s0 = 9 (dotted lines) and different extinction times τ (see Eq. (9)).For large extinction times τ , the trajectories for s0 = 1 and s0 = 9 become equal for large times: In the critical state the population dynamics is dominated by stochasticity, and hence obtains a universal form independent of the initial conditions.