An interplay of geometry and signaling enables robust lung branching morphogenesis

To achieve a large area of gas exchange within a limited thorax volume, lung development must be tightly controlled (Weibel, 1991). Embryonic development of the lung is indeed stereotyped to an extent that random development is unlikely. The lung tree is built by the sequential use of mostly three branching modes, namely lateral branching, planar bifurcations and orthogonal bifurcations (Metzger et al., 2008), and by rare trifurcations (Blanc et al., 2012). The sequence of these branching events is fixed in embryos from the same genetic background, and few errors are observed in wild-type littermates (Metzger et al., 2008), although a certain level of variation in the position and direction of branch outgrowth has been documented (Blanc et al., 2012; Short et al., 2012).

A number of models have been proposed to explain the control of the branching events. The earliest models focused on physical forces (reviewed by Lubkin, 2008), and combinations of experimental and computational studies have since confirmed that mechanical stress and internal pressure influence branching morphogenesis (Gjorevski and Nelson, 2010,, 2012; Kim et al., 2013; Nelson and Gleghorn, 2012; Unbekandt et al., 2008; Varner and Nelson, 2014). WNT signaling affects the epithelial shape of new lung buds, but WNT signaling is not essential for lung branching morphogenesis and is thus not part of the core regulatory network (Kadzik et al., 2014).

Other models focus on diffusion-based signaling effects because the secreted, diffusible proteins fibroblast growth factor 10 (FGF10) and sonic hedgehog (SHH) are necessary for lung branching morphogenesis (Abler et al., 2009; Bellusci et al., 1997a; Chiang et al., 1996; Peters et al., 1994; Weaver et al., 2000). As FGF10 signaling is necessary for the outgrowth of branches, it is generally assumed that FGF10 must vanish at the tip and concentrate on the sides adjacent to the lung bud tip to induce a splitting of the tip during bud outgrowth, thus resulting in bifurcating outgrowth. Lateral branching then requires the emergence of FGF10 signaling in multiple spots along the lung bud, while maintaining FGF10 signaling at the tip to support further outgrowth. It is a long-standing question how the expression and signaling patterns arise in the developing lung. A number of geometry-driven mechanisms have been proposed to explain how bifurcations may be directed. One such proposed mechanism is based on the distance between the FGF10-producing (i.e. the sub-mesothelial mesenchyme) and the FGF10-sensing (i.e. the epithelium) tissue (Bellusci et al., 1997b). It was noted that the closer the two tissues then the steeper the diffusional gradient should be, if the concentration were homogeneous in a given tissue layer (Clément et al., 2012a). If cells were responding to the FGF10 gradient rather than to the FGF10 concentration this mechanism would support bifurcations (Clément et al., 2012a,,b). In an alternative model it was proposed that FGF10 accumulates at the sides and drives bifurcating outgrowth because the tip of the Shh-expressing lung bud epithelium grows closer to the mesothelium. FGF10 and SHH engage in a negative feedback in that FGF10 signaling induces Shh expression (Abler et al., 2009), while SHH signaling represses the expression of Fgf10 (Bellusci et al., 1997b). The higher local SHH concentration would then inhibit Fgf10 production in the Fgf10-expressing sub-mesothelial mesenchyme (Bellusci et al., 1997b). In a computational implementation of the model, the mesothelium had to act as a diffusion barrier and lower SHH concentrations had to induce rather than inhibit Fgf10 expression (Hirashima and Iwasa, 2009), both of which remain to be demonstrated.

A further mechanism has been proposed based on the observation that the expression of key ligands is restricted to parts of the tissue, and certain patterns emerge when the soluble signaling proteins diffuse away from the producing tissue (Nelson et al., 2006). Using 3D shapes of lung bud epithelia extracted from early developing chicken lungs, visual inspection suggested that the areas where branching of secondary bronchi is inhibited coincide with where a steady-state diffusion model would predict high concentrations of diffusible inhibitory molecules when these are secreted from the epithelium into a large computational bounding box (Gleghorn et al., 2012).

We have recently shown that a ligand-receptor-based Turing mechanism can result in FGF10 signaling patterns that correspond to either lateral branching or to bifurcations, and recapitulates even counterintuitive mutant phenotypes such as the abrupt increase in the spacing between buds in the Fgf10 allelic sequence as the Fgf10 expression levels fall below a threshold (Celliere et al., 2012; Menshykau et al., 2012). Using a simplified geometry, others have since shown that FGF10-dependent Turing mechanisms can, in principle, also support the outgrowth of branches (Guo et al., 2014a,,b). Interestingly, although Fgf10 expression is necessary for branching morphogenesis (Abler et al., 2009), it has recently been shown that branching is still observed when Fgf10 is expressed homogenously, although the branching pattern is then different (Volckaert et al., 2013). Similarly, it is well established that lung epithelium will branch in vitro in the absence of mesenchyme if FGF is provided (Nogawa and Ito, 1995). Both experimental results contradict earlier models that were based on the distance of the Fgf10-expressing domain and the epithelium, while they are in agreement with a ligand-receptor-based Turing mechanism, because Turing mechanisms can yield patterns from a homogenous (noisy) distribution of the components without any need for a pre-pattern.

Here, the terms on the left-hand side of Eqs 1 and 2 are the time derivatives. The first term on the right-hand is the diffusion term, with D>1, since the ligand diffuses faster than its receptor. Although the mode of transport for morphogens is still a matter of debate (Müller et al., 2013) and often only a small fraction of morphogens may diffuse freely in the extracellular matrix (Zhou et al., 2012), we note that diffusion-based transport in combination with a realistic description of the receptor dynamics and ligand turnover has previously been shown to faithfully recapitulate the observed patterning process in a number of different developmental systems (Fried and Iber, 2014; Lopez-Rios et al., 2014; Nahmad and Stathopoulos, 2009). The non-dimensionalized reaction kinetics have three parameters: γ, a and b. γ is a scaling factor that influences the number of spots that can emerge on a domain; a and b are the constitutive production rates of ligand and receptor, respectively. −L and −R describe the linear decay of ligand and receptor, while −R²L describes receptor-dependent decay of ligand. The term +R²L describes the combined effects of ligand-triggered receptor turnover and ligand-induced receptor expression. In deriving this formula, a quasi-steady-state approximation was made for the concentration of the ligand-receptor complex, and the signaling complex was approximated by R²L (Badugu et al., 2012; Menshykau and Iber, 2013; Menshykau et al., 2012). We have previously shown that both the FGF10-receptor and SHH-receptor interactions are well described by Eqs 1 and 2 (Kurics et al., 2014; Menshykau et al., 2012).

A sequence of 3D lung geometries at different murine embryonic stages has recently been published (Blanc et al., 2012) and now permits for the first time the testing of the different proposed mechanisms with murine embryonic growth data. To that end we determined the displacement fields between four subsequent developmental stages. We then simulated the different models on the embryonic domains and checked whether the predicted signaling domains would coincide with where the lung bud actually grows out. We repeated the same procedure with time-lapse data for embryonic lung explants. In both cases, we show that only for a particular ligand-receptor Turing mechanism, but not for any of the alternative mechanisms studied, the signaling patterns coincide with the areas of growth. We further show in simulations that morphogen distributions that arise from the tissue-specific expression of morphogens (Nelson et al., 2006) are unstable under deforming outgrowth because the curvature would change as buds start to grow out at the sides. This mechanism thus creates split concentration profiles, but does not support the emergence of bifurcations of the bud tip as the bud is growing out. The Turing mechanism, by contrast, permits the emergence of such bifurcations in 3D simulations. The type of pattern that emerges via Turing mechanisms is typically very dependent on the (noisy) initial conditions. Importantly, the tissue-specific expression of ligand and receptor ensures patterning robustness of the ligand-receptor-based Turing mechanism and thus enables robust branching in 3D.

We used a published sequence of 3D lung geometries at four different murine embryonic stages between E11.25 and E11.75 [stages 1 to 4 (S1-S4) in Fig. 1A-D] (Blanc et al., 2012) to evaluate how well alternative mechanisms would predict the regions of lung bud outgrowth; image segmentation as well as surface and volume mesh generation were performed using the commercial software package Amira (Iber et al., 2015). In a first step, we calculated the displacement fields (Fig. 1E-G) between subsequent stages of the published lung geometries, using the landmark-based Bookstein algorithm (Bookstein, 1989) as implemented in Amira. We found that the lung buds predominantly grow at the tips and shrink in other places (red and blue arrows, respectively, in Fig. 1E-G). We focus on the signaling-dependent control of the outward movement of the epithelium, and do not attempt to predict the extent of shrinkage in the neighborhood, which may occur because of cell migration, deformation and rearrangements during the budding process (Kim et al., 2013). Accordingly, we set the length of all inward-pointing displacement vectors to zero (Fig. 1H-J). In the following, we refer to the resulting displacement vector field as the embryonic growth field. For every mesh point in the epithelium we next calculated the minimal distance to the mesenchyme surface using Amira. The distance between the epithelium and the mesenchyme surface is lowest at the tips of the lung buds (Fig. 1K; supplementary material Fig. S1A,B), but the distances do not vary across the lung bud as much as the embryonic growth field (Fig. 1L; supplementary material Fig. S1C,D).

To evaluate the alternative models that we discussed in the Introduction, and which are summarized in supplementary material Table S1, we determined the match between the embryonic growth fields and the predicted signaling domains on lung buds of the different stages over the likely physiological range of the parameter sets. At three different developmental stages (Fig. 1A-C) we subsequently calculated the deviation Δ (Eq. 4 in the Materials and Methods) of the normalized concentration of the ligand-receptor complex from the normalized length of the embryonic growth field shown in Fig. 1H-J. In other words, Δ quantifies the ‘overall’ deviation between the embryonic growth field and the predicted signaling patterns. We note that comparable results were obtained with different normalizations (supplementary material Text 1.1, Figs S2 and S3). The minimal Δ that corresponds to a good match between the predicted signaling field and the embryonic growth field varies between the stages: for the first transition from S1 (46 somites) to S2 (51 somites), the different models and parameter sets resulted in values of Δ of 0.62 and higher. Visual inspection showed that, for this stage, values of Δ up to 0.7 correspond to a good match between signaling pattern and growth field (supplementary material Fig. S4A). Note that the deviation of the embryonic displacement field to a signaling model with no signal would be equal to 1.

Both the FGF10 and SHH ligand-receptor systems can form independent Turing modules (Kurics et al., 2014). We will analyze the two Turing modules separately, starting with the FGF10-based Turing mechanism, without including the SHH module. In this case, the ligand FGF10 is represented by L, and its receptor FGFR2b is represented by R in Eqs 1 and 2. FGF10 signaling triggers outgrowth of the lung bud (Weaver et al., 2000), and the outgrowth of lung branches (but not of the main bud) is strictly dependent on FGF10 signaling (Peters et al., 1994). Accordingly, we expect branch outgrowth to occur in places where the concentration of the FGF10-FGFR2b signaling complex, R²L, is highest. Fgf10 is expressed only in the mesenchyme (blue areas in the cartoon in Fig. 2), while the FGF10 receptor Fgfr2b is expressed only in the epithelium (gray areas in the cartoon in Fig. 2) (Bellusci et al., 1997b). Accordingly, when we solved Eqs 1 and 2, the ligand expression rate, a, was non-zero only in the mesenchyme and the constitutive receptor expression rate, b, was non-zero only in the epithelium; receptor diffusion was also limited to the epithelium such that there are no FGFR2b receptors in the mesenchyme (see supplementary material Table S1, case T1 for details).

We next sampled the parameter values for a, b, d and γ and compared the simulated patterns with the measured embryonic displacement field. We obtained a wide range of patterns for the different parameter sets with a deviation (Δ, Eq. 4) between the predicted signaling levels and the embryonic growth field of 0.7 and above (Fig. 2A; supplementary material Fig. S4A). The best matching patterns could be further optimized with a local optimization algorithm to Δ=0.62 (Fig. 2A), and the predicted signaling domain (solid) overlays perfectly with the embryonic growth field (arrows). Differential growth and cellular responses can be expected only for a sufficient concentration difference within the domain. We therefore also analyzed the quality of the pattern as the maximal concentration difference within the domain, and we note that those parameters that yield the best match between the model and the embryonic growth field also result in the largest concentration difference inside and outside the signaling spots (Fig. 2, black spots).

Given the domain-dependent expression pattern of receptors and ligands, we could not carry out a linear stability analysis to define the parameter sets that give rise to a Turing instability, i.e. that are part of the Turing space. However, we note that such large concentration differences are typical for Turing patterns, and the parameter sets that yield such large concentration differences also exhibit other properties that are characteristic of Turing patterns (see supplementary material Text 1.2 and Fig. S4B-D).

The best (and the worst) matches between the signaling model and the physiological displacement field were thus all obtained for parameter sets that lay within this (inferred) Turing space (Fig. 2A, black dots). The best fit for a parameter set outside this (inferred) Turing space is considerably worse (Δ=0.8, supplementary material Fig. S4E) than the best fit obtained with parameters within the (inferred) Turing space (Δ=0.62, Fig. 2A).

Shh is essential for the development of the respiratory system and Shh-deficient lungs do not form branches (Pepicelli et al., 1998). We therefore examined whether an SHH-based mechanism could also predict the embryonic growth patterns. We have shown previously that the interaction of the ligand SHH (L in Eqs 1 and 2) with its receptor PTCH1 (R in Eqs 1 and 2) results in the ligand-receptor-based Turing model given in Eqs 1 and 2 (Menshykau et al., 2012) (see supplementary material Table S1, case T2 for details). The expression pattern of Shh/Ptch1 is the inverse of that of Fgf10/Fgfr2b in that the ligand Shh is expressed in the epithelium, whereas the receptor Ptch1 is expressed in the mesenchyme (Bellusci et al., 1997a). We therefore switched the position of the ligand and receptor expression domains such that the ligand would be produced in the epithelium (Fig. 2B, gray in the cartoon) and the receptor would be expressed in the mesenchyme (Fig. 2B, blue in the cartoon). The smallest deviation (Δ=0.9) is now obtained with a non-Turing pattern, but the overall match is poor (Fig. 2B, lower panel). Further local optimization barely improved the patterns. However, we note that, unlike in the case of FGF10, SHH signaling inhibits lung bud outgrowth as it represses the expression of Fgf10 (Bellusci et al., 1997b). Outgrowth should therefore be strongest where the concentration of the SHH-PTCH1 complex, R²L, is lowest, i.e. where the level of 1/R²L is highest. When we compared the level of 1/R²L with the embryonic growth field we obtained a match that was almost as good as for the FGF10-based Turing model as long as the parameters were within the (inferred) Turing space (Fig. 2C, black dots, Δ≥0.76). To judge the relevance of the two distinct tissue layers we carried out a further parameter screen in which we expressed both the receptor and the ligand in the epithelium (supplementary material Table S1, case T3). In this scenario, the best match of the Turing patterns is worse than that of the non-Turing patterns (Fig. 2D) and the overall best match is rather poor (Δ≥0.9, supplementary material Fig. S4F).

We repeated the analysis of FGF10 and SHH signaling models for the next two stages of lung development and obtained similar results (Fig. 2E-H), even though the lung geometry is more complicated at these later stages (Fig. 2E-H, lower panels). Finally, we acquired several movies of branching lung buds over 36 h of culture (Fig. 3A; supplementary material Fig. S5A), segmented the images and calculated the growth fields (see Materials and Methods for details). Again, we find that ligand-receptor-based Turing models, both for FGF10 (Fig. 3B-D; supplementary material Fig. S5B-D) and for SHH (Fig. 3E-G; supplementary material Fig. S5E-G) best predict the growth fields over time, even though the 2D branching patterns differ substantially from those observed in 3D in the embryo. This confirms the robustness of the proposed patterning mechanism. We conclude that both the FGF10-based Turing signaling model and the SHH-based Turing signaling model present excellent candidate mechanisms to control growth in the developing lung bud as long as the genes are expressed in their physiological domains.

Although the ligand-receptor-based Turing mechanism predicts the embryonic growth fields very well, the Turing spaces (i.e. the part of the parameter space where Turing patterns can be observed) for the core ligand-receptor mechanisms (FGF10 and SHH) are tiny (Fig. 4A). We have recently shown that negative feedbacks and the coupling of two Turing systems, as is the case for FGF10 and SHH, can massively increase the size of the Turing space (Kurics et al., 2014). We now examined whether this would also apply in case of the physiological model, in which ligands and receptors are expressed in different tissue layers. We note that the presence of subdomains in itself does not lead to a substantial change in the size of the Turing space (supplementary material Fig. S6) (Fujita and Kawaguchi, 2013). Owing to computational limitations when solving the models on the embryonic 3D domains, we could not explore the combined Turing space of the FGF10/SHH network, and we had to restrict our study to the effects of negative feedbacks within the FGF10 or SHH module that arise because of the SHH-FGF10 regulatory interaction. Accordingly, we studied one set of models with the receptor expressed in the epithelium, as characteristic for FGF10 (Fig. 4A,B; supplementary material Table S1, cases T1, TF1), and one set of models with the receptor expressed in the mesenchyme, as characteristic for SHH (Fig. 4C-F; supplementary material Table S1, cases T2, TF2-4).

In the first case (FGF10), the ligand-receptor complex triggers branch outgrowth, whereas in the second case (SHH) the complex prevents outgrowth. Since we only consider a single Turing system, we could implement a negative feedback only in the layer that expresses the receptor, and we therefore cannot include a negative feedback on the ligand production rate when the receptor is expressed in the epithelium. When the receptor is expressed in the mesenchyme, such negative feedback can be included though, because we approximate the thin epithelium by an infinitely thin layer (Menshykau and Iber, 2012). In both cases, we confirm the increase in the size of Turing space in the presence of additional negative feedbacks. Thus, the additional negative feedback on the receptor expression rate, a, results in an approximately 10-fold larger maximal receptor expression rate [compare the size of the inferred Turing space along the log(a)-axis in B and D with those in A and C, respectively, in Fig. 4]; the minimal receptor expression rate is zero. Similarly, the negative feedback on the ligand expression rate permits Turing patterns to emerge for a 10-fold larger range of the ligand expression rate, b [compare the size of the inferred Turing space along the log(b)-axis in C and E in Fig. 4]. A negative feedback on both the ligand and receptor expression rates results in a Turing space that is enlarged in the direction of both the ligand and receptor expression rates (compare the size of the inferred Turing space in C and F in Fig. 4). Based on our previous results (Kurics et al., 2014), we expect that substantially larger Turing spaces could still be obtained by coupling the FGF10 and SHH regulatory modules. Importantly, the Turing space not only widens, as observed previously (Kurics et al., 2014), but the physiological growth field can now be reproduced over a wider parameter range (i.e. the size of the red/yellow areas increases when negative feedbacks are introduced, Fig. 4). The quality of the fit changes more or less continuously as the parameter values are varied such that mutation and selection processes could have gradually improved the pattern by evolving the system into an optimal part of the parameter space from a first patterning solution.

To test the distance-based mechanism we need to restrict the expression of the ligand L to the outer boundary of the mesenchyme (by setting b in Eq. 3 accordingly; supplementary material Table S1, case A1). The ligand can, in principle, either trigger or inhibit branch outgrowth. If the ligand L triggers lung bud outgrowth the highest ligand concentration should coincide with the strongest displacement. For this case we obtain a deviation of Δ≥0.81 between the signaling model and the embryonic growth field as we screen the physiological range of D; visual inspection reveals a bad fit (Fig. 5A, lower panel). In particular, we always just observed two signaling spots in the simulations, whereas in the embryonic data there are three areas of outgrowth. If L were to act as an inhibitor, then we need to evaluate the match of 1/L and the displacement field, and this match is even worse (Fig. 5B, Δ≥0.9). It has been proposed that lung buds might respond to the local gradient rather than to the local concentration (Clément et al., 2012a). A gradient-based readout mechanism (supplementary material Table S1, case A2) does not improve the best match between model and embryonic growth field if the ligand is activating (Fig. 5C) and worsens the match when the ligand is inhibiting (Fig. 5D) lung bud outgrowth. We note that, in the case of a gradient-based readout, the quality of the match is independent of the parameter value (Fig. 5C,D). Based on the available data we can therefore rule out this mechanism.

We next tested the potential of tissue-specific protein expression to generate ligand patterns that would match the embryonic growth field. To this end we solved the model given by Eq. 3 with ligand production in the epithelium (b=1), but not in the mesenchyme (b=0) (supplementary material Table S1, case A3); ligand could diffuse everywhere. If the ligand L triggered lung bud outgrowth, the highest ligand concentration should coincide with the strongest displacement. However, we did not obtain a good match (Δ≥0.9) between signaling model and embryonic growth field as we screened the physiological range of D. Visual inspection of the best fit confirmed the bad match in that the strength of the growth field (colored arrows) and the intensity of the signaling field (solid colors) do not coincide (Fig. 5E, lower panel). If we assume that the ligand L acts as an inhibitor of growth and we evaluate the match of 1/L and the displacement field we get a slightly better match (Fig. 5F), but the deviation is still very high (Δ≥0.86) and visual inspection confirms a bad fit (Fig. 5F, lower panel). Including a positive feedback (Fig. 5G; supplementary material Table S1, case A4) or a negative feedback (Fig. 5H; supplementary material Table S1, case A5) of L on its own production does not improve the fit, even though we have a further parameter.

Tissue-specific protein expression is also incorporated in the previously analyzed ligand-receptor-based Turing mechanisms (Eqs 1 and 2). Outside the Turing space, the two models are thus very similar, except that in the Turing model a receptor is included. The deviations in the two models are indeed very similar outside the inferred Turing space, i.e. Δ≥0.9 with receptor (Eqs 1 and 2) and Δ≥0.9 without receptor (Eq. 3) if the ligand acts as activator, and Δ≥0.85 with receptor and Δ≥0.86 without receptor if the ligand acts as inhibitor, even though four parameter values rather than one parameter value were optimized in the Turing model. The higher number of parameters thus does not improve the fit.

Similarly, when we analyzed the alternative models with the 2D time-lapse data, we found that the global deviation Δ_γ for the alternative mechanisms (Fig. 5A′-F′) is always higher than that for the ligand-receptor-based Turing-type mechanisms (Fig. 3A,C).

The comparison of the predicted signaling domains and the embryonic growth fields supports a ligand-receptor-based Turing model (Eqs 1 and 2) to define the points of bud outgrowth. We next considered whether any of the proposed mechanisms would also support the outgrowth of a bud. We simulated a single tissue layer that is embedded in a 3D environment and that grows at the sites of strongest signaling in the direction perpendicular to its surface (for details see Materials and Methods). We started with the patterning mechanism that is based on the tissue-specific expression of ligand. As a result of the tissue-specific expression, the diffusible ligand will be lost at the edges and the highest concentration will therefore accumulate in the center of the domain (supplementary material Text 2 and Fig. S7A-C). A positive feedback can further enhance, whereas a negative feedback diminishes, this effect (supplementary material Fig. S7D,E). If this factor drives outgrowth of the domain, then a bud will emerge. More ligand is lost as the curvature of a domain increases (supplementary material Fig. S7F). As the length of the stalk increases, relatively more ligand will be lost at the curved tip, and the highest ligand concentration will therefore be found at the sides, resulting in a split localization of the diffusible ligand (supplementary material Fig. S7G,H). We tested a range of different growth functions to see whether the bifurcating ligand profile would support bifurcating outgrowth. However, none of these led to bifurcating outgrowth (Fig. 6A). The likely reason is that, as the curvature increases, more ligand is lost by diffusion (supplementary material Fig. S7F). As the bud would start to bifurcate, the curvature would increase locally and the concentration would diminish. As a result, outgrowth would come to an immediate halt. Furthermore, as shown in supplementary material Fig. S7B,G,H and Fig. 6A, the ligand distribution on a surface embedded into a 3D domain has rotational symmetry. Thus, tissue-specific ligand expression does not provide a mechanism to break radial symmetry, which is necessary for branching to occur. We conclude that tissue-specific ligand expression alone cannot drive the branching of a domain.

We next explored whether the Turing-type signaling mechanisms would support bifurcating outgrowth of the domain. To that end, we solved the ligand-receptor-based Turing model (Eqs 1 and 2) on a domain in the shape of a thin 2D disk, which was allowed to deform and grow within a 3D domain (see Materials and Methods for details). The Turing mechanism can indeed control the branching of a domain, and yields both bifurcating (Fig. 6B) and trifurcating (Fig. 6C) branch points. Bifurcating outgrowth can be observed on a closed domain (where ligand cannot diffuse away from the ligand-producing tissue layer), but the patterning mechanism is robust only on an open domain (where the ligand can diffuse away). On a closed domain, different patterns can emerge with the same parameter set depending on the noisy initial conditions, and only some of these patterns will support bud formation (Fig. 6D). The reason for the robust patterning on the open domain is the impact of the geometry, which concentrates the ligand in the center of the domain initially and thus strongly biases the Turing mechanism to this particular pattern.

We tested this observation on the embryonic lung geometry. When we remove the mesenchyme and express ligand and receptor both in the epithelium, the solution depends on the initial conditions, and for the same parameter values a wide range of different patterns is observed (Fig. 6E; supplementary material Table S1, case T4). The presence of the mesenchyme is thus important to stabilize the pattern. Interestingly, when lung epithelium is cultured in the absence of mesenchyme, random budding is indeed observed (Ohtsuka et al., 2001), whereas normal branching patterns are observed when both lung epithelium and mesenchyme are maintained in culture (Carraro et al., 2010). The Matrigel that surrounds the epithelium in the mesenchyme-free culture (Ohtsuka et al., 2001) is unlikely to stabilize the budding process because of the higher diffusion constant (Ciocan and Ciocan, 2009), which our analysis shows to reduce the geometry effect (supplementary material Fig. S7C).

It is a long-standing question how branching is controlled during development. Given the remarkably stereotyped nature of branching in the developing lung of wild-type littermates (Metzger et al., 2008), the branching mechanism must yield robust patterning in the same genetic background while allowing for differences in animals from a different genetic background without resulting in the failure of the entire branching program. Moreover, given the differences in the branching programs of different organs, such mechanism should be sufficiently flexible to permit organ-specific differences in branching programs. Finally, given the differences in the signaling networks that control branching in the different organs (Iber and Menshykau, 2013), it must be possible to implement such branching mechanisms with different signaling proteins and networks.

We have shown here that (unlike other proposed mechanisms) a ligand-receptor-based Turing mechanism in combination with tissue-restricted gene expression allows us to predict the growth fields of developing lungs buds. FGF10 or SHH or the combined FGF10/SHH network can each constitute the core Turing patterning mechanism in the developing lung. Given that the negative feedback between these two Turing systems greatly increases the parameter range for which Turing patterns are observed (the Turing space), as well as the range for which we can reproduce the embryonic growth pattern, it is, however, likely that the core mechanism is based on both signaling proteins, FGF10 and SHH. We note that the restriction of receptors to single cells as well as cooperative binding also increase the size of the Turing space (Kurics et al., 2014). The expression of ligand and receptors in different tissue layers is important to obtain the geometry effect that ensures robust pattern formation in spite of molecular noise. We therefore propose that the combination of geometry and signaling enables robust pattern selection and morphogenesis.

We note that a Turing mechanism alone can control the different branching modes (domain branching, planar and orthogonal bifurcations, trifurcations) and sequences (Iber and Menshykau, 2013), such that no ‘hierarchical and modular program that combines a small number of basic operations’, as previously conjectured by Krasnow and colleagues (Metzger et al., 2008), would be required. Given that many ligand-receptor pairs can give rise to a Turing system (Kurics et al., 2014), we propose that ligand-receptor-based Turing mechanisms together with tissue-restricted gene expression constitute a general mechanism to robustly control stereotyped branching as well as other patterning processes during morphogenesis (Badugu et al., 2012).

Shh-GFP-Cre (Harfe et al., 2004) was used to conditionally activate EGFP expression [β-actin-EGFP (Jagle et al., 2007)] in the lung epithelium. The mouse experiments were approved by the legally required regional commission in strict accordance with Swiss law. All studies were classified as grade zero implying minimal suffering of animals.

To follow lung branching morphogenesis in culture, E11.5 mouse embryonic lung rudiments were dissected and imaged as previously described; the culture media was additionally supplemented with bovine fetal serum (Carraro et al., 2010). Lung buds were imaged every 60 min using a Nikon Ti-E epifluorescence inverted microscope with a 4× lens.

We segmented the lung epithelium (GFP) and the mesenchyme (wide field) of the left lobe of the embryonic lung using standard MATLAB functions (see supplementary material Text 1.3 for details). The displacement fields between consecutive movie frames (separated by 2 h) were calculated as the set of vectors that are normal to the epithelium boundary in the current movie frame and that intersect the boundary in the next movie frame (Schwaninger et al., 2014). The growth fields were obtained by setting all vectors pointing inward (shrinkage) to zero.

The partial differential equations (PDEs) were solved on the imported 2D geometries and 3D computational meshes using COMSOL Multiphysics 4.3× as previously described (Menshykau and Iber, 2012; Vollmer et al., 2013). Several independent studies confirm that COMSOL provides accurate solutions to reaction-diffusion equations both on constant (Cutress et al., 2010; Kurics et al., 2014) and growing (Carin, 2006; Thummler and Weddemann, 2007; Weddemann and Thummler, 2008) domains.

Here, the bar indicates the average value in the domain.

Typically, 1000-10,000 parameter sets were first randomly sampled from a log-uniform distribution in the ranges: log₁₀(a): [−1 … 1], log₁₀(b): [−1 … 2], log₁₀(γ): [−4 … −1], log₁₀(D): [1 … 3], or similar. For the cases depicted in Fig. 3, between 10³ to 10⁴ parameter values were sampled. Sampled parameter sets with a minimal value of Δ were used as a starting point for the local minimization with the gradient-free coordinate search algorithm (Conn et al., 2009) as implemented in COMSOL Multiphysics 4.3a (Menshykau et al., 2013). We run this algorithm on the best 10-20 parameter sets obtained with random sampling. In all cases the coordinate search algorithm minimized Δ only by 0.02 or less, except for the case depicted in Fig. 3A, where Δ was reduced from 0.70 to 0.62.

We solved the PDE models on a deforming domain, where the deformation was normal to the surface and proportional to the local concentration, such that the velocity field was given by ⁠. Here, is the normal surface vector, v_g is the growth speed, and m accounts for any possible non-linear dependence on the signal concentration, c. In the case of geometry-based mechanisms c refers to the ligand concentration L, whereas in the case of the Turing mechanism c refers to the concentration of the receptor-ligand complex R²L. To incorporate noisy initial conditions we set the initial conditions to L=1+ξ(x,y,z,θ), where ξ(x,y,z,θ) is a normally distributed random function with a mean value of zero and half width θ.

We thank Rolf Zeller for access to lab space, mouse strains, consumables and support; Erica Montani and Thomas Horn for assistance with live imaging; Xin Sun for discussion; and Markus Affolter for critical reading of the manuscript.