Early branching events during lung development are stereotyped. Although key regulatory components have been defined, the branching mechanism remains elusive. We have now used a developmental series of 3D geometric datasets of mouse embryonic lungs as well as time-lapse movies of cultured lungs to obtain physiological geometries and displacement fields. We find that only a ligand-receptor-based Turing model in combination with a particular geometry effect that arises from the distinct expression domains of ligands and receptors successfully predicts the embryonic areas of outgrowth and supports robust branch outgrowth. The geometry effect alone does not support bifurcating outgrowth, while the Turing mechanism alone is not robust to noisy initial conditions. The negative feedback between the individual Turing modules formed by fibroblast growth factor 10 (FGF10) and sonic hedgehog (SHH) enlarges the parameter space for which the embryonic growth field is reproduced. We therefore propose that a signaling mechanism based on FGF10 and SHH directs outgrowth of the lung bud via a ligand-receptor-based Turing mechanism and a geometry effect.

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.

Turing mechanisms permit the self-organized emergence of a wide range of different patterns based on a diffusion-driven instability (Turing, 1952). They emerge for a particular network architecture (Gierer and Meinhardt, 1972; Prigogine, 1967; Prigogine and Lefever, 1968) and typically require at least two interacting factors that diffuse at substantially different rates, as is naturally the case for receptor-ligand systems. If ligands L and receptors R interact cooperatively and ligand-receptor binding results in an increased emergence of receptor on the membrane (by increased transcription, translation, recycling, less constitutive removal or similar) one obtains the standard Schnakenberg-type or activator-depleted substrate reaction kinetics for Turing patterns (Badugu et al., 2012; Gierer and Meinhardt, 1972; Kurics et al., 2014; Menshykau and Iber, 2013; Menshykau et al., 2012; Prigogine, 1967; Prigogine and Lefever, 1968; Schnakenberg, 1979; Tanaka and Iber, 2013), which take the form:
formula
1
formula
2

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 −R2L describes receptor-dependent decay of ligand. The term +R2L 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 R2L (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.

Model evaluation: comparison of predicted and real areas of growth on 3D lung bud shapes

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.

A ligand-receptor-based Turing mechanism correctly predicts the lung growth fields

An FGF10-based Turing mechanism

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, R2L, 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).

An SHH-based Turing mechanism

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, R2L, is lowest, i.e. where the level of 1/R2L is highest. When we compared the level of 1/R2L 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.

Negative feedbacks enlarge the Turing space for the control of lung outgrowth

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.

Alternative mechanisms for the control of branching morphogenesis

A number of alternative mechanisms have been proposed, as recently reviewed (Iber and Menshykau, 2013) (supplementary material Table S1, cases A1-A5). One of the earliest proposed mechanisms suggested that the distance between the Fgf10-expressing distal mesenchyme and the Shh-expressing epithelium would induce a pattern that could control branch point selection, because the repressive effect of SHH on Fgf10 expression would be stronger at a shorter distance (Bellusci et al., 1997b; Hirashima and Iwasa, 2009). More recently, the tissue-specific expression of ligands in either mesenchyme or epithelium has been shown to result in patterns, and this has been suggested to control the selection of branch points (Gleghorn et al., 2012; Nelson et al., 2006). In the following, we will test both mechanisms with the embryonic dataset. In both cases, we can represent the models by a single non-dimensional equation for the ligand concentration L:
formula
3
where the ligand production rate b is non-zero only in part of the domain, ΩL. When testing the suitability of the models, we varied the single parameter value, i.e. the non-dimensional diffusion coefficient D, over four orders of magnitude; the parameter range was adjusted such that at the low end the diffusion length scale is much smaller than the domain, whereas at the high end it is much larger.

Distance-based mechanisms

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.

Pattern emergence because of tissue-specific expression domains

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).

Branch outgrowth

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).

Mouse strains and ethics statement

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.

Two-dimensional time-lapse imaging

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.

Analysis of 2D lung movie data

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.

Numerical computations

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.

Deviation of predicted signaling patterns and measured displacement fields

The PDE models were solved for a wide range of parameter values. To evaluate the quality of the model predictions the L2 distance (Euclidean distance), Δ, between the computed signaling field and the registered displacement fields (areas of growth) was calculated using the formula:
formula
4
Here, EM refers to all mesh points in the interface between epithelium and mesenchyme. refers to all outward-pointing, normalized vectors of the measured lung displacement field. denotes the length of the vector. The displacement field vectors were normalized by the average vector length, such that the average length of all vectors is 1. Sn refers to the normalized computational signal; the signal was normalized by the average signal, such that the average of Sn is 1. In the case of an activating signal:
formula
5
while in case of an inhibiting signal:
formula
6

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: log10(a): [−1 … 1], log10(b): [−1 … 2], log10(γ): [−4 … −1], log10(D): [1 … 3], or similar. For the cases depicted in Fig. 3, between 103 to 104 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.

Deforming outgrowth of branches

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, vg 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 R2L. 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.

Author contributions

D.I. and D.M. conceived the study; D.M. carried out all computational analysis as well as 2D lung culturing and imaging; P.B. and V.S. provided the 3D embryonic data; E.U. performed all mouse work related to lung culturing. All authors approved the final version of the manuscript.

Funding

The authors acknowledge funding from a SystemsX and Sinergia SNF grant and an ETH Fellowship to D.M.

Abler
,
L. L.
,
Mansour
,
S. L.
and
Sun
,
X.
(
2009
).
Conditional gene inactivation reveals roles for Fgf10 and Fgfr2 in establishing a normal pattern of epithelial branching in the mouse lung
.
Dev. Dyn.
238
,
1999
-
2013
.
Badugu
,
A.
,
Kraemer
,
C.
,
Germann
,
P.
,
Menshykau
,
D.
and
Iber
,
D.
(
2012
).
Digit patterning during limb development as a result of the BMP-receptor interaction
.
Sci. Rep.
2
,
991
.
Bellusci
,
S.
,
Furuta
,
Y.
,
Rush
,
M. G.
,
Henderson
,
R.
,
Winnier
,
G.
and
Hogan
,
B. L.
(
1997a
).
Involvement of Sonic hedgehog (Shh) in mouse embryonic lung growth and morphogenesis
.
Development
124
,
53
-
63
.
Bellusci
,
S.
,
Grindley
,
J.
,
Emoto
,
H.
,
Itoh
,
N.
and
Hogan
,
B. L.
(
1997b
).
Fibroblast growth factor 10 (FGF10) and branching morphogenesis in the embryonic mouse lung
.
Development
124
,
4867
-
4878
.
Blanc
,
P.
,
Coste
,
K.
,
Pouchin
,
P.
,
Azaïs
,
J.-M.
,
Blanchon
,
L.
,
Gallot
,
D.
and
Sapin
,
V.
(
2012
).
A role for mesenchyme dynamics in mouse lung branching morphogenesis
.
PLoS ONE
7
,
e41643
.
Bookstein
,
F. L.
(
1989
).
Principal warps: thin-plate splines and the decomposition of deformations
.
Pattern Anal. Mach. Intell.
11
,
567
-
585
.
Carin
,
M.
(
2006
).
Numerical Simulation of moving boundary problems with the ALE method: validation in the case of a free surface and a moving solidification front
. In
Proceedings of COMSOL Conference 2006
.
Paris
.
Carraro
,
G.
,
del Moral
,
P.-M.
and
Warburton
,
D.
(
2010
).
Mouse embryonic lung culture, a system to evaluate the molecular mechanisms of branching
.
J. Vis. Exp.
40
,
ii: 2035
.
Celliere
,
G.
,
Menshykau
,
D.
and
Iber
,
D.
(
2012
).
Simulations demonstrate a simple network to be sufficient to control branch point selection, smooth muscle and vasculature formation during lung branching morphogenesis
.
Biol. Open
1
,
775
-
788
.
Chiang
,
C.
,
Litingtung
,
Y.
,
Lee
,
E.
,
Young
,
K. E.
,
Corden
,
J. L.
,
Westphal
,
H.
and
Beachy
,
P. A.
(
1996
).
Cyclopia and defective axial patterning in mice lacking Sonic hedgehog gene function
.
Nature
383
,
407
-
413
.
Ciocan
,
E.
and
Ciocan
,
R.
(
2009
).
Optimized numerical pharmacokinetics model for optical molecular probes based on diffusion coefficients in matrigel measured using fluorescence imaging
. In
Conf. Proc. IEEE Eng. Med. Biol. Soc.
2009
,
4925
-
4928
.
Clément
,
R.
,
Blanc
,
P.
,
Mauroy
,
B.
,
Sapin
,
V.
and
Douady
,
S.
(
2012a
).
Shape self-regulation in early lung morphogenesis
.
PLoS ONE
7
,
e36925
.
Clément
,
R.
,
Douady
,
S.
and
Mauroy
,
B.
(
2012b
).
Branching geometry induced by lung self-regulated growth
.
Phys. Biol.
9
,
066006
.
Conn
,
A. R.
,
Scheinberg
,
K.
and
Vicente
,
L. N.
(
2009
).
Introduction to Derivative-Free Optimization. Philadelphia
:
Society for Industrial and Applied Mathematics
.
Cutress
,
I. J.
,
Dickinson
,
E. J. F.
and
Compton
,
R. G.
(
2010
).
Analysis of commercial general engineering finite element software in electrochemical simulations
.
J. Electroanal. Chem.
638
,
76
-
83
.
Fried
,
P.
and
Iber
,
D.
(
2014
).
Dynamic scaling of morphogen gradients on growing domains
.
Nat. Commun.
5
,
5077
.
Fujita
,
H.
and
Kawaguchi
,
M.
(
2013
).
Pattern formation by two-layer Turing system with complementary synthesis
.
J. Theor. Biol.
322
,
33
-
45
.
Gierer
,
A.
and
Meinhardt
,
H.
(
1972
).
A theory of biological pattern formation
.
Kybernetik
12
,
30
-
39
.
Gjorevski
,
N.
and
Nelson
,
C. M.
(
2010
).
Endogenous patterns of mechanical stress are required for branching morphogenesis
.
Integr. Biol.
2
,
424
-
434
.
Gjorevski
,
N.
and
Nelson
,
C. M.
(
2012
).
Mapping of mechanical strains and stresses around quiescent engineered three-dimensional epithelial tissues
.
Biophys. J.
103
,
152
-
162
.
Gleghorn
,
J. P.
,
Kwak
,
J.
,
Pavlovich
,
A. L.
and
Nelson
,
C. M.
(
2012
).
Inhibitory morphogens and monopodial branching of the embryonic chicken lung
.
Dev. Dyn.
241
,
852
-
862
.
Guo
,
Y.
,
Chen
,
T.-H.
,
Zeng
,
X.
,
Warburton
,
D.
,
Bostrom
,
K. I.
,
Ho
,
C.-M.
,
Zhao
,
X.
and
Garfinkel
,
A.
(
2014a
).
Branching patterns emerge in a mathematical model of the dynamics of lung development
.
J. Physiol.
592
,
313
-
324
.
Guo
,
Y.
,
Sun
,
M.
,
Garfinkel
,
A.
and
Zhao
,
X.
(
2014b
).
Mechanisms of side branching and tip splitting in a model of branching morphogenesis
.
PLoS ONE
9
,
e102718
.
Harfe
,
B. D.
,
Scherz
,
P. J.
,
Nissim
,
S.
,
Tian
,
H.
,
McMahon
,
A. P.
and
Tabin
,
C. J.
(
2004
).
Evidence for an expansion-based temporal Shh gradient in specifying vertebrate digit identities
.
Cell
118
,
517
-
528
.
Hirashima
,
T.
and
Iwasa
,
Y
. (
2009
).
Mechanisms for split localization of Fgf10 expression in early lung development
.
Dev. Dyn.
238
,
2813
-
2822
.
Iber
,
D.
and
Menshykau
,
D.
(
2013
).
The control of branching morphogenesis
.
Open Biol.
3
,
130088
.
Iber
,
D.
,
Tanaka
,
S.
,
Fried
,
P.
,
Germann
,
P.
and
Menshykau
,
D.
(
2015
).
Simulating tissue morphogenesis and signaling
.
Methods Mol. Biol.
1189
,
323
-
338
.
Jagle
,
U.
,
Gasser
,
J. A.
,
Müller
,
M.
and
Kinzel
,
B.
(
2007
).
Conditional transgene expression mediated by the mouse beta-actin locus
.
Genesis
45
,
659
-
666
.
Kadzik
,
R. S.
,
Cohen
,
E. D.
,
Morley
,
M. P.
,
Stewart
,
K. M.
,
Lu
,
M. M.
and
Morrisey
,
E. E.
(
2014
).
Wnt ligand/Frizzled 2 receptor signaling regulates tube shape and branch-point formation in the lung through control of epithelial cell shape
.
Proc. Natl. Acad. Sci. USA
111
,
12444
-
12449
.
Kim
,
H. Y.
,
Varner
,
V. D.
and
Nelson
,
C. M.
(
2013
).
Apical constriction initiates new bud formation during monopodial branching of the embryonic chicken lung
.
Development
140
,
3146
-
3155
.
Kurics
,
T.
,
Menshykau
,
D.
and
Iber
,
D.
(
2014
).
Feedbacks, receptor clustering, and receptor restriction to single cells yield large Turing spaces for ligand-receptor based Turing models
.
Phys. Rev. E Stat. Nonlin. Soft Matter Phys.
90
,
022716
.
Lopez-Rios
,
J.
,
Duchesne
,
A.
,
Speziale
,
D.
,
Andrey
,
G.
,
Peterson
,
K. A.
,
Germann
,
P.
,
Unal
,
E.
,
Liu
,
J.
,
Floriot
,
S.
,
Barbey
,
S.
, et al. 
(
2014
).
Attenuated sensing of SHH by Ptch1 underlies evolution of bovine limbs
.
Nature
511
,
46
-
51
.
Lubkin
,
S. R.
(
2008
).
Branched organs: mechanics of morphogenesis by multiple mechanisms
.
Curr. Top. Dev. Biol.
81
,
249
-
268
.
Menshykau
,
D.
and
Iber
,
D.
(
2012
).
Simulating organogenesis with COMSOL: interacting and deforming domains
. In
Proceedings of COMSOL Conference 2012
.
Milan
.
Menshykau
,
D.
and
Iber
,
D.
(
2013
).
Kidney branching morphogenesis under the control of a ligand-receptor-based Turing mechanism
.
Phys. Biol.
10
,
046003
.
Menshykau
,
D.
,
Kraemer
,
C.
and
Iber
,
D.
(
2012
).
Branch mode selection during early lung development
.
PLoS Comput. Biol.
8
,
e1002377
.
Menshykau
,
D.
,
Shrivastsan
,
A.
,
Germann
,
P.
,
Lemereux
,
L.
and
Iber
,
D.
(
2013
).
Simulating organogenesis in COMSOL: parameter optimization for PDE-based models
. In
Proceedings of COMSOL Conference 2013
.
Rotterdam
.
Metzger
,
R. J.
,
Klein
,
O. D.
,
Martin
,
G. R.
and
Krasnow
,
M. A.
(
2008
).
The branching programme of mouse lung development
.
Nature
453
,
745
-
750
.
Müller
,
P.
,
Rogers
,
K. W.
,
Yu
,
S. R.
,
Brand
,
M.
and
Schier
,
A. F.
(
2013
).
Morphogen transport
.
Development
140
,
1621
-
1638
.
Nahmad
,
M.
and
Stathopoulos
,
A.
(
2009
).
Dynamic interpretation of hedgehog signaling in the Drosophila wing disc
.
PLoS Biol.
7
,
e1000202
.
Nelson
,
C. M.
and
Gleghorn
,
J. P.
(
2012
).
Sculpting organs: mechanical regulation of tissue development
.
Annu. Rev. Biomed. Eng.
14
,
129
-
154
.
Nelson
,
C. M.
,
VanDuijn
,
M. M.
,
Inman
,
J. L.
,
Fletcher
,
D. A.
and
Bissell
,
M. J.
(
2006
).
Tissue geometry determines sites of mammary branching morphogenesis in organotypic cultures
.
Science
314
,
298
-
300
.
Nogawa
,
H.
and
Ito
,
T.
(
1995
).
Branching morphogenesis of embryonic mouse lung epithelium in mesenchyme-free culture
.
Development
121
,
1015
-
1022
.
Ohtsuka
,
N.
,
Urase
,
K.
,
Momoi
,
T.
and
Nogawa
,
H.
(
2001
).
Induction of bud formation of embryonic mouse tracheal epithelium by fibroblast growth factor plus transferrin in mesenchyme-free culture
.
Dev. Dyn.
222
,
263
-
272
.
Pepicelli
,
C. V.
,
Lewis
,
P. M.
and
McMahon
,
A. P.
(
1998
).
Sonic hedgehog regulates branching morphogenesis in the mammalian lung
.
Curr. Biol.
8
,
1083
-
1086
.
Peters
,
K.
,
Werner
,
S.
,
Liao
,
X.
,
Wert
,
S.
,
Whitsett
,
J.
and
Williams
,
L.
(
1994
).
Targeted expression of a dominant negative FGF receptor blocks branching morphogenesis and epithelial differentiation of the mouse lung
.
EMBO J.
13
,
3296
-
3301
.
Prigogine
,
I.
(
1967
).
On symmetry-breaking instabilities in dissipative systems
.
J. Chem. Phys.
46
,
3542
-
3550
.
Prigogine
,
I.
and
Lefever
,
R.
(
1968
).
Scitation: symmetry breaking instabilities in dissipative systems. II
.
J. Chem. Phys.
48
,
1965
.
Schnakenberg
,
J.
(
1979
).
Simple chemical reaction systems with limit cycle behaviour
.
J. Theor. Biol.
81
,
389
-
400
.
Short
,
K.
,
Hodson
,
M.
and
Smyth
,
I.
(
2012
).
Spatial mapping and quantification of developmental branching morphogenesis
.
Development
140
,
471
-
478
.
Schwaninger
,
C. A.
,
Menshykau
,
D.
and
Iber
,
D.
(
2014
).
Simulating organogenesis: algorithms for the image-based determination of growth fields. ACM Transactions on Modeling and Computer Simulation (TOMACS) (in press).
Tanaka
,
S.
and
Iber
,
D.
(
2013
).
Inter-dependent tissue growth and Turing patterning in a model for long bone development
.
Phys. Biol.
10
,
056009
.
Thummler
,
V.
and
Weddemann
,
A
. (
2007
).
Computation of space-time patterns via ALE methods
. In
Proceedings of COMSOL Conference 2007
.
Grenoble
.
Turing
,
A. M.
(
1952
).
The chemical basis of morphogenesis
.
Phil. Trans. R. Soc. Lond. B
237
,
37
-
72
.
Unbekandt
,
M.
,
del Moral
,
P.-M.
,
Sala
,
F. G.
,
Bellusci
,
S.
,
Warburton
,
D.
and
Fleury
,
V.
(
2008
).
Tracheal occlusion increases the rate of epithelial branching of embryonic mouse lung via the FGF10-FGFR2b-Sprouty2 pathway
.
Mech. Dev.
125
,
314
-
324
.
Varner
,
V. D.
and
Nelson
,
C. M.
(
2014
).
Cellular and physical mechanisms of branching morphogenesis
.
Development
141
,
2750
-
2759
.
Volckaert
,
T.
,
Campbell
,
A.
,
Dill
,
E.
,
Li
,
C.
,
Minoo
,
P.
and
De Langhe
,
S.
(
2013
).
Localized Fgf10 expression is not required for lung branching morphogenesis but prevents differentiation of epithelial progenitors
.
Development
140
,
3731
-
3742
.
Vollmer
,
J.
,
Menshykau
,
D.
and
Iber
,
D.
(
2013
).
Simulating organogenesis in COMSOL: cell-based signaling models
. In
Proceedings of COMSOL Conference 2013
.
Rotterdam
.
Weaver
,
M.
,
Dunn
,
N. R.
and
Hogan
,
B. L.
(
2000
).
Bmp4 and Fgf10 play opposing roles during lung bud morphogenesis
.
Development
127
,
2695
-
2704
.
Weddemann
,
A.
and
Thummler
,
V.
(
2008
).
Stability analysis of ALE-methods for advection-diffusion problems
. In
Proceedings of COMSOL Conference 2008
.
Hannover
.
Weibel
,
E.
(
1991
).
Fractal geometry: a design principle for living organisms
.
Am. J. Physiol.
261
,
L361
-
L369
.
Zhou
,
S.
,
Lo
,
W.-C.
,
Suhalim
,
J. L.
,
Digman
,
M. A.
,
Gratton
,
E.
,
Nie
,
Q.
and
Lander
,
A. D.
(
2012
).
Free extracellular diffusion creates the Dpp morphogen gradient of the Drosophila wing disc
.
Curr. Biol.
22
,
668
-
675
.

Competing interests

The authors declare no competing financial interests.

Supplementary information