Image Anal Stereol 2002;21:37-41 Original Research Paper A QUANTITATIVE METHOD FOR ANALYSING 3-D BRANCHING IN EMBRYONIC KIDNEYS: DEVELOPMENT OF A TECHNIQUE AND PRELIMINARY DATA Gabriel Fricout1, Luise A CULLEN-McEWEN2, Ian S Harper2, Dominique Jeulin1 and John F Bertram2 1Centre de Morphologie Mathématique, Ecole des Mines de Paris, France, 2Department of Anatomy and Cell Biology, Monash University, Clayton, Victoria, 3800, Australia e-mail: fricout@cmm.ensmp.fr, jeulin@cmm.ensmp.fr, Luise.Cullen@med.monash.edu.au, Ian.Harper@med.monash.edu.au, John.Bertram@med.monash.edu.au (Accepted. February 13, 2002) ABSTRACT The normal human adult kidney contains between 300,000 and 1 million nephrons (the functional units of the kidney). Nephrons develop at the tips of the branching ureteric duct, and therefore ureteric duct branching morphogenesis is critical for normal kidney development. Current methods for analysing ureteric branching are mostly qualitative and those quantitative methods that do exist do not account for the 3-dimensional (3D) shape of the ureteric “tree”. We have developed a method for measuring the total length of the ureteric tree in 3D. This method is described and preliminary data are presented. The algorithm allows for performing a semi-automatic segmentation of a set of grey level confocal images and an automatic skeletonisation of the resulting binary object. Measurements of length are automatically obtained, and numbers of branch points are manually counted. The final representation can be reconstructed by means of 3D volume rendering software, providing a fully rotating 3D perspective of the skeletonised tree, making it possible to identify and accurately measure branch lengths. Preliminary data shows the total length estimates obtained with the technique to be highly reproducible. Repeat estimates of total tree length vary by just 1-2%. We will now use this technique to further define the growth of the ureteric tree in vitro, under both normal culture conditions, and in the presence of various levels of specific molecules suspected of regulating ureteric growth. The data obtained will provide fundamental information on the development of renal architecture, as well as the regulation of nephron number. Keywords: 3-dimensional, branching, distance function, kidney development, skeleton, tree. INTRODUCTION Branching morphogenesis is an important event in the development of many organs including the lung, salivary gland, prostate gland, mammary gland and kidney. The pattern of branching within each of these organs is highly organised and plays a significant part in determining the final architecture of the organs. The normal human adult kidney contains between 300,000 and 1 million nephrons (the functional units of the kidney) (Nyengaard and Bendtsen, 1992). Nephrons develop at the tips of the branching ureteric duct. Evidence suggests that reduced nephron number is associated with the development of essential hypertension and chronic renal failure (Brenner et al., 1988; Brenner and Chertow, 1994) as well as the long-term success of renal allografts (Brenner and Milford 1993). Branching morphogenesis of the ureteric duct during kidney development is therefore not only critically important to normal kidney development, but may also underlie much subsequent renal pathology and abnormal physiology. Understanding the normal growth of this ureteric “tree” and the molecules that regulate its growth is therefore important in both basic and clinical renal science. Current methods for analysing ureteric branching are mostly qualitative. Those available quantitative methods do not take into account the 3-dimensional (3D) shape of the ureteric “tree”. We have developed a method for measuring the total length of the ureteric tree that takes into account the 3D structure. This method, based on the combination of 2D skeletons, is described and preliminary data are presented. 37 Fricout G et al: Analysing 3D branching: development of a new technique METHODS Embryonic mouse kidneys (metanephroi) were dissected at embryonic day 12 and cultured for up to 48 hours as previously described (Clark et al, 2001). Embryonic bodyweight was restricted to 0.060-0.080 g to reduce variation in kidney development prior to culture. At the end of the culture period the epithelial ureteric trees of whole kidneys were immunostained using an immunofluorescence technique (Fig. 1) that employs an Alexa 488 labelled antibody directed against Calbindin-D28K (Sigma-Aldrich Pty. Ltd., Castle Hill, Australia). Metanephroi were then optically sectioned with a 10x lens on a Leica TCS NT confocal microscope. Since the lens has a measured axial resolution (FWHM) of 10 µm, sections were taken at intervals of 4-5 µm to satisfy the requirements for Nyquist sampling. The 2D optical slices of the trees were 8 bit grey level images, typically 512 x 512 pixels. Between 30 and 60 images were required to sample each metanephros (Fig. 1A). As the skeletonisation process works on binary images, the algorithm first segments each confocal image (Fig. 1B) using classical morphological tools such as the watershed line of the gradient image (Serra, 1982; Meyer and Beucher, 1990). Fig. 1. A) four consecutive confocal sections through a metanephros cultured for 24 h B) corresponding binary images. Bar = 250 jdm. Some parts of the grey level images seem not to be segmented despite their important grey level value. This is due to the segmentation process which takes into account the grey level at the same location in all images in order to remove the overlapping information between frames coming from the confocal acquisition: the structures are only segmented on the frames where their grey level is the most important. The second and main part of the algorithm involves construction of a 3D skeleton of the tree, which is a representation of the original tree using thin lines connected through the centre of each branch. The definition of the skeleton will be here less restrictive than the standard definition of the skeleton as the set of centres of maximal balls included in the object. In this paper, a skeleton will be simply a connected set of pixels centred in the 3D object. The definition of a centred point will appear later in the paper, but it already appears that more than one set of pixels can match with this definition of a skeleton, in particular, a connected subset of a skeleton is still a skeleton in our definition. This problems deals with the difference between “real biological branches” and some noise coming from surface irregularities for example, so the relative flexibility in the definition of the skeleton used in this paper enables us to get finally a skeleton matching with the biological interpretation of the images. The method to construct this skeleton would be to find firstly some centred points in the tree and then to connect these points by a path that is also centred. Because of the definition of the confocal microscope, the scale of the pictures in the X and Y directions is very different from the Z scale of the 3D object, therefore the calculation of the centred points is made according to the X and Y direction on the one hand, and the Z direction on the other hand. A few definitions: I is the total 3D image, O is the bunary object obtained after the segmentation process, O is outside of the binary object obtained, F is used to designate one frame of 3D object. In all the following, O is supposed to be connected, if not, each step of the algorithm could be applied for each connected component of O. A path U between two voxels A and B is a finite part of IN p0 p1 pn such as p0 = A and pn = B and Vie[0..n-1\,pmeV(pi) where V(pi) designates the neighbours of pi according to a given connectivity. Then, the length of a path is the total number of points included in this path and is noted l(U). The set of all paths between A and B is called path(A, B). 38 Image Anal Stereol 2002;21:37-41 The distance function inside the tree is then defined by: \/MeOnF,dF(M) = min(l(U),Uepath(X,M),XeOf)F). (1) If the graph used corresponds to a vertical connectivity, that means, if the only neighbours of a voxel are the ones above and below, one obtains a Z distance function: \/M(x,y,z)eO,dz(M) = min(l(U),Ue path(X,M),XeOf)A), (2) where A is the line which has the equation {X = x and Y = y}. Then, the geodesic distance function inside the tree between a set A and a point M is: dg (M) = min(l(U),M0 e A,Ue path{M0,M)Wpe U,pe O). (3) The minimum distance between each point of the object and the outside of the tree on one frame F is calculated to give a function dF called the distance function of the object: the maxima of this function are called centred points in the x and y direction. The set of such maxima is called CF : let CF={MeOr\F/VXeV(M),dF(M)>dF(X)}. (4) As CF is usually not connected, one enlarges the definition of the centred point to the crests lines CLF of the distance function which link the previous maxima. let CLF={MeOnF/M is part of a crest line of dF}. (5) The union of these points gives a connected set for each connected component of each 2D frame of the 3D object (Fig. 2A). This set is in fact a 2D skeleton SF of the considered frame. This has been obtained according to the Zhou et al. (1998) algorithm. let SF=CF\JCLF. (6) One then obtains a 3D object by considering these 2D skeletons for all the frames: this 3D set is formed by all the points centred in the X and Y direction. let CXY = U F SF the union of all the SF. (7) The next operation follows logically. For each (x,y) coordinate, one calculates the points of the tree centred in the Z direction according to the distance dZ. The obtained set is usually composed of surfaces centred in the tree. let CZ={M(x,y,z)eO/X(x,y,t)eV(M)^dZ(M)>dZ(X)}, (8) then CXYZ = CXY[\CZ is the set of the centred points in O. (9) When taking the intersections between these two sets, one obtains the set of points centred in the final 3D object (Fig. 2B) which is a good start to construct a 3D skeleton. This set is, exactly like in 2D, usually not connected, and this problem has to be solved to be able to perform measurements on the tree. The way to process this is as follows. Starting from one point M0 of CXYZ: let M0 e CXYZ and K0 ={M0}. (10) One calculates for each point M of O the length of the minimum path inside O between M0 and M. This produces a map of the geodesic distance to the point M0. In fact this map is computed recursively by propagation starting from M0: the distance between M0 and its neighbours is determined, then the distance between M0 and the neighbours of the already computed points. The full distance map is not processed and one stops the calculation when the closest points M1 to M0 in CXYZ is reached. M1e{MeCXYZ/XeCXYZ^dg(K0,X)>dg(K0,M)}. (11) 39 Fricout G et al: Analysing 3D branching: development of a new technique The calculation of the path of minimum length is then processed by backpropagation in the obtained distance image: the neighbour of M1 whose distance to M0 has been processed and is minimum is kept as part of the path of minimum length between M0 and M1. This step is repeated from neighbours to neighbours until M0 is reached. This path, including M0 and M1 is used as a new starting set to compute a new geodesic distance map and find M2 in the same way as M1. 1- 0Upath g{ 0, 1j ; This iterative process enables us to connect all the points in a set called K (Fig. 2C). K could already be called a “3D skeleton of O”, but is still not a “good object” for obtaining meaningful measurements in the kidney because of two main problems: firstly, K is sometimes made of surfaces instead of single lines, and secondly, the bulbous ending observed in ureteric tree branches can induce some additional branches which are mostly biologically inaccurate. The next two steps of the algorithm try to correct these two problems. As said previously, a connected subset of a skeleton is still a skeleton and therefore, the aim of the next part is to find a connected subset K' of K in which noisy branches has been removed and the surfaces has been replaced by thin lines. Starting from one point M0 of K, the geodesic distance to M0 according to the mask K is performed where the maxima of this function are located at places where there are no points further from M0 in the neighbourhood. Therefore, a maximum of the geodesic distance to M0 corresponds to a non way path in the set K and this is the common definition used for the extremity of a branch. The same recursive process used for connecting the points of Cxyz according to O, is used for linking all these maxima according to the set K: the path of minimum length between an extremity and M0 is computed, then the path of minimum length between the obtained set and another extremity is computed and this is done recursively until all the maxima previously detected are linked to the point M0 by a single path. This process enables us to obtain a set K'cK in which the surfaces are removed and replaced by thin lines. The final process on K’ to obtain the skeleton S starts by the detection of branch points in K’: a point is defined as a branch point if it has more than three neighbours in K’. Usually, these branch points are grouped in connected sets of a few voxels that will be called branching sets. Thin lines called branches connect these branching sets. Starting from one branching set (called the root), it is possible to recursively define a tree structure where the internal nodes are branching sets and the edges are the branches which could be weighted according to their length. Two types of branches can now be defined: the branches connecting two internal nodes and the terminal branches which have only one extremity formed by a branching set whereas the other extremity is only a voxel with only one neighbour. These latter branches are the inaccurate ones due to the shape of the ureteric tree and we can remove them from the set K’ to obtain the skeleton S (Fig. 2D) in which the previously defined tree structure enables us to apply any kind of length measurement. However, some of the terminal branches are real biological branches and would therefore be ‘excessively pruned’. The algorithm has therefore been modified to keep all terminal branches if their length is greater than a limit fixed by the user. Fig. 2. A) 2D skeleton (points centred in x and y direction) of one frame superimposed on the corresponding binary image, B) intersection of the 2D skeletons and the maximum in z direction, C) the connected 3D skeleton D) the 3D skeleton after the cleaning process. Measurements are taken from this image. Bar = 250µm. RESULTS Fig. 3 shows the resulting skeleton superimposed on the original confocal images. The skeleton represents 40 Image Anal Stereol 2002;21:37-41 the three-dimensional ureteric tree; it is thin and centred following the branches. Preliminary data shows that length estimates obtained with the technique are highly reproducible. Repeated measures (n = 3) made on several embryonic day 12 mouse kidneys cultured for 24 h showed estimates of total tree length to vary by just 1 or 2% (Table 1). Fig. 3. Skeleton superimposed on the original confocal images. Bar = 250µm. Table 1. Repeat measurements on the same data set for two individual kidneys. Units are µm. Kidney Attempt 1 Attempt 2 Attempt 3 1 5345 5376 5295 2 6339 6476 6401 When applied to 10 individual embryonic kidneys, the length of the ureteric tree after 12 hours of culture was 2188 ± 390µm (mean ± SD), and at 24 h, ureteric tree length was 5,344 ± 696 µm. DISCUSSION The technique described allows the ureteric tree of the developing kidney to be analysed in 3D. The estimates of total tree length are highly reproducible. Preliminary use of this skeletonising technique, however, has identified various problems. Firstly, binarisation of overlapping branches, which either touch or lie close to each other, results in the joining of these branches in the skeleton. Presently, these branches need to be manually separated on the binary images by placing blank areas between the two. Secondly, as the grey level is not constant in all branches, some parts of the tree are sometimes removed and have to be manually filled in. Thirdly, when kidneys are only cultured for short periods (eg 6 hours) the branch tips are very dilated. When skeletonised, these dilated tips tend to form a loop, with the skeleton following the periphery of the dilated tip instead of stopping at the centre. To correct these problems, a new algorithm is being developed, which will operate directly in 3D, not on individual frames. This technique will now be used to further define the growth of the ureteric tree in vitro, under both normal culture conditions, and in the presence of increased and decreased levels of specific molecules suspected of regulating ureteric growth. The data obtained will not only provide fundamental information on the development of renal architecture, but will also provide understanding of the regulation of nephron number. ACKNOWLEDGMENT This research was supported by the Australian Kidney Foundation. REFERENCES Brenner BM, Milford EL (1993). Nephron underdosing: a programmed cause of chronic renal allograft failure. Am J Kidney Dis 21:66-72. Brenner BM, Chertow GM (1994). Congenital oligonephropathy and the etiology of adult hypertension and progressive renal injury. Am J Kidney Dis 23:171-5. Brenner BM, Garcia DL, Anderson S (1988). Glomeruli and blood pressure. Less of one, more the other? Am J Hypertens 4:335-47. Clark AT, Young RJ, Bertram JF (2001). In vitro studies on the roles of transforming growth factor-beta1 in rat metanephric development. Kidney Int 59:1641-53. Nyengaard JR, Bendsten TF (1992). Glomerular number and size in relation to age, kidney weight, and body surface in normal man. Anat Rec 232:194-201. Serra J (1982). Image Analysis and mathematical morphology. London: Academic Press. Meyer F, Beucher S (1990). Morphological Segmentation. J Visual Comm Image Represent 1:21-46. Zhou Y, Kaufman A, Toga AW (1998). Three-dimensionnal skeleton and centerline generation based on an approximate minimum distance field. The Visual Computer 14:303-14. 41