{'fraction_non_alphanumeric': 0.029268885110517368, 'fraction_numerical': 0.01462229778965266, 'mean_word_length': 4.890828444698813, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Kidney branching morphogenesis under the control of a ligand-receptor based Turing mechanism

The main signalling proteins that control early kidney branching have been defined. Yet the underlying mechanism is still elusive. We have previously shown that a Schnakenberg-type Turing mechanism can recapitulate the branching and protein expression patterns in wildtype and mutant lungs, but it is unclear whether this mechanism would extend to other branched organs that are regulated by other proteins. Here we show that the GDNF-RET regulatory interaction gives rise to a Schnakenberg-type Turing model that reproduces the observed budding of the ureteric bud from the Wolffian duct, its invasion into the mesenchyme, and the observed branching pattern. The model also recapitulates all relevant protein expression patterns in wild-type and mutant mice. The lung and kidney models are both based on a particular receptor- ligand interaction and require: (1) cooperative binding of ligand and receptor, (2) a lower diffusion coefficient for the receptor than for the ligand, and (3) an increase in the receptor concentration in response to receptor-ligand binding (by enhanced transcription, more recycling or similar). These conditions are met also by other receptor-ligand systems. We propose that ligand-receptor based Turing patterns represent a general mechanism to control branching morphogenesis and other developmental processes.


Introduction
Many organs, such as the lung, kidney, and glands are heavily branched structures. At the organ level, the branching pattern is defined by several parameters, such as the site and type of branching (terminal versus lateral branching, bifurcation versus trifurcation), branch angles, rates of elongation, and changes in tubular diameter. All of these processes are controlled to yield an overall pattern unique to each organ. The branched tree in the lung is generated by the sequential, non-random use of three geometrically simple modes of branching (lateral branching, planar and orthogonal bifurcation), which occur in defined routines [1]. Trifurcations have also been documented in the lung [2], but these are much more prevalent in the ureteric bud of the kidney. Culture experiments revealed that most branching events in the kidney are terminal bifurcations and to a lesser extent trifurcations, and only 6% of all branching events are lateral branching events [3][4][5].
It is a long-standing question how branching is controlled during development and whether the mechanism is the same in the different organs in spite of the differences in the architecture and in the regulatory proteins. We have recently proposed a model to explain the branching processes in the lung [6]. The model focused on the two key signalling factors that have been identified experimentally in the developing lung, fibroblast growth factor (FGF)10 and Sonic Hedgehog (SHH) as well as the SHH receptor PTCH1. We showed that the reported biochemical interactions give rise to Schnakenberg-type Turing patterns [7] that result in a distribution of FGF10 as would be expected for lateral branching and bifurcations. The model also reproduced all published mutant phenotypes, including the counter-intuitive widening of the clefts between buds as Fgf10 expression is reduced in an allelic sequence [8,9]. In further simulations we showed that the same core regulatory network together with FGF9 is capable of controlling the emergence of smooth muscles in the clefts between growing lung buds, and Vegfa expression in the distal sub-epithelial mesenchyme [10]. In the model the experimentally observed upregulation of Fgf10 expression by FGF9 [11] promoted lateral branching over the bifurcation mode of branching. Regulatory networks involving FGF10 and SHH also control the branching in the prostate [12,13], the salivary gland [14], and in the pancreas [15,16]. However, the regulatory proteins are different in the developing kidney [17,18], and the organs all differ in their branching pattern and geometry such that it is an open question whether the underlying regulatory principles are the same or different.
The ureteric bud forms at the beginning of embryonic day (E)10.5 at the level of the hindlimb in a specialized region called the metanephric mesenchyme [5]. The metanephric mesenchyme produces inductive signals that causes the nephric duct (also called Wolffian duct) to evaginate and to form a single ureteric bud near its caudal end [5]. After invading the metanephric mesenchyme, the ureteric bud undergoes about ten generations of repeated branching [19], followed by a period of elongation, then by one to two rounds of branching before birth [19]. Similar to the lung, the kidney collecting ducts form via branching and elongation of an epithelial cell layer [20]. Growth factors and components of the extracellular matrix are required, while the mesenchyme is dispensable [21]. The mesenchyme, however, influences the pattern of branching. One striking example came from tissue recombination experiments in which lung mesenchyme induced branching of the ureteric bud with a pattern characteristic of lung epithelium, i.e. with increased lateral branching [22].
Based on these observations the local expression of growth promotors and inhibitors has long been considered to determine the branching pattern. In particular, given the chemoattractive properties of the TFG-beta family protein Glial cell line-derived neurotrophic factor (GDNF) [23,24], it has been suggested that branching of the ureteric bud is caused by the attraction of tips toward local sources of GDNF [25]. GDNF signals to the cell via the receptors RET and GPI-anchored subunit GDNF family receptor alpha-1 (GFRα1). Ureteric bud outgrowth fails in Gdnf −/− and Ret −/− , Gfr α −/− mutant mice [5,20,[26][27][28]. Moreover, beads soaked with GDNF induce the outgrowth of extra ureteric buds in kidney culture explants [29]. Conditional knock-outs further revealed a positive feedback between GDNF/RET and WNT11 signalling as part of the core branching mechanism [20]. Expression of Wnt11 is activated in the epithelial tip of the ureteric bud and WNT11 signalling is in turn required to propagate mesenchymal GDNF signalling, which results in the establishment of an autoregulatory epithelial-mesenchymal feedback signalling loop [20]. FGF signalling, on the other hand, has only a supporting function in the ureteric bud. Thus, branching is still observed in Fgfr2 -mutant kidneys, though at a reduced level [30]. Moreover, enhanced FGF signalling in knock-out mice of the antagonist Sprouty can rescue Gdnf -/-and Ret-/-mutants, which otherwise fail to develop kidneys [31]. Interestingly, GDNF/RET and FGF10/FGFR have the same downstream target, the ETV4/ETV5 transcription factors [32], and it is therefore plausible that they may serve similar functions.
The role of bone morphogenetic proteins (BMPs) in kidney development is less clear. Appreciable levels of phosphorylated SMADs are first detectable at E12.5 in collecting duct trunks and loss of SMAD1 transcriptional activation of WIF1 is associated with reduced Wif1 expression and increased WNT/βcatenin signalling activity in lung epithelia, resulting in specific fetal lung abnormalities [33]. Addition of BMPs enhances branching at low concentrations and inhibits branching at high concentrations [34][35][36]. Equally, removal or misexpression of the BMP receptor Alk3 affects the branching pattern. Thus, expression of constitutively active Alk3 in collecting ducts using the HoxB7 enhancer-promoter reveals an inhibitory effect of BMP pathway activation on collecting ducts, confirming these in vitro findings [37]. Moreover, removal of the BMP antagonist Gremlin (Grem)1 in mice results in excessive BMP activity in the metanephric mesenchyme around the ureteric bud and disrupts the invasion of the mutant metanephric mesenchyme by the ureteric bud and the concurrent establishment of the autoregulatory GDNF/WNT11 feedback signalling loop [38]. In spite of these demonstrated effects of BMPs, removal of the common signal transducer Smad4 within this tissue has little impact on kidney development and branching morphogenesis up to E16.5. The signals of the TGFβ superfamily must therefore be transduced through SMAD4-independent pathways in the collecting duct [39]. To that end it was observed that BMP-2 antagonizes WNT signalling in osteoblast progenitors by promoting an interaction between SMAD1 and Dishevelled-1(DVL-1) that restricts β-catenin activation [40].
Given its detailed characterization and obvious differences to the lung we have focused on branching morphogenesis in the ureteric bud to define common principles in branching morphogenesis. While the signalling circuits in the developing kidney are still incompletely characterized, the available data is sufficient to build a first computational model of the regulatory interactions. We show in the following that the reported biochemical interactions between GDNF and its receptor can give rise to a similar Schnakenberg-type Turing mechanism as observed for the SHH-PTCH1 interaction in the lung. We show that the model is consistent with all available wildtype and mutant data and that it reproduces the main modes of branching, i.e. bifurcations and trifurcations, that are observed during kidney branching morphogenesis. Finally we show that the same regulatory network is capable of controlling the epithelial invasion of the metanephric mesenchyme. We argue that the coupling of a receptor-ligand based Turing patterning mechanism with the impact of the domain geometry may constitute a general mechanism to control branching morphogenesis.

A Computational Model for the Regulatory Network
Based on available data we develop and analyse a parsimonous 3-component model for branch point selection during kidney development. We focus on GDNF/RET signalling, as the key driver of branching morphogenesis, as well as on the GDNF/WNT positive feedback loop. The effect of BMP/GREMLIN1 signalling and of FGF/SPROUTY1 signalling will be discussed, but these components will not be explicitly included in the model. Accordingly the model is restricted to three proteins, GDNF, its receptors (RET and GFRα), and WNT11 ( Figure 1A). GDNF and WNT11 (which we denote by G and W ) are secreted proteins. The receptors RET and GFRα (which we denote by R) are membrane bound proteins and therefore diffuse on the surface of a cell with a diffusion coefficient much lower than that of secreted ligands: D R << D W , D G . Earlier studies have successfully described the in vivo distribution of morphogens with continuous reaction-diffusion equations on a domain with a length scale as small as 10 cells [41][42][43], and we therefore expect that the ligands in our model are also adequately described by continuous reaction-diffusion equations considering that the number of epithelial cells in any direction is larger than 10. The receptors are more of a concern as these are restricted to the cell surface and their diffusion is thus limited by the cell boundaries in the tissue. We have previously noted that during the receptor half life t 1/2 = ln(2)/δ R = 350 s receptors can diffuse over distances of less than the diameter of one epithelial cell l = (2D R t 1/2 ) 1/2 = 3.3 µm. Moreover, we previously showed numerically that qualitatively similar patterns are observed on a cellurarized lung bud domain [6]. In light of these previous results we also use continuous reaction-diffusion equations for the receptors.
Gdnf is expressed in the metanephric mesenchyme and signals to its receptor pair RET and coreceptor GFRα1 ( Figure 1A). The receptor Ret and co-receptor Gfr α1 of GDNF are both expressed in the epithelium [29]. The exact stoichiometry and kinetics of the process are still unknown but the GDNF dimer likely first binds to monomeric or dimeric GFRα1, and the GDNF-GFRα1 complex then interacts with two RET receptors. Some of the RET receptors may be weakly associated with GFRα1 before GDNF binding [44]. For simplicity in the first place we only consider one receptor R that represents the rate limiting species, and assume that each GDNF dimer is bound by two such receptor molecules. We note that the model would still work if the number of rate limiting receptor binding events was greater than two or if RET and GFRα1 were considered individually as discussed below. The unbound receptor is lost by degradation at rate d R [R] and by binding to the ligand GDNF at rate 2d C [R] 2 [G]. Receptors are constitutively expressed in the epithelium at rate p R and their expression is upregulated in the presence of GDNF [5,29] as described by the term p C [R 2 G]; we note that instead of the linear term we could also use a Hill function with a sufficiently large Hill constant. For simplicity we apply the quasi-steady state approximation to the concentration of the receptor-ligand complex, i.e.
where K C is the dissociation constant. We previously showed this approximation to be valid in our lung model for a broad range of physiological parameters [6]. We then have for the receptor dynamics: Gdnf is expressed in the mesenchyme at rate ρ G0 . Its expression is further enhanced by WNT11 signalling [17,20]. We describe this regulatory interaction by a Hill-type function, x m 1+x m , where m accounts for possible cooperative effects. We use m = 2 throughout, and we thus have p G [W ] 2 [W ] 2 +KW . GDNF binds to its epithelial receptor at rate d C [R] 2 [G] and is degraded at rate −d G [G]. These interactions result in the following mathematical formulation of the GDNF dynamics: Wnt is expressed constitutively in the epithelium at rate p W 0 and its expression is upregulated in response to GDNF-receptor complex formation and signalling, p W [17,20]. WNT is lost by degradation at rate d W [W]: In a separate model we also included the co-receptor GFRα1 explicitly in the model. Its dynamics can be described bẏ  Table 1 and a linear stability analysis confirms that these parameter values result in a Turing pattern. To observe the emergence of patterns, simulations were started with no species present, i.e. all concentrations were set to zero at time t = 0.
The idealized computational domain is in the shape of an outgrowing ureteric bud as shown in Figure  1B and comprises two tissue layers, epithelium and mesenchyme. Based on published images we determine the following geometric parameters at the time of branching: internal radius r 0 =50 µm, epithelial thickness l ep = 15 µm, radius of the internal cavity r 0 − l ep =35 µm, and radius of the metanephric mesenchyme r 1 =100 µm. The total domain length h 0 + r 1 varies and can reach a length of up to 200 µm [3,32,45,46]. Figures 2-4 were generated on a static domain with h 0 = 0.6, r 1 = 1. Figure 5 was simulated on a growing domain.The implementation of models on a growing domain is discussed in the Methods section.
There is no experimental indication that the mesenchyme and/or epithelium are surrounded by any insulating layer. We therefore assume that the secreted proteins are free to leave the secreting tissue by diffusion. To adequately represent this we embed the ureteric bud in a larger domain into which the ligands can diffuse. The size of this surrounding domain was chosen such that there are no boundary effects; details of the implementation are discussed in the Methods section. Because of the lower computational cost, most of the results presented in this manuscript are calculated with no-flux boundary conditions at the boundary of the domain representing the ureteric bud. However, all key results were also checked and reproduced with the ureteric bud embedded into a virtually infinitely large domain, into which the diffusible ligands can diffuse.
Finally, since the length and time scales of the process have been established we can convert the dimensionless diffusion and degradation parameters back to their dimensional counterparts and compare these to experimental values. These converted dimensional parameters are summarized in Table 2 and all lie well within the experimentally established physiological range. The protein concentrations have not been established.

The GDNF-Receptor Interaction Gives Rise to Modes of Branching Observed During Kidney Morphogenesis
Epithelial outgrowth is induced by GDNF signalling [5,20,[23][24][25][26][27][28][29]. The only known GDNF receptor which transmits the GDNF signal into the cell is RET, with its co-receptor GFRα1 [5,47]. We therefore use the GDNF-RET complex concentration, R 2 G, as a marker for the points of bud outgrowth. When we solve the model (equations 5) for the signalling interactions depicted in Figure 1A on the 3D computational domain depicted in Figure 1B we observe distributions of the GDNF-receptor complex, R 2 G, that can in principle explain the different branching modes observed in the kidney (Figure 2A-D, grey scale with white indicating the highest concentrations). Since quantitative data on absolute protein concentrations and expression levels are not available we have to restrict ourselves to a qualitative discussion of distribution patterns. We therefore deliberately left out scale bars for better readability. Depending on the choice of parameters we observe GDNF-receptor complexes concentrated either at the tip (elongation mode, Figure  2A), in two spots at the side (bifurcation mode, Figure 2B), in three spots on the side (trifurcation mode, Figure 2C), or at the tip and in spots on the side (lateral branching mode, Figure 2D). We note that small changes in the parameter values were sufficient to switch between the different branching modes, i.e. v = 2.5 for the trifurcation mode was changed to v = 1.5 to obtain the elongation mode in panel A, and to v = 2 to obtain the bifurcating pattern in panel B. The lateral branching pattern was obtained with a model that considered both the receptor RET and its co-receptor GFRα1 explicitly, and required as parameter values h 0 = 2, v 1 = 1.15, and ρ R1 = 1. The implications of such sensitivity on both, pattern robustness and regulatory potential of other factors, is discussed below. As part of the regulatory mechanism shown in Figure 1A signalling by the GDNF-RET receptor complex induces expression of Ret and Wnt11. As expected from the model formulation, the expression of these two genes in the bud domain is indeed positively correlated with the level of GDNF-RET signalling. Thus, when we measured the expression levels of Ret and Wnt11 expression as well as the concentration of the GDNF-RET receptor complex (R 2 G) at each mesh point we found that these are positively correlated ( Figure 2E). The co-localisation of the strongest levels of receptor and Wnt11 expression with the GDNF-receptor complex agrees well with the experimentally observed expression patterns [20]. WNT11 signalling in turn upregulates the expression of Gdnf in the mesenchyme. The Gdnf expressing zones are indeed adjacent to the epithelial GDNF-RET receptor signalling patches (Gdnf expression levels are shown in rainbow colour code in Figure 2A-D).

Branch Mode Preferences in Lung and Kidney
While the same types of branching events have been observed in lungs and kidneys their frequency of use differs greatly. Thus, lateral branching has been found mainly in the lung, but rarely in the kidney, while trifurcations are rare in the lung, but more prevalent in the kidney [1][2][3]. The dominant branching mode in the ureteric bud is terminal bifurcations [3]. Similarly, in our models we observe patterns that correspond to all observed modes of branching (Figure 2A-D), but their frequency differs. Thus, in our previous model for the lung we observed lateral branching events for a large part of the parameter space [6], while in the model for the ureteric bud lateral branches are observed only rarely, and if observed, they are very sensitive to changes in the parameter values ( Figure 3). Thus, adding as little as 5% Gaussian noise to the parameter values used for Figure 2D caused the pattern to change away from the lateral branching mode as shown for one example in Figure 3A, where the spot at the tip of the bud (that would support further bud elongation during lateral domain branching) vanished. On the contrary we could add 50% noise to the parameter values without affecting the trifurcation pattern ( Figure 3B). The kidney model would thus predict lateral branching events to be rare, while bifurcations and trifurcations are observed for a wide range of parameter values ( Figure 4A). While we cannot exclude that these pattern preferences are the result of our particular parameter choices in the two models for the lung and kidney, we note that the particular choice of parameter values is heavily constrained by the mutant phenotypes. The models for the lung and the kidney may appear mathematically similar at first sight but there are important differences. For one the FGF10/SHH interaction results in a negative feedback while the GDNF/WNT11 interaction results in a positive feedback. Furthermore, in case of the ureteric bud the outgrowth inducing ligand (GDNF) is produced in the mesenchyme and its receptor (RET/GFRα) in the epithelium, while in our lung model the outgrowth-inhibiting factor (SHH) is produced in the epithelium and its receptor (PTCH1) in the mesenchyme. However, the latter difference may be less significant as a similar Schnakenberg-type Turing model for the lung could also be constructed based on the FGF10receptor interaction. In that case, much as in the ureteric bud, the outgrowth-inducing factor, FGF10, would be produced in the mesenchyme. Such alternative model has not yet been analysed by us in detail and we can therefore not comment on the branch type preference.
To further establish the size of the parameter space for which we obtain bifurcations, trifurcations, and elongation patterns we varied each parameter one by one until we observed a mode change. As can be seen from Figure 4A trifurcation patterns are the most robust patterns. At the same time almost each parameter can be employed to switch the pattern between bifurcation and trifurcation modes of branching. Trifurcations and bifurcations thus appears to be robust to small variations in parameter values, yet sensitive to regulation. The many additional interactions that have been described, but that were not included in this parsimonious model thus can exhibit their documented impact on branching morphogenesis by altering one or several of these parameters. Thus, FGF signalling appears to control similar downstream targets as RET signalling [31] and Fgfr2 -mutant kidneys exhibit reduced branching as observed in simulations with less Ret expression ( Figure 4B). BMP signalling appears to antagonize WNT, and addition of high concentrations of BMPs inhibits branching and removal of the BMP antagonist Grem1 inhibits ureteric bud outgrowth from the Wolffian duct [33][34][35][36]38,40]. Likewise in the simulations lower levels of W production results in less branching ( Figure 4B). Interestingly, a model where RET and GFRα are considered individually is more robust to parameter variations than the simplified model with only one receptor, R ( Figure 4A). Thus, assembly of a complex from individually regulated subunits increases the robustness to variations in the concentrations of the subunits.

Mutants
An important test for the suitability of a mathematical model is its consistency with a wide range of independent experimental observations. Mutants represent a perturbation of the parameter set of the original system behavior and if a wide range of such perturbations can be reproduced correctly by the model then this provides strong support in favor of a model. Below we discuss the relevant mutants reported in the literature. To compare in silico mutants with those generated in vivo we used the following read-outs: the Gdnf expression level (Figure 4B), the GDNF-receptor signalling activity ( Figure  4C), and the extent of branching ( Figure 4D). GDNF is the main signalling protein controlling ureteric mesenchyme outgrowth and GDNF signals through its receptor RET and co-receptor GFRα. Accordingly, the GDNF-receptor activity was computed as the total concentration of R 2 G in the epithelium. GDNF signalling induces Wnt and Ret expression, and GDNF signalling levels therefore also serve as a proxy for the Wnt and Ret expression levels. The extent of branching was measured via the frequency of the Turing pattern which is inversely proportional to the distance between two spots. Smaller values of the frequency therefore correspond to kidneys with reduced branching and longer branches. Such correlation between the extent of branching and branch length is indeed observed in embryos. Thus, the Tgf β2 +/− mutant exhibits reduced branching and branches are of increased length [48].
Homozygous mutants were implemented by setting the protein production rates to zero (e.g. ρ W 0 and ρ W in case of Wnt11 ). Heterozygous mutants retain one active copy of the gene and accordingly the production rates were halved, even though we acknowledge that net changes may be different in mutants because of further feedbacks that were not considered in the model.

Gdnf Mutants
GDNF is an important ligand controlling kidney branching morphogenesis, and in the Gdnf −/− mutant the ureteric epithelium fails to invade the mesenchyme and no kidney forms [27,47,49]. In our computational model GDNF is part of the core patterning mechanism and in its absence no pattern or branching is observed. Gdnf +/− mutants have reduced branching and reduced expression levels of Gdnf, Ret and Wnt11 [27,29]. The computational model reproduces these mutants and exhibits reduced Gdnf expression levels ( Figure 4B), reduced GDNF signalling ( Figure 4C) (which correlates with the Ret and Wnt expression levels (Figure 2), and an increased wavelength of the Turing pattern ( Figure 4D), as characteristic for reduced branching. The Gdnf expression levels are below 50% in the simulations because of the lower levels of positive WNT feedback in the heterozygous mutant.

Ret and Gfr α Mutants
Ret −/− and Gfr α −/− mutant mice fail to develop kidneys [47,50,51] while Ret +/− mutant mice develop kidneys of normal size. In our computational model the receptors RET and GFRα1 are part of the core patterning mechanism and in their absence no pattern or branching is observed. The model predicts for Ret +/− mutants a small decrease in the expression and signalling levels ( Figure 4B,C) as well as slightly reduced branching ( Figure 4D). Experiments show that heterozygous Ret mutants demonstrate normal [20] or slightly reduced branching [52] (Figure 4 E). Note that Ret and Gfr α mutants are computed according to a model, which considers explicitly both RET and GFRα1 receptors (Equation 6).
To simulate the Wnt11 −/− mutants we set ρ W and ρ W0 to zero. Much as in the embryo, the extent of branching is reduced in Wnt11 mutants ( Figure 4D). Unlike in the lung where the distance between buds increases abruptly as the Fgf10 expression levels fall below a threshold in the Fgf10 allelic sequence [9], there is a gradual decrease of branching in the Ret/Wnt11 allelic sequence (4E) [20]. Both the lung [6] and the kidney phenotype ( Figure 4D) are reproduced by the respective models.
The model considers only one WNT ligand, WNT11, which enhances Gdnf expression. The other WNTs are likely to also enhance Gdnf expression, an effect that is incorporated implicitly in the WNT11independent Gdnf expression rate. We find that reducing either only the Wnt11 expression rate ρ W , only the WNT11-independent Gdnf expression rate ρ G0 , or both can all alter the branching mode in a synergistic manner ( Figure 4F). We further find that setting ρ W to zero reduces the Gdnf expression ( Figure 4G) and the GDNF-receptor complex concentration by about 20% (Figure 4H). A parallel reduction in ρ W and ρ G0 further reduces Gdnf expression, such that for a 20% reduction in ρ G0 and ρ W Gdnf expression and GDNF-receptor complex concentration are reduced to about 50% ( Figure 4G,H).

Ureteric Bud Outgrowth
Genetic analysis shows that the network that controls mesenchyme invasion by the ureteric bud is similar to the one controlling branching [5]. The ureteric buds of homozygous Gdnf, Ret or Gfrα1 knockouts fail to invade the metanephric mesenchyme [5,20,[26][27][28][29]. WNT11 seems to be less important because an invasion and a first round of branching are still observed in homozygous Wnt11 knockouts [20]. Figure 5A shows the computational domain that was used to represent the Wolffian duct, from which the ureteric bud branches out. The core regulatory network depicted in Figure 1A is sufficient to control ureteric bud outgrowth and a first round of branching in silico ( Figure 5B). Recently it was reported that the thickness of the mesenchyme of the embryonic lung plays an important role during branching morphogenesis, in that it has an impact on the outgrowth of the epithelium [2]. The thickness of the mesenchyme is indeed also important in our simulations of kidney branching morphogenesis: the thicker the mesenchyme, the shorter the stalk, and the earlier the first round of branching is observed ( Figure  5C). If the mesenchyme thickness exceeds some threshold value then the mesenchyme will be invaded in more than one place, similar to the case shown in Figure 5D. The impact of the mesenchyme thickness can be rationalized in the following way: the dimensionless diffusional length scale of GDNF in our model is l ∼ (2D G ln(2)/δ G ) 1/2 = 3.5. The thickness of the mesenchyme is around 0.5. Therefore, the thicker the mesenchyme the more GDNF can be supplied to the epithelium-mesenchyme border where GDNF binds to its receptors. An increased mesenchyme thickness has therefore an effect similar to an enhanced GDNF production rate.
Interestingly, mutants with reduced expression levels of the inhibitor of GDNF signalling, Sprouty, develop multiple kidneys [31,56]. Removal of an inhibitor of GDNF signalling corresponds to a reduction of the signalling threshold K R2G in our model. In the non-dimensional model a reduction in K R2G corresponds to an increase in the Ret and Gdnf production rates, ρ R , ρ G0 , and ρ G respectively (Equation 5).To address the Sprouty mutants in our simulation we increased the production rate of Gdnf, ρ G0 . Much as in the experiments, the simulations with an increased production rate of Gdnf indeed resulted in the invasion of the mesenchyme by the epithelium in several places ( Figure 5D). Both the thicker mesenchyme and the removal of Sprouty thus correspond to increased Gdnf production in our non-dimensional model. Sprouty has not been included into the model explicitly because it is an intracellular protein, while the model is focused on spatio-temporal patterning and therefore includes only secreted ligands and their receptors.

Discussion
While many regulatory components and local interactions have been defined an integrated understanding of the regulatory networks that control the different branching processes is lacking. We have shown here that the biochemical interactions between the core regulatory factor GDNF and its receptor result in a Schnakenberg-type Turing mechanism. The model is similar to the one that we previously obtained for the interaction of SHH and PTCH1 in the developing lung [6]. An important difference lies in the feedback architecture in that GDNF and WNT11 engage in a positive feedback while SHH forms a negative feedback with FGF10 in the lung. Moreover, while the diffusible ligands and receptors are expressed in separate tissue layers in both developing organs, Shh is expressed in the epithelium while Gdnf is expressed in the mesenchyme. Interestingly, while both mechanisms give rise to the full range of patterns that would correspond to the observed branching types (lateral branching, bifurcations, trifurcations) the GDNF/WNT11 based mechanism favours bifurcations and trifurcations while the SHH/FGF10 mechanism favours lateral branching events. This is in good agreement with the observed branching patterns in the two organs.
Unlike our model for lung development we solved the kidney model on growing domains that deformed in 3D in response to the local GDNF-receptor concentration. As the domain length increased further, patterns emerge, as is expected for a Turing mechanism. We further noticed that not only the length of the bud, but also the thickness of the mesenchyme affects the distance between two branching events in that the length of the outgrowing stalk is shorter and branching starts earlier if the mesenchyme is thicker. This is in good agreement with experimental observations in the lung where side-branching was noted to occur when sufficient space becomes available around the circumference of a parent branch [2].
It is so far unknown how the sequence of branching events is controlled during organogenesis, but given the stereotyped nature of early branching events, stochastic effects are unlikely to play a role and the branching sequence should emerge from the patterning mechanism in a deterministic fashion. We have previously shown in our model of lung branching morphogenesis that the signalling spots that emerge during bud outgrowth do not all have the same intensity [6]. We further showed that the growth speed can affect the branch mode [6]. If the growth speed depended on the strength of the signal then this could, in principle, result in a sequence of different branching events during bud outgrowth. In this manuscript we also show that the thickness of mesenchyme can influence branch mode selection. As kidney branching progresses the amount of mesenchyme surrounding the epithelial tip changes. In principle, this could influence branch mode selection. More work is clearly needed, both experimentally and computationally, to address the mechanism of branch mode selection.
Turing mechanisms have been proposed for many other biological patterning phenomena, and reproduce the size and geometry-dependence of biological patterns of various complexity [7,57]. However, it has remained difficult to firmly establish their use in biological pattern formation [58], and in several cases Turing-type mechanisms have been wrongly assigned to patterning processes such as e.g. the mechanism by which the stripy expression pattern of pair-rule genes emerge during Drosophila development [59]. These failures reveal the importance of a careful and comprehensive analysis of the underlying molecular interactions before proposing a Turing mechanism. The emerging Turing patterns are highly sensitive to the particular type of biochemical interactions and the parameter values. This aspect can be used to thoroughly test proposed Turing mechanisms with data from mutants. We have previously shown that the SHH/FGF10-based Turing mechanisms recapitulates even counterintuitive mutant phenotypes such as the abrupt increase in the spacings between buds in the Fgf10 allelic sequence as the Fgf10 expression levels fall below a threshold [6], and we have shown here that also the GDNF/RET-based Turing mechanism recapitulates all observed mutant phenotypes in the ureteric bud, including the gradual decrease of branching in the Ret/Wnt11 allelic sequence.
A number of alternative models have been proposed to explain the control of branching in other organs, that focused mainly on geometric rather than signalling effects, i.e. [60][61][62][63]. It is likely that both geometric and signalling effects work together to establish the observed patterns, but it is to be expected that signalling presents the dominating regulatory control.
In summary, Turing mechanisms offer a reliable mechanism for symmetry breaking that allows the outgrowing bud to change from its approximately cylindrical symmetry to the rotational symmetry after undergoing branching. The various branching events, i.e. bifurcations, trifurcation or lateral branching, can all result from the same regulatory interactions, thus permitting the implementation of complex branching pattern with a limited number of proteins. Different feedback architectures can modulate the branching type and sequence. At the same time different signalling systems can be used to implement such ligand-receptor based Turing mechanisms, such that branching in the lung can be explained with a receptor-ligand based Turing mechanism based on FGF10 and / or SHH, and branching in the ureteric bud can be explained with receptor-ligand based Turing mechanism based on GDNF. Unlike classical activator-inhibitor Turing mechanisms, ligand-receptor-based Turing mechanisms can be implemented by a single ligand, as long as the ligand-receptor interaction is cooperative, ligand-receptor binding results in an increased emergence of receptors on the membrane, and ligand diffuses faster than its receptor. This makes them a versatile mechanism for spontaneous pattern formation during development.

Non-Dimensionalization of the Model
1-3 were non-dimensionalized to reduce the total number of parameters: where, τ = tD G /r 2 1 ; The non-dimensional equation for the co-receptor GFRα1 readṡ and in 5 we then need to write R × G × G instead of R 2 × G for the concentration of the ligand-receptor complex. The parameters and their values are summarized in Table 1. In the limit of W 2 >> 1 and W 2 << 1 the partial differential equations (PDEs) that describe the dynamics of the variables G and R are independent of the variable W and the system of equations 5 reduces to the Schnakenberg model.

Linear Stability Analysis
To confirm that the observed patterns result from a Turing mechanism we performed a linear stability analysis at the steady state for the parameter values in Table 1. To that end we determined the eigenvalues of the Jacobian J in the absence and presence of diffusion and show that all real parts of the eigenvalues are negative in the absence of diffusion so that the steady state is stable in the absence of diffusion. In the presence of diffusion we obtain at least on eigenvalue with positive real part such that we have a diffusion-driven instability. The Jacobian for the PDE system given by 5 at the steady-state (G=0.434, R=3.13, W=34.4) in the absence of diffusion is given by: All eigenvalues of this Jacobian have a negative real part: λ 1,2 = −2.15 ± 2.27i, λ 3 = −0.500. In the presence of diffusion the Jacobian has at least one eigenvalue with positive real part, λ = 0.111, with wavenumber k = 6.

Numerical Solution of PDEs
The partial differential equations (PDEs) were solved with a finite element method as implemented in COMSOL Multiphysics 4.2a and 4.3. Details on how these models are implemented in COMSOL have been described by us recently [64,65]. COMSOL Multiphysics is a well-established software package and several studies confirm that COMSOL provides accurate solutions to reaction-diffusion equations both on constant [66] and growing domains [67][68][69]. We have previously used COMSOL to solve a similar mathematical set of PDEs on a bud shaped domain [6]. Table 3 shows that the converged solution differs by at most 5% as the mesh is further refined. When simulations were performed on an open domain the COMSOL 'bulk solution condition' was implemented at a distance 6 √ Dτ max from the mesenchyme, where τ max is the maximum time of model evaluation. It has previously been shown that beyond this, the effects of diffusion is not important on the experimental time scale [70]. When simulating patterns on a growing domain, growth was prescribed to be normal to the boundary and proportional to the local level of signalling v g = nGR 2 , where n denotes unit vector normal to the surface [71].

Robustness to Parameter Variability.
To estimate the robustness to spatial parameter variability we used a similar approach as in [72]. Thus, parameter values were assumed to be given by the formula k = k 0 × (1 + ξ θ (x, y, z))), where ξ θ (x, y, z) is normally distributed random function with a mean value of zero and half width θ. The half width of the distribution was equal for all parameters, except geometrical which were not varied. Tables   Table 1. Values of Dimensionless Parameters used for the Simulations. The domain geometry parameters were: l ep = 0.15, r 0 = 0.5, h 0 =0.6. For those parameters marked by star (*) model read-outs were virtually independent of the exact value of these parameters.
Calculations were carried out with the model from Figure 2D. Note that the spot at the tip vanishes upon addition of 5% Gaussian noise. B: The trifurcation branching pattern is robust even to the addition of 50% Gaussian noise to the parameter values. The grey scale depicts the concentration of the R 2 G signalling complex in the epithelium (white highest, black lowest concentration); the mesenchyme is not shown, but was included in the simulations.  Table 1, except ρ G0 =0.38. Colour code as in panel B.