Statistical theory of branching morphogenesis

Branching morphogenesis remains a subject of abiding interest. Although much is known about the gene regulatory programs and signaling pathways that operate at the cellular scale, it has remained unclear how the macroscopic features of branched organs, including their size, network topology and spatial patterning, are encoded. Lately, it has been proposed that, these features can be explained quantitatively in several organs within a single unifying framework. Based on large‐scale organ reconstructions and cell lineage tracing, it has been argued that morphogenesis follows from the collective dynamics of sublineage‐restricted self‐renewing progenitor cells, localized at ductal tips, that act cooperatively to drive a serial process of ductal elongation and stochastic tip bifurcation. By correlating differentiation or cell cycle exit with proximity to maturing ducts, this dynamic results in the specification of a complex network of defined density and statistical organization. These results suggest that, for several mammalian tissues, branched epithelial structures develop as a self‐organized process, reliant upon a strikingly simple, but generic, set of local rules, without recourse to a rigid and deterministic sequence of genetically programmed events. Here, we review the basis of these findings and discuss their implications.


| INTRODUC TI ON
To sustain life, organisms must exchange nutrients and metabolic waste products with the environment. In unicellular organisms, such as bacteria or yeast, where the surface to volume ratio is high, the surface of the cell is large enough to meet these demands. However, in multicellular organisms, such as mammals, strategies must be developed to maximize the area of surfaces where such exchange can occur. In the small intestine, this challenge is met by organizing the epithelium into an array of finger-like protrusions, known as villi, which extend into the gut lumen. In volumnar tissues, such as kidney, lung, mammary gland, prostate and pancreas, exchange surfaces are packed efficiently around ramified branched epithelial networks.
How do these structures form? How do instructions encoded at the molecular and cellular scale translate into the large-scale organization of complex branched epithelia? This is the problem of branching morphogenesis (Iber & Menshykau, 2013) and is exemplified by the pubertal development of the mouse mammary gland epithelium (Sternlicht, 2005).
In mouse, the mammary glands are specified along the ventral epidermis around embryonic day (E)12 as placode-like structures that sprout and invade an adipocyte-rich stroma. At birth, the mammary gland comprises a small rudimentary tree-like structure involving a minimally branched network ( Figure 1a). Then, during puberty, cellular precursors -termed "mammary stem cells" -drive a serial process of ductal bifurcation and elongation, leading to the specification of a complex ramified ductal network that extends to fill the fat pad (Figure 1a). In adult, hormonal changes through the estrous cycle promote bouts of alveoli growth and regression while, in pregnancy, alveoli mature into the milk-producing glands. The ducts form a simple stratified epithelium comprised of an outer layer of myoepithelial basal cells and an inner layer of luminal cells (Figure 1a).
During puberty, measurements based on short-term incorporation of thymidine analogues shows that the ductal growth is driven by actively cycling progenitors positioned at or near ductal tips -known as "terminal end-buds" -localizing the mammary stem cell population to these sites.
What is the molecular identity, sublineage potential, and fate behavior of mammary stem cells during pubertal development?
Are they stem cells at all? How do stem cells and their progeny integrate fate choice with collective cell rearrangements to direct the large-scale patterning of the ductal network? And are these mechanisms conserved in the patterning of other branched epithelia?
Traditionally, to address the mechanisms that regulate mammary gland development, most studies focus on the repertoire of transcription factors and signaling pathways that regulate cell fate behavior in the terminal end-buds (Macias & Hinck, 2012). But, to address factors that regulate the spatio-temporal patterning and large-scale organization of tissue, these may not be the most useful starting variables. Instead, to resolve the factors that control collective cell fate behavior and patterning, it makes sense to start by considering the larger-scale structural organization of the complex ductal network. Recently, by combining lineage tracing strategies with morphometric measurements of the ductal network structure, recent studies have provided evidence of a conserved mechanism of branching morphogenesis in the mouse mammary gland Scheele et al., 2017), kidney  and pancreas (Sznurkowska et al., 2018). Here, we review the basis of these findings and discuss their wider implications.

| THE LARGE-SCALE ORGANIZATION OF THE MAMMARY GLAND DUCTAL NETWORK IS PREDICTED BY A SIMPLE STATISTICAL RULE
First, to define quantitatively the large-scale structure of the mouse mammary gland epithelium, the ductal organization was traced from whole-gland reconstructions of tissue acquired at the end of puberty and stained for the ductal basal cell marker Keratin 14 ( Figure 1b). The results underline a remarkably complex arrangement, with ductal networks adopting a non-stereotypic organization (Lu, Sternlicht, & Werb, 2006), foliating into a diversity of subtrees of variable size and topology: After several rounds of near-symmetrical dichotomous branching, the resulting subtrees were found to be highly variable, with some subtrees terminating after just two or three further rounds of branching while others extended over 20-30 rounds (Figure 1c).
Combining the results of EdU incorporation, as a marker of proliferation, and whole-mount imaging of the whole mammary gland, the relative abundance of "active" terminal end-buds was found to steadily diminish during puberty , suggesting that terminal end-buds progressively and collectively exit cell cycle during this phase. But what underpins such network heterogeneity? Does the complexity arise from the early specification of mammary stem cells with variable proliferative potential, or do mechanical, chemical or other environmental cues influence distinct fate decisions of equipotent mammary stem cell pools? To discriminate between these possibilities, evidence was sought for changes in the potency of terminal endbuds during pubertal growth. However, notably, after the initial specification of the rudimentary ductal tree, over the remaining course of pubertal development, the average length and width of ductal segments remained approximately constant as a function of branching index -the latter defined as the minimum number of branches between a given ductal segment and the origin of the ductal tree. Moreover, the proliferative activity of terminal endbuds that remain in cycle, as assayed by the constituent fraction of EdU+ cells in the end-bud, also remained approximately constant over the developmental time course .
Together, these results suggested that the potency and proliferative activity of cycling mammary stem cells remains largely unchanged during the phase of pubertal growth. So, if terminal end-buds, and their constituent stem and progenitor cells, remain equipotent during puberty, what is the source of ductal network heterogeneity?
During puberty, the localization of cell proliferation to terminal end-bud regions allied with the network topology suggests that the choice between terminal end-bud bifurcation (in which the number of active mammary stem cells is doubled) and "termination" (in which mammary stem cells and their immediate progeny collectively exit cell cycle) is not predetermined in an intrinsic, deterministic fashion, but is made stochastically. To test this conjecture, probabilities can be assigned to these events as a function of branch index, i.e. at a given generation, with probability q, a terminal end-bud becomes inactive (with all cells exited from cell cycle) while, with probability 1-q, an active terminal end-bud undergoes a symmetrical bifurcation, replicating the size, potency and activity of cells in the newly formed end-buds ( Figure 1d). Then, empirically, from the statistical ensemble of mammary gland networks, the variation of the average probability q with branch index was determined. In itself, such an assignment does not challenge a "model" of stochastic growth. To determine whether such a dependence provides the statistical basis of the network organization, one must further check that networks generated from this statistical rule faithfully predict the statistical organization of the reconstructed glands and their constituent subtrees.
Implementing this program, it was found that the termination probability rose monotonically from zero at the lowest branch levels becoming saturated at around q ≈ 1/2 for the later phase of pubertal development ( Figure 1d). In other words, over much of the developmental time window, active terminal end-buds appear to evolve according to a simple statistical rule in which, with approximately equal probability, they duplicate through ductal bifurcation, or terminate through collective cell cycle exit. But can this dependence alone predict the complexity of the mammary ductal network?
Notably, within this framework, an estimate of the frequency of terminal end-buds that remain in cycle agrees well with that measured F I G U R E 1 Embryonic development of mouse mammary gland epithelium. (a) At birth, the mouse mammary gland epithelium forms a rudimentary ductal tree-like structure (upper-left panel). Expansion of the ductal epithelium is driven by proliferative cells at the ductal tips (marked in red), that drive a sequential process of ductal elongation and bifurcation (upper-right panel). As epithelial cells exit the ductal tip -known as the terminal end-bud -they exit cell cycle giving rise to a simple bilayer epithelium comprised of luminal cells and myoepithelial basal cells (bottom panel). (b) Outline of the ductal epithelial network of a mouse at the end of puberty (8 weeks), when it has expanded to fill a fat pad. (c) Topology of the corresponding ductal network showing that some subtrees terminate early while others go through multiple rounds of division. (d) Schematic (upper panel) shows that ductal subtrees can be segmented as a sequence of collective fate decisions in which active terminal end-buds choose stochastically between termination (cell cycle exit), with probability q, and bifurcation, with probability 1-q. Averaging over multiple terminal end-buds, the probability q is shown empirically to converge towards q = 1/2 (points). The line shows the result of a numerical simulation of the model discussed in the main text and Figure 2a. experimentally based on short-term EdU incorporation . Moreover, using this empirical statistical rule to estimate the subtree persistence -the chance that subtrees survive to a given branch index -and the subtree size distribution, the stochastic model provided a good quantitative prediction of the experimental measurements .
Together, these observations suggested that the complex network topology of the mouse mammary gland epithelium is defined by statistical rules that operate at the ductal scale. But what is happening at the cellular scale to affect the collective dynamics of actively proliferating terminal end-buds? What is the size and sublineage potential of the constituent mammary stem cell pool? And how is the near-balance between terminal end-bud bifurcation and termination controlled?

| UNB IA S ED CLONAL LINE AG E TR ACING RE VE AL S THE MULTIPLI CIT Y AND P OTEN C Y OF MAMMARY S TEM CELL S DURING PUB ERT Y
To address the cellular basis of mouse mammary gland development during puberty, it makes sense to deploy a genetic lineage tracing strategy to trace the long-term fate behavior of individual cells (Desai, Brownfield, & Krasnow, 2014;Van Keymeulen et al., 2011). To

| MAMMARY G L AND DUC TAL MORPHOG ENE S IS A S A B R AN CHING -ANNIHIL ATING R ANDOM WALK
To understand the mechanisms that regulate the collective cell dynamics of terminal end-buds, the large-scale spatial organization of the ducts was found to be revealing. Branching morphogenesis of mouse mammary gland takes place in a largely two-dimensional setting, where the frequency of ductal crossovers is low . This observation led to the conjecture that, during pubertal development, terminal end-buds may grow and branch at a constant rate, but terminate when exposed to secreted factors released from maturing ducts (Figure 2a) -a dynamic that, in the language of nonequilibrium statistical physics, constitutes a "branching-annihilating random walk". Since the trajectories of tips are taken to follow a persistent random walk, this provides a framework in which termination events are effectively stochastic at the network level, although they are tightly regulated at the single tip level. To this process, a further condition was imposed that required end-buds to terminate at the boundary of the fat pad. While the general topology of the resulting network mirrored qualitatively the observed structures (Figure 2a), the value of the model lay in its ability to predict quantitatively the statistical organization of the ductal network.
Remarkably, this model, which depends only on one (measurable) key parameter, the ratio of the branching to elongation rate, predicted the statistics of the branched ductal network, from the evolution of branch probability, q, with ductal index, to the distribution of subtree size and persistence (Figures 1d and 2b and Hannezo et al., 2017). Indeed, the stochastic character of the branching probability, a simplifying assumption within the framework of the model, was consistent with the measured distribution of branch lengths, which were seen to fall onto a strikingly exponential distribution . Further analytical insight can be gleaned from the development of a coarse-grained or hydrodynamic theory of network growth. By defining the density of active terminal endbuds and inactive immobile ducts, the branching dynamics takes the form of a reaction-diffusion system . In this framework, it becomes apparent that, during the growth phase, the system becomes self-organized into a soliton-like front -a "Fisher-KPP pulse" -of active terminal end-buds on the periphery of the growing network that travels in a directional manner at constant speed, leaving in their wake a constant density of inactive ducts.
Given this self-organizing behavior, the observed convergence towards balanced fate of ductal growth becomes easy to understand: With active terminal end-buds localized at the boundary, for each ductal bifurcation, on average, one daughter branch pioneers virgin territory while the other collides with the trailing network and terminates. Notably, as well as these qualitative dependences, from measurements of the growing ductal network, this minimal model could predict quantitatively the constant density of maturing ducts and pulse of active tips, including the "universal" ratio of the exponential decay constants between the advancing and trailing edge of the active pulse .
But what about statistical measures relating to fluctuations of the network? As a driven non-equilibrium system, statistical fluctuations of ductal density are expected to be large. In the parlance of statistical physics, these are known as "giant density fluctuations" and take the form of a power-law dependence of density fluctuations on the average ductal density with an exponent larger than one-half, the value expected for a purely random distribution.
Indeed, comparison of the data and the model revealed a strikingly power-law dependence with the same elevated exponent of approximately two-thirds.
Altogether, these findings suggested that the growth of the mammary gland ductal network is consistent with a minimal selforganizing principle based on local rules. By correlating cell cycle exit with exposure to local factors released from maturing ducts (including TGFβ), the sequential process of ductal branching and elongation leads to the specification of a ductal network of constant average density. Could such a model based on the branching-annihilating random walk be a general paradigm of ductal morphogenesis?

| DUC TAL MORPHOG ENE S IS OF THE MOUS E K IDNE Y AND PAN CRE A S
In mouse, the kidney develops as an outgrowth of the nephric duct that invades into the mesenchyme through a process of dichotomous ductal branching. Then, in a complex and cooperative process, the mesenchyme differentiates into epithelial tubes (nephron segments) that connect tips to form the basic filtration units of the kidney. In vitro studies involving the plating in two-dimensions of two embryonic kidney explants next to each other revealed an arrest of tip growth between tips and ducts, preventing collisions and ductal crossovers (Davies, Hohenstein, Chang, & Berry, 2014). The overall structure and dynamics of such two-dimensional kidney explants was thus accurately predicted by the same branching-annihilating random walk theory used for mammary glands . Interestingly, tip arrest in these explants was shown to be dependent on Bmp7, a member of the TGFβ super-family (Davies et al., 2014), which hints that the core findings from mammary gland could be translatable in other organs. Therefore, based on published data on the detailed structure of three-dimensional kidney morphogenesis in vivo (Sampogna, Schneider, & Al-Awqati, 2015), the question of whether there is a statistical basis to the branching network topology was addressed . Notably, in the three-dimensional system, as well as the measured ratio of ductal elongation to branching rates, a second parameter had to be considered -the contact distance within which active tips become terminated against maturing ducts.
Once again, comparison of the ductal network statistics obtained from measurements made during embryonic development with numerical simulations generated from the branching-annihilating random walk paradigm showed remarkably good agreement over a range of developmental time points . In particular, the data indicated that, although branching was seemingly stereotypic early in development, with most tips belonging to similar generations, this changed markedly post E15.5 (Sampogna et al., 2015), with considerable widening of the tip generation numbers, consistent with the predictions of the model. Quantitatively, by adjusting the two parameters of the theory -the ratio of the branching and elongation rates, and the contact radius of tip annihilation, good quantitative agreement was found for the variation of termination probability with branch index, as well as the distribution of branch number, subtree size and subtree persistence across a range of developmental time points. Indeed, these results suggest that the early phase of symmetric branching may not be a distinct phase of growth, as traditionally thought, but may be a natural consequence of the three-dimensional branching dynamics where the chance of meeting a neighbor early in development is diminished.
Interestingly, a concomitant analysis of the branching rules of kidney morphogenesis up to E15.5 (Lefevre et al., 2017) has proposed an alternative hypothesis. In this framework, it has been argued that ductal organization results from asymmetrical branching rules where, upon each branching event, one tip would branch again soon after, while the other, potentially because of its close contact to neighboring tips/ducts, would be delayed in its branching capacity. How such an intrinsic, contact-induced, delay process would result in stereotypically delaying one tip but not the other upon branching remains unclear. Unfortunately, since both models predict very few "tip terminations" pre-E15 (around 10% in the case of the branching-annihilating random walk), it is difficult to discriminate between these two competitive hypotheses using pre-E15 data alone. Frustratingly, around E15, when differences in model predictions would begin to impact, the two different experimental datasets yield different levels of branching heterogeneity, a discrepancy that will need to be resolved. However, branching data from later time points (E15.5-E19.5, Sampogna et al., 2015) reveals a dramatic widening of the subtree size distributions and tip generation numbers, which can be well-fit by the branching-annihilating random walk paradigm . This contrasts with a model based on deterministic tip delay, which would predict that all tips should have comparable generation numbers (within a small variance of three to four generations).
An outstanding question would thus be to understand the molecular basis for such late-stage tip generation heterogeneity. Although nephrogenesis had previously been proposed to suppress branching (Costantini & Kopan 2010;Sweeney, Lindström, & Davies, 2008), and could have been a candidate to track terminations , results from a recent study suggest that the genesis of nephrons does not correlate with delays in branching (Short et al., 2018). This argues that density-dependent feedback on branching/elongation would occur via other mechanisms, leaving open the question of the cellular and molecular nature of tip termination and/or delay in vivo. Analysis of larger data sets, time windows and mouse models will be required to resolve the basis of kidney branching morphogenesis.
In a parallel study, the basis of ductal morphogenesis of mouse pancreas was addressed using a similar approach to the mouse mammary gland (Sznurkowska et al., 2018). In adult, the mouse pancreas forms a quasi-two-dimensional structure in which acinar cells lie anchored to the tips of an intricate and complex ductal network, which is interspersed with islets of Langerhans (Shih, Wang, & Sander, 2013). Mouse pancreas development initiates around E8.5 and evolves through a complex two-stage process involving an early phase of plexus formation followed by an extended phase of plexus remodeling, which has been proposed to result from optimization of fluid-transport rules (Dahl-Jensen et al., 2018), and tip-driven branching morphogenesis from the pancreas periphery (Bankaitis, Bechard, & Wright, 2015 and Figure 2c). To study the cellular basis of pancreas development, a clonal cell lineage tracing assay was combined with morphometric measurements of the large-scale ductal organization.
Notably, measurements of the ductal subtree sizes revealed a distribution that, when rescaled by the average, overlapped closely with that obtained from studies of the mouse mammary gland (Sznurkowska et al., 2018). These findings suggest the basis of ductal morphogenesis in pancreas may involve the same branchingannihilating random walk paradigm. This conclusion was reinforced by unbiased lineage tracing using the R26-CreERT2;R26-Confetti reporter system, which revealed qualitatively the existence of "tree-shaped" acinar and ductal clones (i.e. clones closely tracking a single subtree in a monoclonal manner), and quantitatively a close correspondence between the statistical distribution clone and subtree sizes. This is expected from a branching process driven by selfrenewing "stem-like" cells localized at or near the tips of active ducts.
Although this analysis alone doesn't determine the exact transition point between central plexus remodeling and peripheral branching morphogenesis, a strength of the branching and annihilating framework is that it is self-organized, and thus insensitive to initial conditions such as the starting geometry of the remodeled plexus (Sznurkowska et al., 2018).
At the cellular level, pancreatic precursors become increasingly sublineage-restricted during embryonic development, with tripotent cells giving rise to self-renewing bipotent ductal-islet precursors and lineage-restricted acinar precursors which co-localize at the growing ductal tips, duplicating through serial rounds of ductal branching.
However, in contrast to mammary gland, inspection of clonal imprints on the maturing postnatal day (P)14 tissue suggests that the abundance of self-renewing cells is almost two orders of magnitude smaller than that found in the mammary gland. With confetti labelled ducts showing drift to monoclonality over just a few rounds of ductal branching, the size of self-renewing ductal population was estimated to be just a handful of cells with similar estimates for tiplocalized self-renewing acinar precursors.

| D ISCUSS I ON AND PER S PEC TIVE S
In summary, these studies suggest that, at the organ scale, the In an alternative strategy, a much more efficient pattern can be generated by adjusting the branch length continuously with branch index leading to a more regular fractal-like geometry with a high packing density. Studies of the ductal organization of mouse lung show evidence of a stereotypic pattern of early branching (Metzger, Klein, Martin, & Krasnow, 2008) based on a similar kind of organization, with geometric structures regulated by side-branching patterns leading to power-law scaling of ductal sizes with branch index (Horsfield, 1990).
Whether these stereotypic patterns and rules are conserved down to the finest branch scales, or whether the later phases of branching are governed by a simpler statistical paradigm remains in question.
A striking and conserved feature of the mammary gland and pancreas development is the apparent cooperativity of sublineagerestricted self-renewing cells. In the former, basal-and luminalrestricted mammary stem cells act cooperatively to generate the ductal epithelium, becoming proportionately expanded through ductal bifurcation. Similarly, in the pancreas, self-renewing ductalislet-restricted precursors act in concert with self-renewing acinarcommitted precursors at the ductal tips, duplicating through serial rounds of branching. How is this cooperativity enforced? One possibility is that, in common with the mouse adult trachea or developing lung, basal and luminal cells (viz. club cells in the trachea) act as a niche for each other, allowing the relative size of the compartments to achieve a stable equilibrium (Nabhan, Brownfield, Harbury, Krasnow, & Desai, 2018;Pardo-Saganta et al., 2015).
However, the problem of how sublineage-restricted stem cells control their ratio questions the identity of other niche factors that promote self-renewal activity. One possibility is that stromal cells surrounding the ductal tips -fibroblasts and endothelial cells -secrete factors that inhibit cell cycle exit (or differentiation). However, to ensure the long-term persistence of renewing cells, such factors would have to co-move with the elongating tip and be duplicated or recruited during ductal bifurcation -in essence, "pulling" the terminal end-bud or tip like the proverbial "carrot tied to the donkey".
Alternatively, the maturing ducts themselves may secrete a factor that drives cells to exit cycle (or differentiate) -"pushing" the terminal end-bud forward. Evidence in favor of key secreted factors from both the local tip environment (FGFs, etc.) as well as the maturing ducts (TGFβ, etc.) are present, suggesting that both mechanisms may act in concert.
Finally, to what extent does the branching-annihilating random walk model provide true mechanistic insight? The answer to this question is likely to be subjective. The branching-annihilating random walk model provides a predictive understanding of branching dynamics and so, for many must be considered as mechanistic.
Moreover, by surrendering information at the cellular and molecular scale, it affords a unifying or "universal" description of the largescale dynamics. Put differently, multiple viable mechanisms of molecular or cellular regulation could lead to the same dynamics at the ductal scale, belonging to the class of branching-annihilating random walks. However, for others, without a detailed understanding of the transcription factors and signaling molecules that effect the regulation of cell fate decisions, the model falls short -possibly well-short -of real mechanistic understanding. In the end, this debate is difficult to reconcile.
In general, mechanistic understanding of a given phenomenon must be tailored to the appropriate level of abstraction -physicists speak of "phenomenology". For example, we may understand the nature and action of a morphogenic program without resolving its microscopic basis. And such understanding may provide the means to frame targeted questions into the more microscopic level of description. In this sense, the problem of branching morphogenesis may serve as an exemplar for how to study and define the basis of emergent or collective behaviors in cell biological systems.