- Email: [email protected]
Voronoi tessellation to study the numerical density and the spatial distribution of neurones
Journal of Chemical Neuroanatomy 20 (2000) 83 – 92 www.elsevier.com/locate/jchemneu
Voronoi tessellation to study the numerical density and the spatial distribution of neurones Charles Duyckaerts a,*, Gilles Godefroy b a
Laboratoire de Neuropathologie R. Escourolle, Hoˆpital de La Salpeˆtrie`re, 47 Boule6ard de l’Hoˆpital, 75651 Paris Cedex 13, France b Equipe d’analyse, Uni6ersite´ Paris 6, Boıˆte 186, 4, place Jussieu, 75252 Paris Cedex 05, France
Abstract The conditions of regularity and isotropy, required by standard morphometric procedures, are generally not fulfilled in the central nervous system (CNS) where cells are distributed in a highly complex manner. The evaluation of the mean numerical density of neuronal or glial cells does not take into account the topographical heterogeneity and thereby misses the information that it contains. A local measurement of the density can be obtained by evaluating the ‘numerical density of one cell’, i.e. the ratio 1/(the volume that the cell occupies). This volume is the region of space that is closer to that cell than to any other. It has the shape of a polyhedron, called Voronoi (or Dirichlet) polyhedron. In 2-D, the Voronoi polyhedron is a polygon, the sides of which are located at mid-distance from the neighbouring cells. The Voronoi polygons are contiguous and their set fills the space without interstice or overlap, i.e. they perform a ‘tessellation’ that may yield a density map when the same colours are used to fill polygons of similar sizes. The use of Voronoi polygons allows computing the confidence interval of a mean numerical density that makes statistical comparisons possible. The tessellation also provides information concerning spatial distribution; the areas of the Voronoi polygons do not vary much when the cells are regularly distributed. On the contrary, small and large polygons are found when cellular clusters are present. The coefficient of variation of the polygon areas is an objective measurement of their variability and helps to define ‘regular’, ‘clustered’ and ‘random’ distributions. When cells are clustered, small polygons are contiguous and may be objectively identified by simple algorithms. Voronoi tessellations are easily performed in 2-D. On an average the area of a polygon times the thickness of the section equals the volume of the corresponding polyhedron. 3-D tessellations that are theoretically possible and for which algorithms have been published remain to be adapted to histological works. © 2000 Elsevier Science B.V. All rights reserved. Keywords: Dirichlet tessellation; Voronoi tessellation; Morphometry; Neuronal density; Neuronal loss
1. Introduction Consider a histological view of such complex structures as the thalamus or the cortex. The section of the cells (their ‘profiles’) are grouped in subnuclei or in layers where they appear to be variously packed. One wants to compare the same structure in various cases or, in the same case, to analyse the characteristics of different regions. Ideally, we would like to be able to state that a layer in the cortex or a subnucleus in the thalamus lacks neurones or, on the contrary, is more densely packed than its neighbour. A crude way to * Corresponding author. Tel.: + 33-1-42161891; fax: +33-14239828. E-mail address: [email protected] (C. Duyckaerts).
obtain this information is based on the simple counting of the cellular profiles in a region that has been delineated. Dividing the number of profiles by the area of the subregion provides a mean ‘numerical density’ usually a number of cellular profiles per square mm. Notwithstanding the stereological problems that are raised by such a procedure (which we will consider later), average density values miss significant information; if, for instance, the cortex has a density of 1000 neuronal profiles per mm2 rather than 1200, this provides no information on the local characteristics of the changes, for instance, is layer IV affected more than layer V? Assessing the profile density in smaller areas, for example, in the various layers of the cortex is possible but raises other difficulties; the limits of the neuronal groupings such as layers or subnuclei are often difficult to define objectively since their identifica-
0891-0618/00/$ - see front matter © 2000 Elsevier Science B.V. All rights reserved. PII: S 0 8 9 1 - 0 6 1 8 ( 0 0 ) 0 0 0 6 4 - 8
84
C. Duyckaerts, G. Godefroy / Journal of Chemical Neuroanatomy 20 (2000) 83–92
tion itself relies on a subjective evaluation of numerical density. The error that is made in drawing the borders has a very significant impact on the final value, especially, when the structure is small. To improve the measurement, one would like to be able to evaluate objectively the ‘local density’ of profiles, without having to delineate the layers or the subnuclei whose identification is based on slight differences in packing density of the cellular profiles. We would also like to be able to consider a small group of neurones and immediately tell what their numerical density is. How small can this group be? Could we find, so to speak, the numerical density of one neuron? Its ‘density’ would then be 1/the area in which it stands. What is this area? If each neuronal profile has its own area, their sum should equal the total area of the structure studied. In other words, we have to divide the plane of the structure in small pieces, each one belonging to one profile; this is called a ‘tessellation’, from a Latin word meaning small tile (tessella). Unlike tiles, the pieces that cover the plane do so without overlap or interstice. Consider two neighbouring profiles (Fig. 1), arbitrarily labelled 1 and 2; the portion of space that belongs to the first should be closer to it than to the other. If we draw a line at mid distance between the two profiles (in the configuration of Fig. 1), the plane located on the right side of the line belongs to profile 1, on the left side to profile 2. By repeating this procedure, a polygon is delineated around each neuronal profile that encloses the portion of the plane that is closer to that neuronal profile than to any other one. This polygon is called Voronoi or Dirichlet polygon (sometimes also ‘cell’ but, in our case, it could be confused with a biological cell). Voronoi tessellation has been extensively studied (Okabe et al., 1992) but has been only recently applied to the analysis of histological maps (Famiglieti, 1992;
Marcelpoil and Usson, 1992; Reymond et al., 1993; Seilhean et al., 1993; Duyckaerts et al., 1994; Grignon et al., 1998; Savy et al., 1999). We can now consider that each neuronal profile has a density of 1/(area of the polygon) (Duyckaerts et al., 1994). Several advantages may be gained by using tessellation rather than the usual counting of profiles inside boundaries. Two of them may be emphasised: 1. A mean neuronal density may be immediately qualified by a confidence interval. To assess the density of a group of neurones, one indeed needs first to calculate the mean area of the polygons assigned to each neurone. The means of n polygon areas have a Gaussian distribution when n is sufficiently high by virtue of the central limit theorem and may thus be associated with a symmetrical confidence interval. The standard error of the mean (S.E.M.) is the ratio of S.D./square root (n), n being the number of neurones taken into account in the calculation. As usual, if n is sufficiently large, there is a 0.95 probability that the actual mean is included in the range: estimated mean− 1.96 S.E.M.; estimated mean+ 1.96 SEM (practically 2 S.E.M.) This may be now converted into density values: 1/(mean polygon area) with confidence interval 1/(mean +2 S.E.M.) and 1/(mean − 2 S.E.M.). This interval is not symmetrical any more; the mean density is closer to the lower bound than to the upper bound of the interval; for instance, if the mean area of the polygons is 10 arbitrary units (au) with a confidence interval of 1 au, the confidence interval will be (10− 2) and (10+2) au. In terms of density, the mean is 1/10 and the lower bound is 1/12, i.e. closer to 1/10 than the upper bound (1/8). 2. The Voronoi tessellation is a way of determining the nearest neighbours of a given neurone. This simplifies the topographical analysis of a map. Simple techniques may then help to clarify the underlying structure and its pathology (clusters, boundaries between regions with different numerical densities, alterations of the laminar or columnar arrangement).
2. Practical steps
2.1. Data acquisition Fig. 1. How to draw a Voronoi polygon. Six points are represented. The Voronoi polygon of point numbered 1 in ‘a’ is serially built. A line located at mid-distance between 1 and 2 and perpendicular to line 1 and 2 is first drawn. The same procedure is performed in b. The two lines, just drawn, ends at their intersection (c) and a third side is now drawn. In d and e, the polygon is completed. The lines that were drawn between point 1 and each one of its neighbours are erased in e where the Dirichlet polygon is shown. The point is not necessarily in the centre of the polygon. By definition, there is only one point per polygon.
To perform the tessellation, the contour of a region of interest has to be drawn. The region must include the structures — such as layers or subnuclei — which are to be analysed. The position of each profile has to be mapped as a point within the contour. Several systems (Zeiss IBAS®, Biocom®) are commercially available to acquire the data (Fig. 2a). The microscopic slide is viewed on a video screen. The location of the moving
C. Duyckaerts, G. Godefroy / Journal of Chemical Neuroanatomy 20 (2000) 83–92
85
profile is visible on the first section and invisible on the second. This ‘di-sector’ procedure is treated at length in other articles of this issue. The optical disector consists in counting only those cells, which are visible in the upper part of the slide and become invisible when the focus is made on the lower part of the slide. Only high power objectives with a high numerical aperture are to be used since they have a small focal depth. It is possible but time-consuming to draw maps taking into consideration only the cells that are selected by an optical disector. Fig. 2. Equipment used to acquire the data. The image is grabbed by a video camera plugged in the phototube of a microscope. It is transmitted via the computer (E) to the video-screen (D). Linear transducers (B) and (C) keep track of the location of the moving stage of the microscope. The X and Y co-ordinates are fed to the computer. The computer mouse (F) is used to draw the regions of interest and to point the cells.
stage of the microscope is determined by linear transducers added to the microscope to record movements along the x- and y-axes. Another possibility is to use the co-ordinates of an automatic moving stage. The position of the mouse relative to the screen is also known. The final co-ordinates of an object pointed on the screen is the sum of the co-ordinates of the moving stage and of the mouse. The mouse is used to draw the borders of the region of interest and to point the profiles. The acquisition could theoretically be performed automatically but the present state of the art in image recognition is generally not sufficient to ensure an acceptable quality in the identification of the profiles. The profiles that have to be considered in the mapping have to be carefully identified. Obviously, they should first be recognised as neuronal or glial, possibly with the help of immunohistochemistry but this is not sufficient. When calculating a numerical density each cell, big or small, should count for 1 unit. In other words, each cell should be reduced to a point in the containing volume. If this volume was serially cut, each point should appear only once on the sections. The profile of the sectioned cells appears, on the contrary, in several slices, their number depending on the size of the cell. If any profile is counted as one unit, the large cells will be counted more often than the small ones and the evaluation will be biased. The bias depends on the size of the object, which is counted relative to the section thickness. Counting the nucleoli rather than the nuclei reduces the bias (see (Duyckaerts et al., 1989) for an estimate of the error). Other so-called ‘stereological’ methods avoiding the bias have been proposed; counting the inferior (or superior) pole of the cells completely avoids it (Sterio, 1984). Two contiguous sections are examined, the cell is taken into account only if its
2.2. Data processing The X and Y co-ordinates of the cells and the contour of the region are used to obtain the tessellation. Various algorithms are available, among which is Fortune’s described in detail in de Berg et al. (1997). We have adapted the one proposed by Green and Sibson (1978). Briefly, each cell is represented by a point defined by its x- and y-co-ordinates located within a closed border. Each one of these points is associated with the list of its neighbours in an ordered manner (for example, clockwise). The neighbours may either be another point or a border. Each point is included in only one polygon whose sides are located at mid-distance between two points and are perpendicular to the line drawn between them. The tessellation is obtained sequentially. At step one, only one point is introduced in the map. Its contiguity list contains only the borders. At step two, a second point is introduced. The contiguity list of point 1 is modified; it now includes point 2 and only some of the borders, while the contiguity list of point 2 contains point 1 and the other borders. The tessellation is completed by repeating the procedure adding one point at a time. A new point, when added to the map, will necessarily be in the contiguity list of the point that is nearest to it. This property — that may be designated as ‘neighbour of nearest point’ — was used in the Green and Sibson algorithm to initiate the changes of the contiguity lists. It is valid when the tessellation is performed in convex contours, i.e. in contours in which any line joining two points never crosses a border (Fig. 3). When the tessellation is performed within contours of any shape, the ‘neighbour of nearest point’ property is not valid any more but may be replaced by an alternative property (Fig. 3). When a new point is added, it necessarily falls within one polygon and will necessarily be the neighbour of the point to which the polygon belongs. This property may be designated as ‘neighbour if within polygon’. Rather than find the nearest neighbour, one then has to find in which polygon the new point is included. With the ‘neighbour if within polygon’ property, it is possible to tessellate almost all contours (see Fig. 3).
86
C. Duyckaerts, G. Godefroy / Journal of Chemical Neuroanatomy 20 (2000) 83–92
Fig. 3. Various types of contours in which the points may be located. Panel A shows that any line is completely included in a convex contour which may be intersected only twice (arrows). When the contour is not convex (on the right), a line may intersect the contour more than twice (arrows) and is not fully included in its borders. In such a non-convex region, the nearest point is not necessarily a neighbour, point 2 is closer to 1 but 3 is its neighbour.
Practically, the points are acquired with any commercially available system that might finally provide a file where the co-ordinates of the points are listed, as well as the co-ordinates of the border. A first program is dedicated to the computation of the tessellation. It provides a new file where the points are listed with their contiguity list and which is finally used in a second program that allows the treatment and the analysis of the data. The two programmes are freely available on request1.
3. Analysis of the data Once the tessellation has been performed, the data may be analysed in several ways.
3.1. Calculation of the numerical density. Comparison of densities The numerical density of neuronal profiles is the ratio: number of profiles/area in which they stand.
Fig. 4. Relationship between the numerical density and the Voronoi polygons. Two points are included within a rectangular contour. The numerical density is 2/area of the rectangle. The area of the Voronoi polygon of point 1 is A1, of point 2, A2. The mean area of the two polygons is (A1 +A2)/2. The mean density, i.e. 1/(mean polygon area) is 2/(A1 + A2); A1+A2=area of the rectangle; 1/mean polygon area = 2/(A1 +A2)= 2/area of the rectangle. The tessellation method and the standard counting procedure yield identical results of numerical density when the global values are considered. 1 The interested reader should ask for the programs at the following e-mail address: [email protected].
When the complete structure is considered the mean densities as calculated directly or by Voronoi tessellation are identical. This is easy to understand in a simple configuration as the one shown in Fig. 4. A rectangle contains two points. The area of the polygon surrounding point 1 is A1, the area surrounding point 2 is A2. A1 +A2 is the total area At of the rectangle. The numerical density computed in the usual way is 2/At. If the tessellation is used, the mean polygon area is (A1+ A2)/2 and the mean density is the inverse of this value, i.e. 2/(A1 + A2)= 2/At. Individual densities cannot be averaged directly (density of point 1+ density of point 2)/2 (mean density of (point 1, 2). But mean density of (point 1, 2) = harmonic mean of density of point 1 and 2= 1/(mean area of (A1, A2)). However, even if the global mean density is considered there are some advantages in using the tessellation. Various colours may be attributed to the Voronoi polygons according to their sizes; in this way, a coloured density map is drawn, where regions of high densities are readily visible (Fig. 5). Moreover, as already alluded to, the values of the polygon areas may be used to compute not only a mean value but also the standard deviation and the standard error. This permits the evaluation of a confidence interval attributed to the mean and provides information concerning the variability of the numerical density inside the region under study. It is intuitively clear that a mean value of density is more reliable in a structure where it does not vary much, i.e. in a region where the neurones are regularly spaced than in a structure comprising several clusters. The confidence interval allows an objective comparison of two numerical densities; if the two confidence intervals do not overlap, the two values are statistically different (as shown in Fig. 6). When the confidence intervals overlap, the test is unable to prove that the values are different. Tessellation allows selecting any group of neurones to determine its numerical density. This permits the analysis of the laminar or columnar arrangements of the neurones. It was done, for instance, in a study concerning the cerebral cortex in AIDS (Seilhean et al., 1993). No statistical difference could be found between patients with and patients without cognitive impairment. The columnar and laminar structure was analysed and found to be similar. The density of the neurones that are located close to the border of the map mostly depends on the way the region has been delineated. Their density is too small, when the borders have been drawn at large or too high when they are located at too short a distance of the peripheral neurones. Tessellation allows a rapid identification of the neurones located at the periphery (they have a border in their contiguity list). They can be removed before the density is calculated in order to decrease the border effects on the final value (Fig. 7).
C. Duyckaerts, G. Godefroy / Journal of Chemical Neuroanatomy 20 (2000) 83–92
87
Fig. 5. Density map of the neurofibrillary tangles in the human amygdala at an early stage of Alzheimer’s disease. Serial sections through the right amygdala were performed at 1 mm intervals. The anterior part is at the top of the page. The neurofibrillary tangles were mapped after tau immunolabelling. The colour scale follows a logarithmic progression and indicates the density of profiles of neurofibrillary tangles per mm2. A high density of neurofibrillary tangles was observed in the entorhinal cortex and in the basal and accessory basal nuclei of the amygdala, that appear in red. Data from Duyckaerts et al. (1999). Bar = 5 mm.
88
C. Duyckaerts, G. Godefroy / Journal of Chemical Neuroanatomy 20 (2000) 83–92
Fig. 6. BRDU positive neurones in the cortex of developing mice (courtesy, Dr Cristina Costa). At that stage, the density of positive cells is clearly higher in the deeper layers of the cortex. To discriminate the two large areas of different density, an automatic search of clusters was performed. It allowed the selection of two clusters with contrasted (and statistically different) densities within the same region.
3.2. Stretching/shrinking of a tessellated map Most of the time, atrophy accompanies neuronal loss, which it masks; shrinkage, indeed brings numerical density values back to normal. It is possible to shrink or to stretch a tessellation since the change of scale does not alter the contiguities. The coefficient of variation of the polygon areas is not modified, contrarily to the numerical densities. A ‘stretching’ of the cortical samples was, for instance, performed to compare the numerical densities observed in the isocortex of patients affected by Alzheimer’s disease; the procedure has been used to standardise the cortical thickness and allowed to uncover the neuronal loss (Grignon et al., 1998) and Fig. 8.
3.3. Characterisation of the distribution (random, regular, clustered) Three main types of distribution of points may be described, random, clustered or regular. In the Poisson distribution, the density of the points is, on an average, the same whatever the scale at which the map is examined and wherever it is observed. The chance of finding n points in a region of area A follows Poisson’s law, Pn,a =((da)ne − da)/n! where Pn,a is the probability of finding n points in an area a, d is the mean density of points and n! is n factorial: 1×2 × 3 ×···× n. Poisson’s distribution in the plane is simulated by taking random numbers for the x- and y-co-ordinates. The distribution of the points may depart from Pois-
son’s in two ways; they may be either too regularly spaced or clustered (Diggle, 1983). Voronoi tessellation may help to determine which spatial arrangement is the most plausible. When they follow the Poisson distribution, the Voronoi polygons are of variable areas, some large and other small; the variability of polygon areas is easily assessed by their variance. The coefficient of variation, CV (standard deviation of the polygon areas/ mean) allows to express the variability in a scale independent manner. When the distribution is regular the CV is low. When clusters are present, the polygons are small in the clusters and high between them. The association of small and large polygons explains why the CV is high. Periodic structures could lead to high CV values although exhibiting, in some sense, regularity; when the CV is high, however, it is a clustered distribution of points that is repeated periodically. To our knowledge the distribution of the polygon areas in a Poisson point process has not been directly calculated. To determine which was the range of variation of the CV of polygon areas, we had to resort to a Monte Carlo simulation. We found that it ranged from 33 to 64% with a probability of 0.99 (Duyckaerts et al., 1994). The mean value was 52.9% (identical with the value obtained by (Mo¨ller, 1994). It is worth noticing here that the CV of the polyhedron volumes was found to be 42.3 in a 3-D Poisson process, a result also obtained by simulation (Mo¨ller, 1994). The CV of the polygon areas was used to study the distribution of the tyrosine hydroxylase-positive neurones in the retina of weaver mice. It could be shown
C. Duyckaerts, G. Godefroy / Journal of Chemical Neuroanatomy 20 (2000) 83–92
89
Fig. 7. Decreasing the border effects by erosion. The picture was taken from a diagram of Nieuwenhuys et al. (1988). It represents the inferior olive and the neurones it contains. The two lower panels are the corresponding density maps. Black polygons are the smallest ones and are located in areas of high density, while the white ones are the largest, located in areas of low density. The size of the peripheral polygons is directly influenced by the drawing of the border. They appear to be filled with a light shade of grey on the left lower panel because the border was drawn too far away from the neurones. To avoid this artefact, all the peripheral polygons, i.e. those touching a border, have been removed in the lower right panel. In this way, a so-called ‘erosion’ was obtained. The size of the polygons that are left is solely determined by the proximity of the cells themselves; this increases the precision of the evaluation. Before erosion, the density was 436 points/arbitrary unit of area (confidence interval, 419–455), it reached 534 (confidence interval, 514–556) after erosion, a significantly different value.
that those neurones were less homogeneously distributed in the mutant mice; they were more clustered (Fig. 9; Savy et al., 1999).
3.4. Automatic clustering The tessellation leads to a natural and objective definition of a cluster of neurones: clustered neurones are contiguous and exhibit similar densities, i.e. their polygons have a similar surface area. When several clusters are present in the map, neurones between them have polygons of a larger area than the neurones inside them (Fig. 10). A small polygon, in a map, necessarily means that several neurones are close together, i.e. are included in a cluster. The smallest polygon of the map is, therefore, necessarily included in a cluster, if a cluster is to be present. This smallest polygon is first to be found. The contiguous polygons are then sequentially examined. They are left in the cluster only if their area is ‘similar’ to the one of the smallest polygon. The similarity is evaluated by calculating the difference in the surface areas of the two polygons or their ratio. The cluster is completed when there is no polygon left in the map, fulfilling the two conditions of contiguity and of area similarity. The two groups of neurones isolated in Fig. 6 were identified by an automatic clustering procedure.
4. From 2-D to 3-D The tessellation that we have used up to now is 2-D. Is it possible to use the information obtained in the plane to reconstruct the biological structure in 3-D? This is the subject of stereology, a word coined on 11 May 1961 when the International Society for Stereology was founded. ‘‘Stereology is the three-dimensional interpretation of flat images, such as sections and projections, by criteria of geometric probability’’ (Helias and Hyde, 1983)2. It is indeed possible to propose 3-D interpretation of a 2-D tessellation. It should first be considered that a 2-D tessellation is not a slice of a 3-D tessellation. Such a slice would indeed include polygons without point, the point being located above or under the slice. This remark may then lead to the next question: can the volume of a Voronoi polyhedron be extrapolated from the area of a Voronoi polygon? This is obviously not possible for individual points but it may be demonstrated that, on average, the volume of the polyhedron is equal to the area of the polygon times the thickness of the section. This might be intuitively 2 The reader, interested in etymology, may notice that the Greek word ‘stereos’ is an adjective meaning ‘solid’, used by Euclid in the geometrical sense of ‘volume’. It might also be worth noticing that the word has been used (for example by Plato) in the sense of ‘harsh’ or ‘uncouth’ — a quality that is fortunately not always to be applied to the statements or reviews produced by enthusiast stereologists.
90
C. Duyckaerts, G. Godefroy / Journal of Chemical Neuroanatomy 20 (2000) 83–92
understood if we consider a slice of a 3-D tessellation with flat top and bottom. We take in this slice the polygons, which contain a point and complete their polyhedron above or under the slice. By so doing, the top and the bottom of the slice become ‘hilly’. Now we consider the polygons, which do not contain a point. Their polyhedron is completed and taken off the slice. This would lead to the formation of ‘valleys’. We can now sum up the volumes of the polyhedra which are above or under the slice and call the result ‘volume of the hills’ and sum the volumes of the polyhedra which, on the contrary, entered the slice and call the result ‘volume of the valleys’. In a Poisson distribution, volume of valleys = volume of the hills and hence, as a mean, area of the polygons ×section thickness =volume of the polyhedra. This has been formerly demonstrated by one of us (GG) in (Duyckaerts et al., 1994). It should also be mentioned here that, in a Poisson point process, the coefficient of variation of the volumes of the polyhedra (in 3-D) is different from the coefficient of variation of the areas of the polygons (in 2-D). The value is 0.423 in 3-D (Mo¨ller, 1994) and 0.529 in 2-D (Duyckaerts et al., 1994; Mo¨ller, 1994).
Fig. 8. Stretching/shrinking a map (courtesy, Dr Yves Grignon). A and B are maps of an isocortical area (area 40) from an intellectually normal case and from a case severely affected by Alzheimer disease. Similar number of neurones were recorded in the two cases. The shade of grey is proportional to the density of the neuronal profiles (see scale on the left side of the picture). In the upper panels, the density appears similar in A and B, but B is atrophic: the height of the cortex is smaller. In C and D, the samples have been ‘stretched’ to obtain the same height as measured on the left side of the map. The stretching of map D uncloses the deficit in the density of neuronal profiles.
5. Further developments It is theoretically possible to acquire the co-ordinates of the cells in 3-D. A linear transducer affixed to the microscope moving stage could determine its position with reference to the z-axis. An objective with a high numerical aperture would have to be used to reduce the focal depth as much as possible. A cell would have to be pointed only if in focus (possibly using a disector type procedure); the x-, y- and z- co-ordinates could then be sent to the computer. Confocal laser microscopy could also allow to acquire the co-ordinates of the structures in 3-D. Algorithms of 3-D tessellation have been devised and are in use in image analysis (Bertin et al., 1993). To the best of our knowledge, they have not been applied yet to the evaluation of cellular density in 3-D.
6. Relation with other methods of statistical analysis of spatial point patterns The various tests using the nearest neighbour distance (Wa¨ssle and Riemann, 1978) rely on the comparison between an observed and the theoretical distribution, which is a function of the observed mean density of points. Only a subset of the significant distances (the short ones) is taken into account. This is the main weakness of the method: it actually amounts to ignoring all the sides of the tessellation polygons but one. Other methods have been used to describe the spatial arrangement of ‘events’ (in this context: neurones) in the plane (in this context the microscopic section) or in the volume. Some of them, such as analyses based on covariance and pair correlation functions, have been named ‘second-order stereology’ (Cruz-Orive, 1989) and have given valuable results in histology (Mayhew, 1999) and pathology (Mattfeldt et al., 1993). Diggle (1983) lists various descriptors of spatial arrangement such as inter-event distances (i.e. the distance between the points), nearest neighbour distances or point to nearest event distances. Several functions have also been devised: K(t)= density − 1 E (number of further events within distance t of an arbitrary event), G(Y)= P (distance from an arbitrary event to the nearest other event is at most Y), F(x)=P (distance from an arbitrary spot to the nearest event is at most x) where E ( ) stands for mean of ( ) and P ( ) for probability that ( ). In each case the distribution of the empirical data are to be compared with the results of simulations using a Poisson point process. Diggle proposes three sets of points to be used to determine the efficiency of the various methods in discriminating random, regular or clustered distributions. We used those sets to test the tessellation method. The CV of polygon
C. Duyckaerts, G. Godefroy / Journal of Chemical Neuroanatomy 20 (2000) 83–92
91
Fig. 9. Comparing the distribution of tyrosine hydroxylase-positive neurones in a normal mouse and in a weaver mutant at post-natal day 10 (courtesy Claudine Savy). The tyrosine hydroxylase positive neurones are more homogeneously distributed in the control mouse (on the left) than in the mutant (on the right). Each point is a neurone. The scale (in number of neurones per mm2) follows a logarithmic progression to cover the whole range of values. There are more dark polygons on the right, an evidence of clustering. The CV was 45% on the left, 51% on the right. For further information see Savy et al. (1999).
geneous Poisson process, contagious distribution, Poisson cluster processes etc…. Its weaknesses are the difficulty to establish a spatial model that fits the experimental data and the necessity of a time-consuming Monte-Carlo procedure to establish the ranges of the values obtained with the model. Tessellation provides robust techniques to visualise and to analyse in a quantitative manner the distribution of neurones, or, for that matter, of any cellular population. It should complement the currently used stereological methods, which do not allow apprehending maps of large populations of neurones.
Acknowledgements Fig. 10. Clusters. This is a diagram of three clusters of points to show that the polygons are small within the clusters (small arrows) and large between (long arrows).
areas permitted to discriminate them readily; the random distribution had a CV of 57% (i.e. comprised between 33 and 64%); the clustered one of 92% (i.e. \ 64%) and the regular one 29% (i.e. B 33%). The tessellation method is thus efficient in categorising the data into one of those three broad types of arrangement. The strength of the methods devised by Diggle is the possibility to test the data not only against the hypothesis of complete spatial randomness but also against other theoretical distributions such as inhomo-
The help of Patricia Gaspar (INSERM U106) is greatly acknowledged. We are grateful to the many people who have now used tessellation and contributed to improve the technique among whom Jean-Jacques Hauw, Danielle Seilhean, Yves Grignon, Malika Bennecib, Claudine Savy.
References Bertin, E., Parazza, F., Chassery, J.M., 1993. Segmentation and measurement based on 3D Voronoi diagram: application to confocal microscopy. Comput. Med. Imaging Graph. 17, 175–182.
92
C. Duyckaerts, G. Godefroy / Journal of Chemical Neuroanatomy 20 (2000) 83–92
Cruz-Orive, L.M., 1989. Second-order stereology: estimation of second moment volume measures. Acta Stereol. 8, 641–646. de Berg, M., van Kreveld, M., Overmars, M., Schwarzkopf, O., 1997. Computational Geometry. Algorithms and Applications. Springer, Berlin. Diggle, P.J., 1983. Statistical Analysis of Spatial Point Patterns. Academic Press, London. Duyckaerts, C., Llamas, E., Delae`re, P., Miele, P., Hauw, J.-J., 1989. Neuronal loss and neuronal atrophy. Computer simulation in connection with Alzheimer’s disease. Brain Res. 504, 94–100. Duyckaerts, C., Godefroy, G., Hauw, J.-J., 1994. Evaluation of neuronal numerical density by Dirichlet tessellation. J. Neurosci. Methods 51, 47 – 69. Duyckaerts, C., Sazdovitch, V., Colle, M.A., Hauw, J.J., 1999. Le`sions du noyau amygdalien dans les stades I et II de la maladie d’Alzheimer. In: Michel, B.F., Touchon, J., Pancrazi, M.P. (Eds.), Affect Amygdale Alzheimer. Solal, Marseille, pp. 243–250. Famiglieti, E.V., 1992. Polyaxonal amacrine cells of rabbit retina: size and distribution of PA1 cells. J. Comp. Neurol. 316, 406– 421. Green, P.J., Sibson, R., 1978. Computing Dirichlet tessellation in the plane, Comput. J. 21, 68–173. Grignon, Y., Duyckaerts, C., Bennecib, M., Hauw, J.-J., 1998. Cytoarchitectonic alterations in the supramarginal gyrus of late onset Alzheimer’s disease. Acta Neuropathol. (Berlin) 95, 395 – 406. Helias, H., Hyde, D.M., 1983. A Guide to Practical Stereology. Karger, Basel. Marcelpoil, R., Usson, Y., 1992. Methods for the study of cellular sociology: Voronoi diagrams and parametrization of the spatial relationships. J. Theor. Biol. 154, 359–369.
.
Mattfeldt, T., Frey, H., Rose, C., 1993. Second-order stereology of benign and malignant alterations of the human mammary gland. J. Microsc. 171, 143 – 151. Mayhew, T.M., 1990. Quantitative description of the spatial arrangement of organelles in a polarised secretory epithelial cell: the salivary gland acinar cell. J. Anat. 194, 279 – 285. Mo¨ller, J., 1994. Lectures on Random Voronoi Tessellations. Springer, New York. Nieuwanhuys, R., Voogd, J., Van Huı¨jzen, C., 1988. The Human Central Nervous System: a Synopsis and Atlas, third ed. Springer, New York. Okabe, A., Boots, B., Sugihara, K., 1992. Spatial Tessellations. Concepts and Applications of Voronoi Diagrams. Wiley, Chichester. Reymond, E., Raphae¨l, M., Grimaud, M., Vincent, L., Binet, J.L., Meyer, F., 1993. Germinal center analysis with the tools of mathematical morphology on graph. Cytometry 14, 848–861. Savy, C., Martin-Martinelli, E., Simon, A., Duyckaerts, C., Verney, C., Adelbrecht, C., Raisman-Vozari, R., Nguyen-Legros, J., 1999. Altered development of dopaminergic cells in the retina of weaver mice. J. Comp. Neurol. 412, 656 – 668. Seilhean, D., Duyckaerts, C., Vazeux, R.M., Bolgert, F., Brunet, P., Katlama, C., Gentilini, M., Hauw, J.-J., 1993. HIV-1 associated cognitive/motor complex: absence of neuronal loss in the cerebral cortex. Neurology 43, 1329 – 1334. Sterio, D.C., 1984. The unbiased estimation of number and sizes of arbitrary particles using the disector. J. Microsc. 134, 127–136. Wa¨ssle, H., Riemann, H.J., 1978. The mosaic of nerve cells in the mammalian retina. Proc. R. Soc. London B Biol. Sci. 200, 441– 461.
Voronoi tessellation to study the numerical density and the spatial distribution of neurones Charles Duyckaerts a,*, Gilles Godefroy b a
Laboratoire de Neuropathologie R. Escourolle, Hoˆpital de La Salpeˆtrie`re, 47 Boule6ard de l’Hoˆpital, 75651 Paris Cedex 13, France b Equipe d’analyse, Uni6ersite´ Paris 6, Boıˆte 186, 4, place Jussieu, 75252 Paris Cedex 05, France
Abstract The conditions of regularity and isotropy, required by standard morphometric procedures, are generally not fulfilled in the central nervous system (CNS) where cells are distributed in a highly complex manner. The evaluation of the mean numerical density of neuronal or glial cells does not take into account the topographical heterogeneity and thereby misses the information that it contains. A local measurement of the density can be obtained by evaluating the ‘numerical density of one cell’, i.e. the ratio 1/(the volume that the cell occupies). This volume is the region of space that is closer to that cell than to any other. It has the shape of a polyhedron, called Voronoi (or Dirichlet) polyhedron. In 2-D, the Voronoi polyhedron is a polygon, the sides of which are located at mid-distance from the neighbouring cells. The Voronoi polygons are contiguous and their set fills the space without interstice or overlap, i.e. they perform a ‘tessellation’ that may yield a density map when the same colours are used to fill polygons of similar sizes. The use of Voronoi polygons allows computing the confidence interval of a mean numerical density that makes statistical comparisons possible. The tessellation also provides information concerning spatial distribution; the areas of the Voronoi polygons do not vary much when the cells are regularly distributed. On the contrary, small and large polygons are found when cellular clusters are present. The coefficient of variation of the polygon areas is an objective measurement of their variability and helps to define ‘regular’, ‘clustered’ and ‘random’ distributions. When cells are clustered, small polygons are contiguous and may be objectively identified by simple algorithms. Voronoi tessellations are easily performed in 2-D. On an average the area of a polygon times the thickness of the section equals the volume of the corresponding polyhedron. 3-D tessellations that are theoretically possible and for which algorithms have been published remain to be adapted to histological works. © 2000 Elsevier Science B.V. All rights reserved. Keywords: Dirichlet tessellation; Voronoi tessellation; Morphometry; Neuronal density; Neuronal loss
1. Introduction Consider a histological view of such complex structures as the thalamus or the cortex. The section of the cells (their ‘profiles’) are grouped in subnuclei or in layers where they appear to be variously packed. One wants to compare the same structure in various cases or, in the same case, to analyse the characteristics of different regions. Ideally, we would like to be able to state that a layer in the cortex or a subnucleus in the thalamus lacks neurones or, on the contrary, is more densely packed than its neighbour. A crude way to * Corresponding author. Tel.: + 33-1-42161891; fax: +33-14239828. E-mail address: [email protected] (C. Duyckaerts).
obtain this information is based on the simple counting of the cellular profiles in a region that has been delineated. Dividing the number of profiles by the area of the subregion provides a mean ‘numerical density’ usually a number of cellular profiles per square mm. Notwithstanding the stereological problems that are raised by such a procedure (which we will consider later), average density values miss significant information; if, for instance, the cortex has a density of 1000 neuronal profiles per mm2 rather than 1200, this provides no information on the local characteristics of the changes, for instance, is layer IV affected more than layer V? Assessing the profile density in smaller areas, for example, in the various layers of the cortex is possible but raises other difficulties; the limits of the neuronal groupings such as layers or subnuclei are often difficult to define objectively since their identifica-
0891-0618/00/$ - see front matter © 2000 Elsevier Science B.V. All rights reserved. PII: S 0 8 9 1 - 0 6 1 8 ( 0 0 ) 0 0 0 6 4 - 8
84
C. Duyckaerts, G. Godefroy / Journal of Chemical Neuroanatomy 20 (2000) 83–92
tion itself relies on a subjective evaluation of numerical density. The error that is made in drawing the borders has a very significant impact on the final value, especially, when the structure is small. To improve the measurement, one would like to be able to evaluate objectively the ‘local density’ of profiles, without having to delineate the layers or the subnuclei whose identification is based on slight differences in packing density of the cellular profiles. We would also like to be able to consider a small group of neurones and immediately tell what their numerical density is. How small can this group be? Could we find, so to speak, the numerical density of one neuron? Its ‘density’ would then be 1/the area in which it stands. What is this area? If each neuronal profile has its own area, their sum should equal the total area of the structure studied. In other words, we have to divide the plane of the structure in small pieces, each one belonging to one profile; this is called a ‘tessellation’, from a Latin word meaning small tile (tessella). Unlike tiles, the pieces that cover the plane do so without overlap or interstice. Consider two neighbouring profiles (Fig. 1), arbitrarily labelled 1 and 2; the portion of space that belongs to the first should be closer to it than to the other. If we draw a line at mid distance between the two profiles (in the configuration of Fig. 1), the plane located on the right side of the line belongs to profile 1, on the left side to profile 2. By repeating this procedure, a polygon is delineated around each neuronal profile that encloses the portion of the plane that is closer to that neuronal profile than to any other one. This polygon is called Voronoi or Dirichlet polygon (sometimes also ‘cell’ but, in our case, it could be confused with a biological cell). Voronoi tessellation has been extensively studied (Okabe et al., 1992) but has been only recently applied to the analysis of histological maps (Famiglieti, 1992;
Marcelpoil and Usson, 1992; Reymond et al., 1993; Seilhean et al., 1993; Duyckaerts et al., 1994; Grignon et al., 1998; Savy et al., 1999). We can now consider that each neuronal profile has a density of 1/(area of the polygon) (Duyckaerts et al., 1994). Several advantages may be gained by using tessellation rather than the usual counting of profiles inside boundaries. Two of them may be emphasised: 1. A mean neuronal density may be immediately qualified by a confidence interval. To assess the density of a group of neurones, one indeed needs first to calculate the mean area of the polygons assigned to each neurone. The means of n polygon areas have a Gaussian distribution when n is sufficiently high by virtue of the central limit theorem and may thus be associated with a symmetrical confidence interval. The standard error of the mean (S.E.M.) is the ratio of S.D./square root (n), n being the number of neurones taken into account in the calculation. As usual, if n is sufficiently large, there is a 0.95 probability that the actual mean is included in the range: estimated mean− 1.96 S.E.M.; estimated mean+ 1.96 SEM (practically 2 S.E.M.) This may be now converted into density values: 1/(mean polygon area) with confidence interval 1/(mean +2 S.E.M.) and 1/(mean − 2 S.E.M.). This interval is not symmetrical any more; the mean density is closer to the lower bound than to the upper bound of the interval; for instance, if the mean area of the polygons is 10 arbitrary units (au) with a confidence interval of 1 au, the confidence interval will be (10− 2) and (10+2) au. In terms of density, the mean is 1/10 and the lower bound is 1/12, i.e. closer to 1/10 than the upper bound (1/8). 2. The Voronoi tessellation is a way of determining the nearest neighbours of a given neurone. This simplifies the topographical analysis of a map. Simple techniques may then help to clarify the underlying structure and its pathology (clusters, boundaries between regions with different numerical densities, alterations of the laminar or columnar arrangement).
2. Practical steps
2.1. Data acquisition Fig. 1. How to draw a Voronoi polygon. Six points are represented. The Voronoi polygon of point numbered 1 in ‘a’ is serially built. A line located at mid-distance between 1 and 2 and perpendicular to line 1 and 2 is first drawn. The same procedure is performed in b. The two lines, just drawn, ends at their intersection (c) and a third side is now drawn. In d and e, the polygon is completed. The lines that were drawn between point 1 and each one of its neighbours are erased in e where the Dirichlet polygon is shown. The point is not necessarily in the centre of the polygon. By definition, there is only one point per polygon.
To perform the tessellation, the contour of a region of interest has to be drawn. The region must include the structures — such as layers or subnuclei — which are to be analysed. The position of each profile has to be mapped as a point within the contour. Several systems (Zeiss IBAS®, Biocom®) are commercially available to acquire the data (Fig. 2a). The microscopic slide is viewed on a video screen. The location of the moving
C. Duyckaerts, G. Godefroy / Journal of Chemical Neuroanatomy 20 (2000) 83–92
85
profile is visible on the first section and invisible on the second. This ‘di-sector’ procedure is treated at length in other articles of this issue. The optical disector consists in counting only those cells, which are visible in the upper part of the slide and become invisible when the focus is made on the lower part of the slide. Only high power objectives with a high numerical aperture are to be used since they have a small focal depth. It is possible but time-consuming to draw maps taking into consideration only the cells that are selected by an optical disector. Fig. 2. Equipment used to acquire the data. The image is grabbed by a video camera plugged in the phototube of a microscope. It is transmitted via the computer (E) to the video-screen (D). Linear transducers (B) and (C) keep track of the location of the moving stage of the microscope. The X and Y co-ordinates are fed to the computer. The computer mouse (F) is used to draw the regions of interest and to point the cells.
stage of the microscope is determined by linear transducers added to the microscope to record movements along the x- and y-axes. Another possibility is to use the co-ordinates of an automatic moving stage. The position of the mouse relative to the screen is also known. The final co-ordinates of an object pointed on the screen is the sum of the co-ordinates of the moving stage and of the mouse. The mouse is used to draw the borders of the region of interest and to point the profiles. The acquisition could theoretically be performed automatically but the present state of the art in image recognition is generally not sufficient to ensure an acceptable quality in the identification of the profiles. The profiles that have to be considered in the mapping have to be carefully identified. Obviously, they should first be recognised as neuronal or glial, possibly with the help of immunohistochemistry but this is not sufficient. When calculating a numerical density each cell, big or small, should count for 1 unit. In other words, each cell should be reduced to a point in the containing volume. If this volume was serially cut, each point should appear only once on the sections. The profile of the sectioned cells appears, on the contrary, in several slices, their number depending on the size of the cell. If any profile is counted as one unit, the large cells will be counted more often than the small ones and the evaluation will be biased. The bias depends on the size of the object, which is counted relative to the section thickness. Counting the nucleoli rather than the nuclei reduces the bias (see (Duyckaerts et al., 1989) for an estimate of the error). Other so-called ‘stereological’ methods avoiding the bias have been proposed; counting the inferior (or superior) pole of the cells completely avoids it (Sterio, 1984). Two contiguous sections are examined, the cell is taken into account only if its
2.2. Data processing The X and Y co-ordinates of the cells and the contour of the region are used to obtain the tessellation. Various algorithms are available, among which is Fortune’s described in detail in de Berg et al. (1997). We have adapted the one proposed by Green and Sibson (1978). Briefly, each cell is represented by a point defined by its x- and y-co-ordinates located within a closed border. Each one of these points is associated with the list of its neighbours in an ordered manner (for example, clockwise). The neighbours may either be another point or a border. Each point is included in only one polygon whose sides are located at mid-distance between two points and are perpendicular to the line drawn between them. The tessellation is obtained sequentially. At step one, only one point is introduced in the map. Its contiguity list contains only the borders. At step two, a second point is introduced. The contiguity list of point 1 is modified; it now includes point 2 and only some of the borders, while the contiguity list of point 2 contains point 1 and the other borders. The tessellation is completed by repeating the procedure adding one point at a time. A new point, when added to the map, will necessarily be in the contiguity list of the point that is nearest to it. This property — that may be designated as ‘neighbour of nearest point’ — was used in the Green and Sibson algorithm to initiate the changes of the contiguity lists. It is valid when the tessellation is performed in convex contours, i.e. in contours in which any line joining two points never crosses a border (Fig. 3). When the tessellation is performed within contours of any shape, the ‘neighbour of nearest point’ property is not valid any more but may be replaced by an alternative property (Fig. 3). When a new point is added, it necessarily falls within one polygon and will necessarily be the neighbour of the point to which the polygon belongs. This property may be designated as ‘neighbour if within polygon’. Rather than find the nearest neighbour, one then has to find in which polygon the new point is included. With the ‘neighbour if within polygon’ property, it is possible to tessellate almost all contours (see Fig. 3).
86
C. Duyckaerts, G. Godefroy / Journal of Chemical Neuroanatomy 20 (2000) 83–92
Fig. 3. Various types of contours in which the points may be located. Panel A shows that any line is completely included in a convex contour which may be intersected only twice (arrows). When the contour is not convex (on the right), a line may intersect the contour more than twice (arrows) and is not fully included in its borders. In such a non-convex region, the nearest point is not necessarily a neighbour, point 2 is closer to 1 but 3 is its neighbour.
Practically, the points are acquired with any commercially available system that might finally provide a file where the co-ordinates of the points are listed, as well as the co-ordinates of the border. A first program is dedicated to the computation of the tessellation. It provides a new file where the points are listed with their contiguity list and which is finally used in a second program that allows the treatment and the analysis of the data. The two programmes are freely available on request1.
3. Analysis of the data Once the tessellation has been performed, the data may be analysed in several ways.
3.1. Calculation of the numerical density. Comparison of densities The numerical density of neuronal profiles is the ratio: number of profiles/area in which they stand.
Fig. 4. Relationship between the numerical density and the Voronoi polygons. Two points are included within a rectangular contour. The numerical density is 2/area of the rectangle. The area of the Voronoi polygon of point 1 is A1, of point 2, A2. The mean area of the two polygons is (A1 +A2)/2. The mean density, i.e. 1/(mean polygon area) is 2/(A1 + A2); A1+A2=area of the rectangle; 1/mean polygon area = 2/(A1 +A2)= 2/area of the rectangle. The tessellation method and the standard counting procedure yield identical results of numerical density when the global values are considered. 1 The interested reader should ask for the programs at the following e-mail address: [email protected].
When the complete structure is considered the mean densities as calculated directly or by Voronoi tessellation are identical. This is easy to understand in a simple configuration as the one shown in Fig. 4. A rectangle contains two points. The area of the polygon surrounding point 1 is A1, the area surrounding point 2 is A2. A1 +A2 is the total area At of the rectangle. The numerical density computed in the usual way is 2/At. If the tessellation is used, the mean polygon area is (A1+ A2)/2 and the mean density is the inverse of this value, i.e. 2/(A1 + A2)= 2/At. Individual densities cannot be averaged directly (density of point 1+ density of point 2)/2 (mean density of (point 1, 2). But mean density of (point 1, 2) = harmonic mean of density of point 1 and 2= 1/(mean area of (A1, A2)). However, even if the global mean density is considered there are some advantages in using the tessellation. Various colours may be attributed to the Voronoi polygons according to their sizes; in this way, a coloured density map is drawn, where regions of high densities are readily visible (Fig. 5). Moreover, as already alluded to, the values of the polygon areas may be used to compute not only a mean value but also the standard deviation and the standard error. This permits the evaluation of a confidence interval attributed to the mean and provides information concerning the variability of the numerical density inside the region under study. It is intuitively clear that a mean value of density is more reliable in a structure where it does not vary much, i.e. in a region where the neurones are regularly spaced than in a structure comprising several clusters. The confidence interval allows an objective comparison of two numerical densities; if the two confidence intervals do not overlap, the two values are statistically different (as shown in Fig. 6). When the confidence intervals overlap, the test is unable to prove that the values are different. Tessellation allows selecting any group of neurones to determine its numerical density. This permits the analysis of the laminar or columnar arrangements of the neurones. It was done, for instance, in a study concerning the cerebral cortex in AIDS (Seilhean et al., 1993). No statistical difference could be found between patients with and patients without cognitive impairment. The columnar and laminar structure was analysed and found to be similar. The density of the neurones that are located close to the border of the map mostly depends on the way the region has been delineated. Their density is too small, when the borders have been drawn at large or too high when they are located at too short a distance of the peripheral neurones. Tessellation allows a rapid identification of the neurones located at the periphery (they have a border in their contiguity list). They can be removed before the density is calculated in order to decrease the border effects on the final value (Fig. 7).
C. Duyckaerts, G. Godefroy / Journal of Chemical Neuroanatomy 20 (2000) 83–92
87
Fig. 5. Density map of the neurofibrillary tangles in the human amygdala at an early stage of Alzheimer’s disease. Serial sections through the right amygdala were performed at 1 mm intervals. The anterior part is at the top of the page. The neurofibrillary tangles were mapped after tau immunolabelling. The colour scale follows a logarithmic progression and indicates the density of profiles of neurofibrillary tangles per mm2. A high density of neurofibrillary tangles was observed in the entorhinal cortex and in the basal and accessory basal nuclei of the amygdala, that appear in red. Data from Duyckaerts et al. (1999). Bar = 5 mm.
88
C. Duyckaerts, G. Godefroy / Journal of Chemical Neuroanatomy 20 (2000) 83–92
Fig. 6. BRDU positive neurones in the cortex of developing mice (courtesy, Dr Cristina Costa). At that stage, the density of positive cells is clearly higher in the deeper layers of the cortex. To discriminate the two large areas of different density, an automatic search of clusters was performed. It allowed the selection of two clusters with contrasted (and statistically different) densities within the same region.
3.2. Stretching/shrinking of a tessellated map Most of the time, atrophy accompanies neuronal loss, which it masks; shrinkage, indeed brings numerical density values back to normal. It is possible to shrink or to stretch a tessellation since the change of scale does not alter the contiguities. The coefficient of variation of the polygon areas is not modified, contrarily to the numerical densities. A ‘stretching’ of the cortical samples was, for instance, performed to compare the numerical densities observed in the isocortex of patients affected by Alzheimer’s disease; the procedure has been used to standardise the cortical thickness and allowed to uncover the neuronal loss (Grignon et al., 1998) and Fig. 8.
3.3. Characterisation of the distribution (random, regular, clustered) Three main types of distribution of points may be described, random, clustered or regular. In the Poisson distribution, the density of the points is, on an average, the same whatever the scale at which the map is examined and wherever it is observed. The chance of finding n points in a region of area A follows Poisson’s law, Pn,a =((da)ne − da)/n! where Pn,a is the probability of finding n points in an area a, d is the mean density of points and n! is n factorial: 1×2 × 3 ×···× n. Poisson’s distribution in the plane is simulated by taking random numbers for the x- and y-co-ordinates. The distribution of the points may depart from Pois-
son’s in two ways; they may be either too regularly spaced or clustered (Diggle, 1983). Voronoi tessellation may help to determine which spatial arrangement is the most plausible. When they follow the Poisson distribution, the Voronoi polygons are of variable areas, some large and other small; the variability of polygon areas is easily assessed by their variance. The coefficient of variation, CV (standard deviation of the polygon areas/ mean) allows to express the variability in a scale independent manner. When the distribution is regular the CV is low. When clusters are present, the polygons are small in the clusters and high between them. The association of small and large polygons explains why the CV is high. Periodic structures could lead to high CV values although exhibiting, in some sense, regularity; when the CV is high, however, it is a clustered distribution of points that is repeated periodically. To our knowledge the distribution of the polygon areas in a Poisson point process has not been directly calculated. To determine which was the range of variation of the CV of polygon areas, we had to resort to a Monte Carlo simulation. We found that it ranged from 33 to 64% with a probability of 0.99 (Duyckaerts et al., 1994). The mean value was 52.9% (identical with the value obtained by (Mo¨ller, 1994). It is worth noticing here that the CV of the polyhedron volumes was found to be 42.3 in a 3-D Poisson process, a result also obtained by simulation (Mo¨ller, 1994). The CV of the polygon areas was used to study the distribution of the tyrosine hydroxylase-positive neurones in the retina of weaver mice. It could be shown
C. Duyckaerts, G. Godefroy / Journal of Chemical Neuroanatomy 20 (2000) 83–92
89
Fig. 7. Decreasing the border effects by erosion. The picture was taken from a diagram of Nieuwenhuys et al. (1988). It represents the inferior olive and the neurones it contains. The two lower panels are the corresponding density maps. Black polygons are the smallest ones and are located in areas of high density, while the white ones are the largest, located in areas of low density. The size of the peripheral polygons is directly influenced by the drawing of the border. They appear to be filled with a light shade of grey on the left lower panel because the border was drawn too far away from the neurones. To avoid this artefact, all the peripheral polygons, i.e. those touching a border, have been removed in the lower right panel. In this way, a so-called ‘erosion’ was obtained. The size of the polygons that are left is solely determined by the proximity of the cells themselves; this increases the precision of the evaluation. Before erosion, the density was 436 points/arbitrary unit of area (confidence interval, 419–455), it reached 534 (confidence interval, 514–556) after erosion, a significantly different value.
that those neurones were less homogeneously distributed in the mutant mice; they were more clustered (Fig. 9; Savy et al., 1999).
3.4. Automatic clustering The tessellation leads to a natural and objective definition of a cluster of neurones: clustered neurones are contiguous and exhibit similar densities, i.e. their polygons have a similar surface area. When several clusters are present in the map, neurones between them have polygons of a larger area than the neurones inside them (Fig. 10). A small polygon, in a map, necessarily means that several neurones are close together, i.e. are included in a cluster. The smallest polygon of the map is, therefore, necessarily included in a cluster, if a cluster is to be present. This smallest polygon is first to be found. The contiguous polygons are then sequentially examined. They are left in the cluster only if their area is ‘similar’ to the one of the smallest polygon. The similarity is evaluated by calculating the difference in the surface areas of the two polygons or their ratio. The cluster is completed when there is no polygon left in the map, fulfilling the two conditions of contiguity and of area similarity. The two groups of neurones isolated in Fig. 6 were identified by an automatic clustering procedure.
4. From 2-D to 3-D The tessellation that we have used up to now is 2-D. Is it possible to use the information obtained in the plane to reconstruct the biological structure in 3-D? This is the subject of stereology, a word coined on 11 May 1961 when the International Society for Stereology was founded. ‘‘Stereology is the three-dimensional interpretation of flat images, such as sections and projections, by criteria of geometric probability’’ (Helias and Hyde, 1983)2. It is indeed possible to propose 3-D interpretation of a 2-D tessellation. It should first be considered that a 2-D tessellation is not a slice of a 3-D tessellation. Such a slice would indeed include polygons without point, the point being located above or under the slice. This remark may then lead to the next question: can the volume of a Voronoi polyhedron be extrapolated from the area of a Voronoi polygon? This is obviously not possible for individual points but it may be demonstrated that, on average, the volume of the polyhedron is equal to the area of the polygon times the thickness of the section. This might be intuitively 2 The reader, interested in etymology, may notice that the Greek word ‘stereos’ is an adjective meaning ‘solid’, used by Euclid in the geometrical sense of ‘volume’. It might also be worth noticing that the word has been used (for example by Plato) in the sense of ‘harsh’ or ‘uncouth’ — a quality that is fortunately not always to be applied to the statements or reviews produced by enthusiast stereologists.
90
C. Duyckaerts, G. Godefroy / Journal of Chemical Neuroanatomy 20 (2000) 83–92
understood if we consider a slice of a 3-D tessellation with flat top and bottom. We take in this slice the polygons, which contain a point and complete their polyhedron above or under the slice. By so doing, the top and the bottom of the slice become ‘hilly’. Now we consider the polygons, which do not contain a point. Their polyhedron is completed and taken off the slice. This would lead to the formation of ‘valleys’. We can now sum up the volumes of the polyhedra which are above or under the slice and call the result ‘volume of the hills’ and sum the volumes of the polyhedra which, on the contrary, entered the slice and call the result ‘volume of the valleys’. In a Poisson distribution, volume of valleys = volume of the hills and hence, as a mean, area of the polygons ×section thickness =volume of the polyhedra. This has been formerly demonstrated by one of us (GG) in (Duyckaerts et al., 1994). It should also be mentioned here that, in a Poisson point process, the coefficient of variation of the volumes of the polyhedra (in 3-D) is different from the coefficient of variation of the areas of the polygons (in 2-D). The value is 0.423 in 3-D (Mo¨ller, 1994) and 0.529 in 2-D (Duyckaerts et al., 1994; Mo¨ller, 1994).
Fig. 8. Stretching/shrinking a map (courtesy, Dr Yves Grignon). A and B are maps of an isocortical area (area 40) from an intellectually normal case and from a case severely affected by Alzheimer disease. Similar number of neurones were recorded in the two cases. The shade of grey is proportional to the density of the neuronal profiles (see scale on the left side of the picture). In the upper panels, the density appears similar in A and B, but B is atrophic: the height of the cortex is smaller. In C and D, the samples have been ‘stretched’ to obtain the same height as measured on the left side of the map. The stretching of map D uncloses the deficit in the density of neuronal profiles.
5. Further developments It is theoretically possible to acquire the co-ordinates of the cells in 3-D. A linear transducer affixed to the microscope moving stage could determine its position with reference to the z-axis. An objective with a high numerical aperture would have to be used to reduce the focal depth as much as possible. A cell would have to be pointed only if in focus (possibly using a disector type procedure); the x-, y- and z- co-ordinates could then be sent to the computer. Confocal laser microscopy could also allow to acquire the co-ordinates of the structures in 3-D. Algorithms of 3-D tessellation have been devised and are in use in image analysis (Bertin et al., 1993). To the best of our knowledge, they have not been applied yet to the evaluation of cellular density in 3-D.
6. Relation with other methods of statistical analysis of spatial point patterns The various tests using the nearest neighbour distance (Wa¨ssle and Riemann, 1978) rely on the comparison between an observed and the theoretical distribution, which is a function of the observed mean density of points. Only a subset of the significant distances (the short ones) is taken into account. This is the main weakness of the method: it actually amounts to ignoring all the sides of the tessellation polygons but one. Other methods have been used to describe the spatial arrangement of ‘events’ (in this context: neurones) in the plane (in this context the microscopic section) or in the volume. Some of them, such as analyses based on covariance and pair correlation functions, have been named ‘second-order stereology’ (Cruz-Orive, 1989) and have given valuable results in histology (Mayhew, 1999) and pathology (Mattfeldt et al., 1993). Diggle (1983) lists various descriptors of spatial arrangement such as inter-event distances (i.e. the distance between the points), nearest neighbour distances or point to nearest event distances. Several functions have also been devised: K(t)= density − 1 E (number of further events within distance t of an arbitrary event), G(Y)= P (distance from an arbitrary event to the nearest other event is at most Y), F(x)=P (distance from an arbitrary spot to the nearest event is at most x) where E ( ) stands for mean of ( ) and P ( ) for probability that ( ). In each case the distribution of the empirical data are to be compared with the results of simulations using a Poisson point process. Diggle proposes three sets of points to be used to determine the efficiency of the various methods in discriminating random, regular or clustered distributions. We used those sets to test the tessellation method. The CV of polygon
C. Duyckaerts, G. Godefroy / Journal of Chemical Neuroanatomy 20 (2000) 83–92
91
Fig. 9. Comparing the distribution of tyrosine hydroxylase-positive neurones in a normal mouse and in a weaver mutant at post-natal day 10 (courtesy Claudine Savy). The tyrosine hydroxylase positive neurones are more homogeneously distributed in the control mouse (on the left) than in the mutant (on the right). Each point is a neurone. The scale (in number of neurones per mm2) follows a logarithmic progression to cover the whole range of values. There are more dark polygons on the right, an evidence of clustering. The CV was 45% on the left, 51% on the right. For further information see Savy et al. (1999).
geneous Poisson process, contagious distribution, Poisson cluster processes etc…. Its weaknesses are the difficulty to establish a spatial model that fits the experimental data and the necessity of a time-consuming Monte-Carlo procedure to establish the ranges of the values obtained with the model. Tessellation provides robust techniques to visualise and to analyse in a quantitative manner the distribution of neurones, or, for that matter, of any cellular population. It should complement the currently used stereological methods, which do not allow apprehending maps of large populations of neurones.
Acknowledgements Fig. 10. Clusters. This is a diagram of three clusters of points to show that the polygons are small within the clusters (small arrows) and large between (long arrows).
areas permitted to discriminate them readily; the random distribution had a CV of 57% (i.e. comprised between 33 and 64%); the clustered one of 92% (i.e. \ 64%) and the regular one 29% (i.e. B 33%). The tessellation method is thus efficient in categorising the data into one of those three broad types of arrangement. The strength of the methods devised by Diggle is the possibility to test the data not only against the hypothesis of complete spatial randomness but also against other theoretical distributions such as inhomo-
The help of Patricia Gaspar (INSERM U106) is greatly acknowledged. We are grateful to the many people who have now used tessellation and contributed to improve the technique among whom Jean-Jacques Hauw, Danielle Seilhean, Yves Grignon, Malika Bennecib, Claudine Savy.
References Bertin, E., Parazza, F., Chassery, J.M., 1993. Segmentation and measurement based on 3D Voronoi diagram: application to confocal microscopy. Comput. Med. Imaging Graph. 17, 175–182.
92
C. Duyckaerts, G. Godefroy / Journal of Chemical Neuroanatomy 20 (2000) 83–92
Cruz-Orive, L.M., 1989. Second-order stereology: estimation of second moment volume measures. Acta Stereol. 8, 641–646. de Berg, M., van Kreveld, M., Overmars, M., Schwarzkopf, O., 1997. Computational Geometry. Algorithms and Applications. Springer, Berlin. Diggle, P.J., 1983. Statistical Analysis of Spatial Point Patterns. Academic Press, London. Duyckaerts, C., Llamas, E., Delae`re, P., Miele, P., Hauw, J.-J., 1989. Neuronal loss and neuronal atrophy. Computer simulation in connection with Alzheimer’s disease. Brain Res. 504, 94–100. Duyckaerts, C., Godefroy, G., Hauw, J.-J., 1994. Evaluation of neuronal numerical density by Dirichlet tessellation. J. Neurosci. Methods 51, 47 – 69. Duyckaerts, C., Sazdovitch, V., Colle, M.A., Hauw, J.J., 1999. Le`sions du noyau amygdalien dans les stades I et II de la maladie d’Alzheimer. In: Michel, B.F., Touchon, J., Pancrazi, M.P. (Eds.), Affect Amygdale Alzheimer. Solal, Marseille, pp. 243–250. Famiglieti, E.V., 1992. Polyaxonal amacrine cells of rabbit retina: size and distribution of PA1 cells. J. Comp. Neurol. 316, 406– 421. Green, P.J., Sibson, R., 1978. Computing Dirichlet tessellation in the plane, Comput. J. 21, 68–173. Grignon, Y., Duyckaerts, C., Bennecib, M., Hauw, J.-J., 1998. Cytoarchitectonic alterations in the supramarginal gyrus of late onset Alzheimer’s disease. Acta Neuropathol. (Berlin) 95, 395 – 406. Helias, H., Hyde, D.M., 1983. A Guide to Practical Stereology. Karger, Basel. Marcelpoil, R., Usson, Y., 1992. Methods for the study of cellular sociology: Voronoi diagrams and parametrization of the spatial relationships. J. Theor. Biol. 154, 359–369.
.
Mattfeldt, T., Frey, H., Rose, C., 1993. Second-order stereology of benign and malignant alterations of the human mammary gland. J. Microsc. 171, 143 – 151. Mayhew, T.M., 1990. Quantitative description of the spatial arrangement of organelles in a polarised secretory epithelial cell: the salivary gland acinar cell. J. Anat. 194, 279 – 285. Mo¨ller, J., 1994. Lectures on Random Voronoi Tessellations. Springer, New York. Nieuwanhuys, R., Voogd, J., Van Huı¨jzen, C., 1988. The Human Central Nervous System: a Synopsis and Atlas, third ed. Springer, New York. Okabe, A., Boots, B., Sugihara, K., 1992. Spatial Tessellations. Concepts and Applications of Voronoi Diagrams. Wiley, Chichester. Reymond, E., Raphae¨l, M., Grimaud, M., Vincent, L., Binet, J.L., Meyer, F., 1993. Germinal center analysis with the tools of mathematical morphology on graph. Cytometry 14, 848–861. Savy, C., Martin-Martinelli, E., Simon, A., Duyckaerts, C., Verney, C., Adelbrecht, C., Raisman-Vozari, R., Nguyen-Legros, J., 1999. Altered development of dopaminergic cells in the retina of weaver mice. J. Comp. Neurol. 412, 656 – 668. Seilhean, D., Duyckaerts, C., Vazeux, R.M., Bolgert, F., Brunet, P., Katlama, C., Gentilini, M., Hauw, J.-J., 1993. HIV-1 associated cognitive/motor complex: absence of neuronal loss in the cerebral cortex. Neurology 43, 1329 – 1334. Sterio, D.C., 1984. The unbiased estimation of number and sizes of arbitrary particles using the disector. J. Microsc. 134, 127–136. Wa¨ssle, H., Riemann, H.J., 1978. The mosaic of nerve cells in the mammalian retina. Proc. R. Soc. London B Biol. Sci. 200, 441– 461.