Image-based modeling of kidney branching morphogenesis reveals GDNF-RET based Turing-type mechanism and pattern-modulating WNT11 feedback.

Branching patterns and regulatory networks differ between branched organs. It has remained unclear whether a common regulatory mechanism exists and how organ-specific patterns can emerge. Of all previously proposed signalling-based mechanisms, only a ligand-receptor-based Turing mechanism based on FGF10 and SHH quantitatively recapitulates the lung branching patterns. We now show that a GDNF-dependent ligand-receptor-based Turing mechanism quantitatively recapitulates branching of cultured wildtype and mutant ureteric buds, and achieves similar branching patterns when directing domain outgrowth in silico. We further predict and confirm experimentally that the kidney-specific positive feedback between WNT11 and GDNF permits the dense packing of ureteric tips. We conclude that the ligand-receptor based Turing mechanism presents a common regulatory mechanism for lungs and kidneys, despite the differences in the molecular implementation. Given its flexibility and robustness, we expect that the ligand-receptor-based Turing mechanism constitutes a likely general mechanism to guide branching morphogenesis and other symmetry breaks during organogenesis.


Supplementary Notes 1 The Necessary Conditions for the Turing Mechanism
For the convenience for the reader, we summarise well established details about the Turing mechanism.

The Turing Mechanism
In this section we summarize the criteria for the emergence of a Turing pattern in reaction-diffusion systems with two species. We consider a general case of the ligand-receptor dynamics (Eq 1) systems in the form because they are easy to handle and have a biological interpretation (impermeable boundary). We note, however, that other boundary conditions would not greatly alter the following analysis. A Turing instability appears when a reaction-diffusion system has a stable steady state in the absence of diffusion, which loses its stability in the presence of diffusion such that spatial patterns emerge.

Linear stability in the absence of diffusion
Let R 0 and L 0 denote the steady state of the diffusion-free system of ordinary differential equations is the Jacobian evaluated at the point (R 0 , L 0 ). From now on, we write the partial derivatives evaluated at the steady state without their arguments for brevity. The steady state of the linearized system is

Diffusion-driven instability
Now let us add diffusion to our system of ODEs and consider the reaction-diffusion system linearized about the steady state w = (0, 0) T , which has the form w t = γJw + D∆w,  where D = diag(1, d) is a diagonal matrix containing the diffusion coefficients of the nondimensionalized system (Supplementary Equation-2). We look for a solution of the form where the exponents λ k determine the temporal growth of the solution and the time-independent functions W k are the solutions of the elliptic eigenvalue problem For instance, in one dimension on the interval [0, L] the eigenvalues are k = nπ/L (n = 0, 1, 2, . . .), also called wavenumbers, and the eigenfunctions are W (x) = cos(nπx/L) = cos(kx). The constants k ) T are the Fourier-coefficients of the initial conditions. Inserting equation (Supplementary Equation-6) into equation (Supplementary Equation-7) and using the fact that the set of eigenfunctions of the Laplace operator {W k } forms a complete orthonormal system, we obtain as linearized system for each wavenumber k. Writing where I = I 2 is the 2-by-2 identity matrix, we obtain the eigenvalues λ = λ k of the matrix M = γJ −k 2 D.
Expanding the above determinant, we obtain that λ k is the root of the second order polynomial equation Since we look for unstable solutions, we require that (λ k ) > 0 for some k = 0. This means that either the coefficient of λ and/or the constant term must be negative. Since the steady state is required to be linearly stable in the absence of diffusion (which corresponds to the case k = 0), we must have Hence, to obtain a λ with positive real part in the presence of diffusion we and the minimum value of h is which is negative if the expression in the bracket is negative.
In summary, the well-known conditions (see [2,Sec. 2.3]) for which a reaction-diffusion system with two species exhibits a Turing instability are as follows: where all partial derivatives are evaluated at the steady state (R 0 , L 0 ). We note that it is possible that these conditions are satisfied, but that no pattern emerges. This is the case when h is not negative for any k within the discrete set of wavenumbers, and only takes a negative value in between two of these discrete wavenumbers. The distance between wavenumbers shrinks as γ is increased, and in the limit of infinite γ the spectrum of k is continuous. Since γ is related to the size of the spatial domain, it follows that on small domains pattern formation may not happen, while on a sufficiently increased domain patterns may be observed.

The Necessary Conditions for a Turing Instability
The necessary conditions for a Turing instability follow from (Supplementary Equation-9). Thus, i) A difference in the diffusion coefficients, d > 1.
Substituting d = 1 into the equation (Supplementary Equation-9) we obtain f R + g L > 0 and f R + g L < 0, which cannot be both true. Therefore a difference in the diffusion coefficients is a necessary condition for a Turing instability.
ii) The presence of the positive feedback of the ligand-receptor signalling on the receptor abundance, i.e. v > mµ in Eq 1. In the absence of a sufficiently strong positive feedback Eq 1 can be written as: In this case, the condition df R + g L > 0 cannot be fulfilled as both f R = −R 0 − 2R 0 L 0 ≤ 0 and iii) Cooperative interactions between the receptor and the ligand (m = n when m + n = 2). If the ligand-receptor complex stoichiometry is one to one, then the ligand-receptor interactions are given by the following set of equations: Also, in this case not all the necessary conditions for the Turing instability can be fulfilled simultaneously.

Convergence and Accuracy of the Computational Method
As described in the Methods section we sampled the parameter space for four alternative models (T1-T4) from the log-normal distribution ( Supplementary Figures 3 (wt), 10 (FF) and 11 (FGS)).  The convergence test shows that the deviation, ∆, attains a minimum at b ≈ 10 1.6 (wt), b ≈ 10 2 (FF) and b ≈ 10 1.6 (FGS).

Model Validation with biochemical Perturbations
Uniform addition of GDNF ligand to the embryonic kidney culture results in a widening of the buds (Supplementary Figure 14A OPT Image Processing 1. The median filter [4] was applied to reduce noise in reconstructed OPT images.
2. Filtered images were segmented with the local Otsu filter [5] implemented in the Fiji Auto Local Threshold plugin collection [6]; the value for the radius of the local filter was set to 50 voxels.

Interbud Distance Analysis
Visual inspection of kidney images suggests that buds with minimal bud to bud distance are typically branched off from the same branching point (node) (Supplementary Figure 22D) and therefore reflect the distance buds have grown away from each other, rather than the minimal distance to which buds can approach each other. Given that bifurcations and trifurcations prevail during kidney branching morphogenesis [9], the third and further (Supplementary Figure 22D) shortest distances should represent an adequate measure for the distance by which buds can approach each other. Supplementary Figure   22E, F depicts the distribution and median distance to n-th buds. In an alternative approach the density of ureteric buds was estimated as the ratio of the kidney surface area to the number of epithelial tips.
The kidney surface area was approximated by the area of an ellipsoid that was fitted (ellipsoid fit [10] MATLAB function) to the coordinates of epithelial tips.

Supplementary Tables
Supplementary Table-1. A summary of the evaluated models. E and M stands for the epithelium and for the mesenchyme, accordingly.

Diagram
Equations Description R is expressed and diffuses in the epithelium. L is expressed in the mesenchyme, diffuses everywhere and binds to R in the epithelium.
R is expressed and diffuses in the epithelium. L is expressed in the mesenchyme, diffuses everywhere and binds to R in the epithelium.
R is expressed and diffuses in the epithelium. L is expressed in the mesenchyme, diffuses everywhere and binds to R in the epithelium.  E : R is expressed and diffuses in the epithelium. L is expressed in the mesenchyme, diffuses everywhere and binds to R in the epithelium. L 1 is expressed in the epithelium, diffuses everywhere and increases L expression in the mesenchyme. Supplementary 1 Dotted red line indicates source of GDNF at the outer mesenchymal boundary. n = 2, K = 20. 2 The mesenchyme in the model was virtually infinitely large, so that a further increase in the size of the mesenchyme did not affect the computational result.

Supplementary Figures
Supplementary Figure 1 (Table Supplementary Table-1: T1), red -a model without receptor up-regulation (Table   Supplementary Table-1: T2), green -a model where receptor and ligand diffusion coefficients were set equal (Table Supplementary Table-1: T3), and blue -a model with 1:1 stoichiometry of the ligandreceptor complex (Table Supplementary Table-  Comparison of growth areas as extracted from the imaging data (in grey) and the concentrations of (top row) the ligand-receptor complex, (center row) the receptor, and (bottom) the ligand as predicted by the complete ligand-receptor based model (Table Supplementary Table-1: T1, in black), a model without receptor up-regulation (Table Supplementary Table-1: T2, in red), a model where receptor and ligand diffusion coefficients were set equal (Table Supplementary Table-1: T3, in green) and a model with 1:1 stoichiometry of the ligand-receptor complex (Table Supplementary Table-1 (Table Supplementary Table-1: T1, black), performs better than any of the non-Turing mechanisms (red -a model without receptor up-regulation (Table  Supplementary Table-1: T2), green -a model where receptor and ligand diffusion coefficients were set equal (Table Supplementary Table-1: T3), and blue -model with 1:1 stoichiometry of the ligand-receptor complex (Table Supplementary Table-1: T4)).  (Table Supplementary Table-1: T1, black), performs better than any of the non-Turing mechanisms (red -a model without receptor up-regulation (Table  Supplementary Table-1: T2), green -a model where receptor and ligand diffusion coefficients were set equal (Table Supplementary Table-1: T3), and blue -model with 1:1 stoichiometry of the ligand-receptor complex (Table Supplementary Table-1: T4)).
Supplementary Figure 9. In Silico Morphogenesis, wild type Kidney. The non-Turing ligand-receptor based models (A) without positive feedback, T2, (B) without cooperative interactions, T3, and (C) with equal diffusion coefficients, T4, fail to induce branching morphogenesis. The initial shape of the computational domain qualitatively resembles that of wild type kidney explants and is the same as in the simulation for the ligand-receptor based model, T1 (Fig 3G). The model parameters are summarised in Table Supplementary Table-  field, E, extracted from the experimental data. Each coloured dot represents a deviation ∆, which has been calculated for a specific parameter set and at a given time. Black -the complete ligand-receptor based model (Table Supplementary Table-1: T1), red -a model without receptor up-regulation (Table   Supplementary Table-1: T2), green -a model where receptor and ligand diffusion coefficients were set equal (Table Supplementary Table-1: T3), and blue -a model with 1:1 stoichiometry of the ligandreceptor complex (Table Supplementary Table- (Table Supplementary Table-1: T1), red -a model without receptor up-regulation (Table   Supplementary Table-1: T2), green -a model where receptor and ligand diffusion coefficients were set equal (Table Supplementary Table-1: T3), and blue -a model with 1:1 stoichiometry of the ligandreceptor complex (Table Supplementary Table- , and (C) with equal diffusional coefficients, T4, fail to induce branching morphogenesis. The initial shape of the computational domain qualitatively resembles that of FF (F gf 10 −/− ) kidney explants and is the same as in the simulation for the ligand-receptor based model, T1 (Fig 3G). The model parameters are summarised in Table Supplementary Table- Table  Supplementary Table- Table Supplementary Table- Table Supplementary Table- Solving the ligand-receptor based Turing-type model on a growing domain. The upper and lower rows depict the concentration of the receptor and the ligand, accordingly. The vector fields indicate the displacement fields, the wireframe rendering represents the computational FEM mesh; note that the mesh used to carry out the computations reported in the manuscript was finer than the one depicted here, but finer meshes would appear all black at this resolution. →, and = denote solutions of a model on a growing domain, mapping of the solution on a FEM mesh, and solution mapped on a new mesh, respectively.  Blue circles depict kidney tips, the arrows indicate bud to bud distances. The red arrow indicates a distance between two tips that branch out from the same node.
(E) Histograms of the distances of a bud tip to the next closest, 2nd closest, 3rd closest, 4th closest, and 5th closest bud tip in wild type and mutant kidneys. (F) Median distance to the n-th closest bud.