Journal of Microscopy, Vol. 256, Issue 1 2014, pp. 46–60

doi: 10.1111/jmi.12157

Received 17 April 2014; accepted 21 June 2014

Characterization of squamous cell carcinomas of the head and neck using methods of spatial statistics T. MATTFELDT∗ & F. FLEISCHER†

∗ Institute of Pathology, University of Ulm, Um, Germany

†Medical Data Services/Biostatistics, Boehringer-Ingelheim Pharma GmbH & Co. KG, Biberach, Germany

Key words. Pattern recognition, point process, spatial statistics, squamous cell carcinoma, stereology, tessellation.

Summary In the present study, 53 cases of squamous cell carcinomas of the head and neck were characterized by a quantitative histological texture analysis based on principles of spatial statistics. A planar tessellation of the epithelial tumour component was generated by a skeletonization algorithm. The size distribution of the virtual cells of this planar tessellation, and the size distribution of the profiles of the tumour cell nuclei were estimated in terms of area and boundary length. The intensity, the reduced second moment function (K -function) and the pair correlation function of the point process of the centroids of the profiles of the tumour cell nuclei were also estimated. For both purposes, it is necessary to correct for edge effects, which we consider in this paper in some detail. Specifically, the point patterns of the tumour cell nuclei were considered as realizations of a point process, where the points exist only in the epithelial tumour component (the permitted phase) and not in the stroma (the forbidden phase). The methods allow to characterize each individual tumour by a series of summary statistics. The total set of cases was then partitioned into two groups: 19 cases without lymph node metastases (pN0), and 34 nodal positive cases (pN1 or pN2). Statistical analysis showed no significant differences between the intensities, the mean K -functions and the mean pair correlation functions of the tumour cell nucleus profiles of the two groups. However, there were some significant differences between the sizes of the virtual cells and of the nucleus profiles of the nodal negative cases as compared to the nodal positive cases. In a logistic regression analysis, one of the quantitative nuclear size variables (mean nuclear area) was found to be a significant predictor of lymph node metastasis, in addition to tumour stage. The study shows the potential of methods of spatial statistics for objective quantitative grading of squamous cell carcinomas of the head and neck, and provides an example for modelling histological Correspondence to: Prof. Dr. T. Mattfeldt, Institute of Pathology, Oberer Eselsberg M23, D-89081 Ulm, Germany. Tel.: +49-731-500-56322; fax: +49-731-50056384; e-mail: [email protected]

point patterns as realizations of planar point processes occupying a reference phase which is only a partial component of the total tissue.

Introduction Squamous cell carcinoma is the most frequent malignant tumour of the head and neck region. This category includes squamous cell carcinomas of the oropharynx, hypopharnyx, oral cavity, and larynx (Barnes et al., 2005). Major options for therapy are surgical resection, irradiation, and chemotherapy. Malignancy grading and assessment of prognostic factors are important for therapy planning. Tumour and lymph node status are important prognostic parameters. Moreover, histopathological grading is performed routinely in order to classify the tumours into grades of malignancy. Usually, it is tried to classify the tumours in a system consisting of three grades, where grade I means a high differentiation, grade II means an intermediate differentiation, and grade III means a low grade of differentiation (Barnes et al., 2005). A problem in the histopathological grading of squamous cell carcinomas is the local heterogeneity of the tumour tissue. Another difficulty is the subjective nature of the histopathological grading, which may lead to poor reproducibility of grading results between observers and also to intraobserver variability. Various attempts have been made at an objective and quantitative evaluation of the histological image of squamous cell carcinomas. Interest has been focused largely on the tumour cell nuclei. The main ideas in this direction may be considered under three headings: nuclear morphometry, stereology and evaluation of nuclear texture. With nuclear morphometry we mean quantitation of geometrical properties of profiles of tumour cell nuclei on histological sections, such as areas, boundary lengths, shape factors of nucleus profiles etc., see, e.g. (Bjerregaard et al., 1993; Narvaez et al., 1997; Okon et al., 1998; Dobros et al., 1999; Weyn et al., 2000; Maiolino et al., 2002; Sekine et al., 2003; Lavie et al., 2006; Martano et al.,  C 2014 The Authors C 2014 Royal Microscopical Society Journal of Microscopy 

SPATIAL STATISTICS OF SQUAMOUS CELL CARCINOMAS

2006; Pektas et al., 2006; Huang et al., 2008; Natarajan et al., 2010; Glazer et al., 2011). Stereological methods try to extrapolate three-dimensional properties of constituents of the tumour tissue from the planar sections to space by the help of mathematical methods. Here, we mention, e.g. the estimation of v¯ v , the volume-weighted mean volume of nuclei (Bundgaard et al., 1992; Sorensen et al., 1992; Bjerregaard et al., 1993) and other stereological parameters (Liu et al. 2012). Evaluation of nuclear texture means a quantitation of the distribution of the chromatin inside the tumour cell nucleus profiles at high magnification (Weyn et al., 2000; Dreyer et al., 2001; Natarajan et al., 2010). Most researchers agree that the aforementioned methods provide significant indicators of prognosis, but the relative value of the objective methods as compared to routine grading is debated. In this paper, we suggest another approach based on the following ideas. The methods quoted above provide numerical and objective data on nuclei, but they do not give information on the tumour architecture, i.e. the pattern (arrangement) of the nuclei and the cells within the tumour tissue. However, malignancy grading relies not only on the tumour cell nuclei, but also heavily on the degree of tissue differentiation, i.e. the tumour texture, and on the tumours cells as a whole (Wahi et al., 1971; Barnes et al., 2005). We suggest to quantify the tumour structure using methods of spatial statistics of point processes and tessellations (Illian et al., 2008; Chiu et al., 2013), which also reflect texture and cell size. These are well established, e.g. in geography, general biology and telecommunications (Gloaguen et al., 2006; Fleischer et al., 2008), but have only scarcely been applied to the study of squamous cell carcinomas of the head and neck. The primary aim of the study was to establish our technique methodologically. In practice, we used a set of 53 squamous cell carcinomas of different size, stage and locations, 34 of which had lymph node metastases, and 19 of which were free of lymph node metastases. Estimates of various extracted quantities were compared between nodal-negative and nodal-positive cases. The influence of input variables on the nodal state was evaluated by logistic regression analysis. Materials and methods Cases, histopathology and image acquisition From the files of the Department of Pathology of the University of Ulm, 53 cases of squamous cell carcinomas of the head and neck were recruited (19 pN0, 15 pN1 and 19 pN2 cases). Routine sections stained with Haematoxylin and Eosin were used for the study. At a magnification of 16× at the level of the objective, the images were acquired using a Zeiss light microscope attached to a Kontron IBAS image analysis system. The epithelial phase, the stromal phase and the tumour cell nucleus profiles were segmented interactively (see Figs. 1A, B). Two randomly selected visual fields per case were evaluated.  C 2014 The Authors C 2014 Royal Microscopical Society, 256, 46–60 Journal of Microscopy 

47

With respect to image analysis data are provided as greyscale valued images, where three different greyscale values are assumed (see Figs. 1A, B): a black phase describing the stroma, a white phase describing the continuum of tumour cell tissue, and a grey phase describing the cell nuclei. Image sizes are 512 × 512 pixels which corresponds to 250μm × 250μm at the level of the specimen. Spatial statistics General principles of planar tessellations. A planar tessellation is a space-filling partition of a planar structure into contiguous cells, usually based on a random mechanism (Chiu et al., 2013). A planar tessellation may visually be conceived as a planar mosaic. In the present application, we consider virtual tessellations generated by computer in the phase occupied by tumour cells. On routine sections, one sees only nucleus profiles inside a continuum, where boundaries between neighbouring cells cannot be delineated exactly. The centroids of the nucleus profiles, which can be well delineated, are taken as centres of virtual units which are constructed by image analysis, as described delow. The individual units of a planar tessellation are called ‘cells’. These cells are two-dimensional and virtual, whereas the true tumour cells of the carcinoma are real and three-dimensional. Hence, the tessellation cells are mathematical constructs and must not be confounded with the biological tumour cells. In this approach, we shall be content with a quantitation of the planar pattern of the nuclei and tessellation cells, as we see it in the sections. We use the estimated quantities only for the purpose of prediction and classification, and not for the elucidation of the 3D properties. Hence a ’nonstereological’ perspective appears operationally justified here. Image analysis. The aim of image analysis is to partition the epithelial continuum into separate cells according to the following rules:

r r r r

Each cell should contain one (and only one) nucleus. A nucleus should be contained in one (and only one) cell. Stromal tissue should not be contained in a cell. Cells should be constructed according to a nearest neighbour principle, i.e. a pixel belonging to the continuum of cell tissue should belong to the cell where the distance from the pixel to the associated cell nuclei is minimal.

In order to achieve a partition that is in accordance with the rules stated above, skeletonization algorithms have proven to be very useful tools. Skeletonization by morphological operators. The aim of skeletonization algorithms is to replace given structures or objects by thin skeletons, ideally of width one pixel, while preserving

48

T. MATTFELDT AND F. FLEISCHER

Fig. 1. (A) Original image displaying tumor tissue with epithelial component (darker areas) and stroma component (lighter areas). Haematoxylin-Eosin stain. (B) Triphasic image arising from upper panel after segmentation into three phases: black – stroma; grey – tumor cell nuclei; white – tumour tissue without nuclei.

the main shape characteristics and the homotopy, that means the number of objects and holes. In the following we use skeletonization by morphological operators (Serra, 1988; J¨ahne, 2001; Soille, 2003). Here 8 structural elements M1 , . . . , M8 are used to transform a given binary image into a skeleton structure. These structural elements are given by the following matrices ⎛ ⎞ ⎛ ⎞ ⎛ ⎞ ⎛ ⎞ 0 a a a 1 a a 1 1 a 1 a ⎝0 1 1⎠,⎝0 1 1⎠,⎝a 1 a ⎠,⎝1 1 0⎠, 0 a 1 0 0 a 0 0 0 a 0 0 ⎛

⎞ ⎛ ⎞ ⎛ 1 a 0 a 0 0 0 ⎝1 1 0⎠,⎝1 1 0⎠,⎝a a a 0 a 1 a 1

0 1 1

⎞ ⎛ ⎞ 0 0 0 a a ⎠,⎝0 1 1⎠, a a 1 a

where ‘a’ denotes an arbitrary entry (0 or 1 for binary images). A pixel with entry 1 is set to 0 in the kth step if its corresponding matrix of neighbour entries fullfils all the conditions given by the matrix M(k mod 8)+1 . More precisely in the kth step the matrix M(k mod 8)+1 is applied to the binary image in the following way. Let p k−1 (i, j ) denotes the value of the pixel at position (i, j ) after k − 1 steps of the algorithm. Then p k (i, j ) = 0 if either p k−1 (i, j ) = 0 or if p k−1 (i, j ) fullfils the conditions imposed by M(k mod 8)+1 . Otherwise p k (i, j ) = 1. For example in the 12th step the structural element is given by M(12 mod 8)+1 = M5 . A pixel p 12 (i, j ) is set to 0 if either p 11 (i, j ) = 0 or if p 11 (i, j ) = 1 and the following holds. The pixels p 11 (i − 1, j ) and p 11 (i − 1, j + 1) have to be 1 and the pixels p 11 (i + 1, j − 1), p 11 (i + 1, j ) and p 11 (i + 1, j + 1) must equal 0. The other three pixels in the direct neighbour-

hood, namely p 11 (i − 1, j − 1), p 11 (i, j − 1) and p 11 (i, j + 1) can be of arbitrary value (0 or 1). Notice that if after the k ∗ th step p k ∗ (i, j ) = 0 then p k (i, j ) = 0 for all k ≥ k ∗ . The application of the structuring elements in a rotating fashion to the bi  nary image is performed until i, j p k (i, j ) = i, j p k+8 (i, j ), i.e. the number of pixels that have values equal to 1 has not changed for a whole cycle of the eight structuring elements M(k mod 8)+1 , M((k+1) mod 8)+1 , . . . , M((k+7) mod 8)+1 . Since the number of pixels in the binary image is finite, the algorithm terminates after a finite number of steps and due to the algorithm structure the resulting binary structures have a thickness of one pixel. Algorithm modifications. In order to adjust the skeletonization algorithm described in the previous section to the specific data structure and the aims provided at the beginning of this section, some modifications have to be applied. First of all, the greyscale images have to be transformed into binary images. This is achieved by taking the grey phase of the greyscale image, or in other words the nuclei, as one phase of the binary image, and the stroma and the tumour cell phase without nuclei as the other phase of the binary image. Thereby it is ensured that after skeletonization a cell has exactly one nucleus and vice versa. The skeletonization algorithm itself ensures that cells are constructed according to a nearest neighbour principle with respect to the cell nuclei. What is left to ensure is the property that no stroma tissue belongs to the cell. This is ensured by a slight modification of the skeletonization algorithm. Each pixel belonging to the stroma tissue is identified and it is assured that  C 2014 The Authors C 2014 Royal Microscopical Society, 256, 46–60 Journal of Microscopy 

SPATIAL STATISTICS OF SQUAMOUS CELL CARCINOMAS

the skeletonization procedure does not change pixels next to a pixel that belongs to the tumour-free tissue. Notice that this method makes it necessary to afterwards read the boundary between the stromal tissue and the union of the cell tissue and the cell nuclei to the skeleton structure in order to get bounded cells. Summary of the applied algorithm. The image data given as triphasic greyscale images (Fig. 2A) are transformed into two binary images, one describing the cell nuclei (Fig. 2B) and another one describing the stromal tissue (Fig. 2C). Afterwards a skeletonization algorithm using morphological operators is applied to the binary image describing the cell nuclei, where pixels remain unchanged if they belong or are next to the stromal tissue (as it is marked in the binary image describing the stromal tissue). As a next step the border between the stroma and the cell tissue is added to the skeleton in order to get bounded cells (Fig. 2D). Notice that these cells then fulfill all conditions imposed at the beginning of this section. Finally the skeleton structure is pruned by removing dead ends, i.e. by iteratively removing pixels belonging to the skeleton that only have one neighbouring pixel that also belongs to the skeleton (Fig. 2E). Some examples for results of the algorithm are displayed in Figures 2F and 3B, D, F, where the latter originate from Figures 3A, C and E. Note that it is also possible that no stromal tissue is included in a sample image (see Figs. 3A, B). Statistical analysis of the tessellation data. In order to perform a statistical analysis of the segmented images, various characteristics are estimated. First morphological quantities like the areas and the boundary lengths of the individual segmented cells and the cell nuclei are estimated. This is performed according to the methods described in (Klenk et al., 2006) and (Schmidt & Spodarev, 2005) using the GeoStoch library. GeoStoch is an open-library system based on Java that can be used for stochastic-geometric modelling and spatial statistical analysis of image data (Mayer, 2003; Mayer et al., 2004), see http://www.geostoch.de). Then the following distribution characteristics were extracted: mean, standard deviation, kurtosis and skewness of the areas and of the boundary lengths of the cells and of the nuclei. Point process statistics. The locations of the cell nucleus profiles were analysed by taking the centre of gravity as the midpoint (Figs. 4A, B). The point pattern consisting of these midpoints may be examined by the help of point process characteristics like Ripley’s K -function (Ripley, 1981; Mattfeldt, 2005; Chiu et al., 2013). Intuitively, K (r ) is the mean (= E = expected) number of other points of the process lying within a circle of radius r , centred about a typical point (x, y) of the process, divided by the intensity λ of the process: K (r ) =

E (number of other points with distance ≤ r |point at (x, y)) , (1) λ

 C 2014 The Authors C 2014 Royal Microscopical Society, 256, 46–60 Journal of Microscopy 

49

where the symbol ‘|’ denotes ‘conditional to’. The concept of a typical point of the process means that a point is selected from all points of the process at random with uniform probability. It has to be contrasted against the concept of a uniform random point in the xy-plane where all points (and not just those of the process) are considered. In the evaluation of K (r ), the typical point itself, i.e. the centre of the circle in which the points are counted, is not counted, hence the word ‘other’ in the definition of K (r ) above. As a second point process characteristic also the pair correlation function was considered (see e.g. Mattfeldt et al., 2006; Illian et al., 2008; Chiu et al., 2013). The pair correlation function g(r ) can be considered as the corresponding density function to the K -function. It is defined as: g(r ) =

ρ (2) (r ) , λ2

(2)

where ρ (2) (r ) is the product density of second order. This product density of second order can be explained as follows. If two disks C 1 and C 2 are regarded that have infinitesimal areas F1 and F2 and midpoints x1 and x2 , respectively, the probability for having in each disc at least one point of X is approximately equal to ρ (2) (x1 , x2 )d F1 d F2 . In the case of complete spatial randomness (Poisson process) we have that K P oi (r ) = πr 2 and g P oi (r ) = 1. Edge effects of virtual cells. When looking at the images, it becomes immediately clear that edge-correction techniques have to be used for the construction of the virtual cells, in order to obtain unbiased estimates of characteristics like the mean circumference of a cell. Hence, with respect to morphological characteristics, a minus sampling was performed in the sense that only cells were counted whose midpoints were inside of a reduced sampling window (312 × 312 pixel). Thereby it was assured that only such virtual cells were considered that could be completely observed. Either the virtual cell was lying completely inside the reduced sampling window; or it was transected by the border of this window, and in this case the safety margin of 100 pixels width was sufficient to observe the whole virtual cell. This means that for a pixel belonging to the cell tissue, the minimal distance to the boundary of the tissue is the minimum of the minimal distance to the boundary of the image and the minimal distance to the boundary between tumour cell tissue and stromal tissue. Point processes on partial phases of the tissue. With respect to the point patterns of the centroids of the tumour cell nucleus profiles, one must take into account that the nuclei exist only within the epithelial phase of the tumour, not in the stroma (Figs. 4A, B). In contrast to ordinary stationary point processes, the points lie only in a ‘permitted’ phase whereas there is a ‘forbidden’ phase where the points never occur. Similarly as for stationary point processes that fill the entire R2 , the well-known second-order statistics K (r ) and g(r ) are defined

50

T. MATTFELDT AND F. FLEISCHER

Fig. 2. Image segmentation algorithm. Upper left panel (A): greyscale image. Upper right panel (B): binary image describing cell nuclei. Middle left panel (C): binary image describing epithelial phase and stroma. Middle right panel (D): raw skeleton of cell tissue. Lower left panel (E): skeleton after pruning. Lower right panel (F): greyscale image including skeleton.

for such point processes. For the estimation of both functions, two types of edge corrections are necessary (see below for details). First, one needs minus sampling at the image boundaries to be able to observe the nucleus profiles completely, which is necessary to find the centre of gravity of the nuclei. Second, edge effects here occur not only at the boundaries of the image, but also within the images where the permitted phase is interrupted by the forbidden phase. In order to take also these effects into account, edge-corrected estimators as proposed in (Stoyan & Stoyan, 1994) were used. In order to see the effect of modelling the point process as a point process living in only a restricted reference phase or as a stationary process, the two K -functions and the two pair correlation functions g(r ) were computed for all visual fields of all cases, and averaged (see Figs. 5A, B). Clearly modelling the point patterns related to the restricted reference space of epithelial tumour tissues leads to less deviation of K (r ) and g(r ) at long distances r from the reference functions for a homogeneous planar Poisson point prcocess. Also one sees that the difference between the two curves is considerable and non-negligible.

It was explored whether our patterns could be considered as realizations of ‘interrupted point processes’ (Stoyan, 1979; Kautz et al., 2011). Such processes result from a dependent thinning of a stationary point process  with intensity λ in R2 : A random set A is generated in R2 ; within this set we have the random point process of intensity λ, and outside of A we have no points (dependent thinning). The latter phase is called the interrupting field. In the case of interrupted point processes, it is required that the point process  is independent from the random set A (Stoyan, 1979; Kautz et al., 2011). It was explored whether the intensity of the point process of the tumour cell nuclei per unit area of the epithelial phase and the area fraction of the epithelial phase per unit tissue area were correlated. A highly significant negative correlation between these two parameters was found (r = −0.48, p < 0.001) (Fig. 6). This means that fields of vision containing high fractions of epithelial tumour tissue (and hence relatively few stroma), are on average containing less nuclei per unit area of epithelium than those with low fractions of epithelial tumour tissue and high stromal content. This is biologically plausible because areas deep inside the tumour tissue often tend to higher

 C 2014 The Authors C 2014 Royal Microscopical Society, 256, 46–60 Journal of Microscopy 

SPATIAL STATISTICS OF SQUAMOUS CELL CARCINOMAS

51

Fig. 3. Upper panels: three original samples (A, C, E). Lower panels: same samples with segmented cell borders (B, D, F) .

Fig. 4. (A) Triphasic image with nuclei (grey), epithelial tumour issue without nuclei (white), and stroma (black). (B) Same image, where the centroids of the tumor cell nucleus profiles have been indicated as crosses. The crosses represent a point pattern inhabitating the white phase only. The stroma contains no points. (C) Same image, but showing the safety margin of 50 pixels on each side yielding a reduced observation window of 412 × 412 pixels (minus sampling).

differentiation due to keratinization, which reduces the density of the tumour cells, whereas margins of tumour tissue often appear more cellular. In any case, this result showed that the concept of an interrupted point process (in the strict

 C 2014 The Authors C 2014 Royal Microscopical Society, 256, 46–60 Journal of Microscopy 

sense) could not be used for modelling our data set. Instead, one might consider a primary point process  in R2 , which is then thinned in a dependent manner by taking out all areas which belong to the forbidden phase; but in contrast to the interrupted

52

T. MATTFELDT AND F. FLEISCHER

Fig. 5. (A) Plots of the mean functions K (r ) of the centroids of the nucleus section profiles of all patterns. Blue: modelling as a stationary point process; red: modelling as a point process with restricted reference space; black: reference function K P oi (r ) = πr 2 . (B) Plots of the mean functions g(r ) of the centroids of the nucleus section profiles of all patterns. Blue: modelling as stationary point process; red: modelling as as a point process with restricted reference space; black: reference function and g P oi (r ) = 1.

point process in the strict sense, it is not a priori required that  and the random set A are independent. Such a point process might be called a point process on a restricted part of the reference space, or a restricted point process, for the sake of simplicity.

In practice, the epithelial phase was segmented as white, and the stromal phase as black. Then the centroids of the nucleus profiles were detected. Finally, the resulting point pattern was superponed on the binary image (Fig. 4B). After this, only the white phase is inhabitated by the point process, whereas  C 2014 The Authors C 2014 Royal Microscopical Society, 256, 46–60 Journal of Microscopy 

SPATIAL STATISTICS OF SQUAMOUS CELL CARCINOMAS

53

Fig. 6. Horizontal axis: area fraction of epithelium, vertical axis: intensity of point process of nucleus profiles per unit area of the epithelial phase. The points reflect the 53 cases. Note strongly negative correlation.

the black phase is the forbidden phase. The estimation of the intensity, of K (r ) and of g(r ) was performed with the spatstat package (Baddeley & Turner, 2005), which allows for irregularly delimited fields of vision. In the case of the nuclei which are much smaller than the virtual cells, a minus sampling with a slightly reduced sampling window was considered as sufficient (412 × 412 pixel) (Fig. 4C). With respect to point process characteristics edge-corrected estimators as proposed in (Stoyan & Stoyan, 1994) were used, where special care was spent to the fact that the stromal tissue must also be considered as a boundary for this purpose. The computation of g(r ) was performed similarly as described in (Mattfeldt et al., 2006). Briefly, an estimator for g(r ) is given by  g(r ) =

 (2) (r )  , λ2

(3)

where 1  (2) (r ) =  2πr

X i ,X j ∈Wi = j

kh (r − ||X i − X j ||) |WXi ∩ WX j |

(4)

is an estimator for (2) (r ), the product density of second order.  (2) (r ), the denominator |W ∩ W | For the estimator  Xi Xj ensures edge-correction. The term kh (x) denotes a kernel function, which is used for smoothing. We used the Epanech C 2014 The Authors C 2014 Royal Microscopical Society, 256, 46–60 Journal of Microscopy 

nikov kernel

 x2 (5) 1 − 2 1(−h,h) (x), h √ with a bandwidth h = 0.15/  λ, see (Mattfeldt et al., 2006) for further details. The symbols  λ and λ2 denote denote the estimates of the intensity λ and of the squared intensity λ2 , respectively. kh (x) =

3 4h

Logistic regression analysis Logistic regression is a mathematical modelling approach that can be used to describe the relationship of several independent variables to a binary variable l(z). Such a variable assumes only two values, here it is the nodal state, i.e. pN0 versus positive nodal states pN+ (pN1 or pN2), which were coded as 0 and 1, respectively. The notions pN1 and pN2 relate to the number, size and location of the lymph node metastases. The groups pN1 and pN2 were united to a pN+ group which represents the existence of a lymphatic spread as a whole as an indicator of aggressiveness; moreover, a finer subdivision would have led to very small groups. We consider now the linear function z = f (X) of a k-dimensional input vector X = (x1 , · · · , xk ), where the xi are the k independent variables: z = α + β1 x1 +  β2 x2 + · · · + βk xk = α + ik=1 βi xi . The function z = f (X) can be interpreted as an index of combined risk factors. The logistic model can be written as k

l(z) = 1/(1 + e −z ) = 1/(1 + e −(α+

i =1

βi xi )

),

(6)

54

T. MATTFELDT AND F. FLEISCHER

where the parameters α and βi are obtained (fitted) according to the principle of maximum likelihood from the data (Kleinbaum, 1994). Essentially, this leads to a mapping of a linear combination of the independent variables onto {0, 1}. In the present study we used the procedure of logistic regression as implemented in the statistic software package of SAS/STAT (proc logistic) with the stepwise selection option (SAS Institute, 2000). In this case the independent variables are entered into and removed from the model in such a way that each forward selection step is followed by one or more backward elimination steps. The stepwise selection process terminates if no further variable can be added to the model. Simplifying, this means that significant prognostic (predictive) factors are identified: all other variables are eliminated and only significant predictors remain in the model. Results Basic data The data set included 53 squamous cell carcinomas of the head and neck: 10 cases from the oral cavity, 25 cases from the oropharynx and 18 cases from the hypopharynx. The nodal negative cases included the following tumour stages: 5 pT1, 10 pT2, 1 pT3 and 3 pT4. In the nodal positive cases, the following tumour stages were found: 3 pT1, 10 pT2, 9 pT3 and 12 pT4. On the whole, the median tumour stage was 2 in the nodal negative and 3 in the nodal positive group. The difference between the tumour stages was significant according to a Kruskal-Wallis test ( p < 0.01). Grading ranged from 1 to 3 with a median grade of 2 in both groups (difference according to a Kruskal–Wallis test not significant). Group comparisons: results of spatial statistics Group comparison of point process statistics. Exploratory point process statistics were applied separately to each of the two images per case, which resulted in two estimates of the intensity, of K (r ) and of g(r ) per case. For each case, the mean of these pairs was computed. The intensity of the point process of the nucleus profile centroids with the epithelial tumour component as reference space did not differ significantly between the groups nodal-negative and the nodal-positive group (0.001122 vs. 0.001115, N.S.). The mean area fraction occupied by epithelial tumour tissue also was very similar (nodalnegative and the nodal-positive group (0.78 vs. 0.76, N.S.). Figures 7A– C shows the mean K -functions of the centres of the nucleus profiles in the two groups. The estimated functions are nearly identical. The curves begin with a segment where K (r ) = 0 due to the positive size of the nuclei which cannot overlap. Thereafter the function approaches the theoretical K -function of a homogeneous 2D Poisson point process. Figures 8A– C shows the mean g-functions of the tu-

mour cell nucleus profiles of our two groups of squamous cell carcinomas. Again the hard-core property is evident from the plots, and for high r -values spatial correlation vanishes, hence the mean functions approach the value 1 in both groups. Group comparisons of parameters of tessellation cells. For both groups, the frequency distributions of the areas and of the boundary lengths of the tessellation cells were estimated for the set of all virtual cells on the two images per case. From these data, the mean value, the S D , the skewness and the kurtosis were computed. Then the group mean values of these characteristics were tested for significant differences by means of the nonparametric Kruskal-Wallis test (Sachs, 2004). The results are summarized in Table 1. For the mean, SD, skewness and kurtosis of the area distributions of the cells, no significant difference between the nodal negative and the nodal positive groups could be found. There was a significant increase of the mean boundary length of the cells in the nodal positive group ( p < 0.05). The other distributional characteristics of the cell boundary length did not differ significantly between the groups. Group comparisons of parameters of tumour cell nucleus profiles. In the same manner as described in the previous section, estimates of the mean value, the S D , the skewness and the kurtosis of the area and the boundary length of the tumour cell nucleus section profiles were computed. The Kruskal-Wallis test was used to test for significant differences between group mean values. There were significant increases of the mean nuclear area ( p < 0.01) and of the standard deviation of the nuclear area ( p < 0.05) in the nodal positive group (see Table 2). Skewness and kurtosis of the nuclear area, and all distributional characteristics of the nuclear boundary length showed no significant differences between the groups. Results of logistic regression analysis Applying the aforementioned stepwise algorithm for logistic regression to our data, we used the lymph node state pN0 and pN+ = pN1 or pN2 as dependent variable l(z) ∈ {0, 1}. When tumour grade, pT-stage, and mean, standard deviation, kurtosis and skewness of the areas and of the boundary lengths of the cells and of the nuclei were used as predictor variables, the pT-stage and the variable ’mean area of the tumour cell nucleus profiles’ met the 5% significance level for entry into the logistic regression model (see eq. 6). In general, the coefficients of the logistic regression are changed in such a way that the fit (the likelihood) is maximized. The likelihood ratio test computes the ratio of the log likelihood of a certain model as compared to the log likelihood corresponding to a global null hypothesis, where all values of the βi coef 2 has a χ 2 (chi square) ficients are zero. The test statistic χ

 C 2014 The Authors C 2014 Royal Microscopical Society, 256, 46–60 Journal of Microscopy 

SPATIAL STATISTICS OF SQUAMOUS CELL CARCINOMAS

55

Fig. 7. (A) Estimated mean function K (r ) (green) for the nodal-negative group with 95% confidence interval (borders: blue, red); black: reference function K P oi (r ) = πr 2 corresponding to a Poisson point process. (B) Estimated mean function K (r ) for the nodal-positive group (green) with 95% confidence interval (blue, red); black: reference function. (C) Estimated mean functions K (r ) for nodal negative (red) and nodal positive (blue) groups. Black: reference function.  C 2014 The Authors C 2014 Royal Microscopical Society, 256, 46–60 Journal of Microscopy 

56

T. MATTFELDT AND F. FLEISCHER

Fig. 8. (A) Estimated mean function g(r ) for the nodal-negative group (green) with 95% confidence interval (borders: blue, red); black: reference function g P oi (r ) = 1 corresponding to a Poisson point process. (B) Estimated mean function g(r ) for the nodal-positive group with confidence interval and reference function; colours as in upper panel. (C) Estimated mean functions g(r ) for nodal negative (red) and nodal positive (blue) groups. Black: reference function.  C 2014 The Authors C 2014 Royal Microscopical Society, 256, 46–60 Journal of Microscopy 

SPATIAL STATISTICS OF SQUAMOUS CELL CARCINOMAS

57

Table 1. Group comparisons distribution characteristics of virtual cell size. Group Number of cases Variable

Nodal neg. cases

Nodal pos. cases

Mean

19 SD

Mean

34 SD

p-value of Kruskal-Wallis-test

Mean cell area SD of cell area Skewness of cell area Kurtosis of cell area

786 397 1.43 3.03

228 141 0.38 1.97

846 404 1.35 3.51

240 155 0.87 7.82

0.481 0.941 0.099 0.087

Mean cell boundary length SD of cell boundary length Skewness of cell boundary length Kurtosis of cell boundary length

89 29 0.74 0.77

32 5 0.28 1.09

108 29 0.75 0.98

26 6 0.40 1.83

0.027∗ 0.738 0.941 0.867

In both tables, the abbreviation SD indicates standard deviation. The data represent the mean characteristics of the distributions of the cell and nucleus profile areas and boundary lengths, i.e. for each variable the mean values and standard deviations between cases within the two groups are indicated. The Kruskal–Wallis test yields a statistics χ 2 with 1 degree of freedom. The p-values corresponding to this χ 2 -value give the probability of error when the null hypothesis, that the two values stem from the same population, is rejected. The p-values are marked by asterisks if significant: *p < 0.05, **p < 0.01. Table 2. Group comparisons distribution characteristics of nucleus profile size. Group Number of cases Variable Mean nucleus profile area SD of nucleus profile area Skewness of nucleus profile area Kurtosis of nucleus profile area Mean nucleus profile boundary length SD of nucleus profile boundary length Skewness of nucleus profile boundary length Kurtosis of nucleus profile boundary length

Nodal neg. cases

Nodal pos. cases

Mean

19 SD

Mean

34 SD

p-value of Kruskal-Wallis-test

172 106 1.82 3.77

76 81 0.69 5.09

239 141 1.65 4.39

91 55 0.77 5.03

0.009∗∗ 0.024∗ 0.331 0.894

59.39 29.90 1.19 2.65

8.13 17.07 0.64 4.98

64.32 25.48 1.39 3.03

10.06 11.34 0.54 2.93

0.063 0.838 0.988 0.092

In both tables, the abbreviation SD indicates standard deviation. The data represent the mean characteristics of the distributions of the cell and nucleus profile areas and boundary lengths, i.e. for each variable the mean values and standard deviations between cases within the two groups are indicated. The Kruskal–Wallis test yields a statistics χ 2 with 1 degree of freedom. The p-values corresponding to this χ 2 -value give the probability of error when the null hypothesis, that the two values stem from the same population, is rejected. The p-values are marked by asterisks if significant: *p < 0.05, **p < 0.01.

distribution with as many degrees of freedom as the model has variables. For our above-mentioned variables, the corresponding value of the likelihood ratio test with 2 degrees of 2 = 8.58, which means significance at freedom amounted to χ ( p < 0.01). When only tumour grade, mean, standard deviation, kurtosis and skewness of the areas and of the boundary lengths of the cells and of the nuclei were used as predictor variables, only the variable ’mean area of the tumour cell nucleus profiles’ met the 5% significance level for entry into the model. The corresponding value of the likelihood-test 2 = 7.17 ( p < 0.01). with 1 degree of freedom amounted to χ The other variables were eliminated from the model as noncontributory.

 C 2014 The Authors C 2014 Royal Microscopical Society, 256, 46–60 Journal of Microscopy 

Discussion Point process statistics for quantitative evaluation of tumour tissue In the present study, the spatial distribution of the centroids of tumour cell nucleus profiles was described in terms of Ripley’s K -function and the pair correlation function. This approach gives additional information to the intensity of the process, i.e. the number of nucleus profiles per unit area. It gives the opportunity to distinguish between point processes that have the same intensity, but differ by spatial arrangement. In a recent study, similar methods were applied for the study of the spatial distribution of capillary profiles in prostatic

58

T. MATTFELDT AND F. FLEISCHER

tissue (Mattfeldt et al., 2006; Mattfeldt et al., 2007). There was no significant difference between the nodal-negative and the nodal-positive group w.r.t. intensity and both second-order statistics. Apparently for the first time, a point process with a reference space restricted to a certain phase of the tissue has been used for the description of a histological point pattern. It should be noted that such modelling is more natural for many such applications than the usual approach, where it is assumed that the whole area of the histological section is the reference space. For example, blood capillaries only lie within the stromal phase of a tumour; in this case, the epithelial phase is the forbiddden phase. There were systematic differences in the resulting estimates of K (r ) and g(r ) depending on consideration of the point process as stationary or restricted, see Figures 5A, B. The statistics K (r ) and g(r ) lie distinctly higher than their reference function K P oi (r ) = πr 2 and g P oi (r ) = 1 at large r values, if the point patterns are considered as realizations of a stationary point process occupying an uninterrupted reference space. This indicates a clustering of the points at large r . This effect is not found when modelling the point patterns as a process on the restricted reference space. This clustering interaction may stem from the fact that the points in each separate epithelial domain appear as large clusters if the whole tissue is considered as the reference space, see Figure 4B. One cannot say that only one of the two approaches is correct and the other is wrong; the resulting statistics only reflect two different views. As the two components of the process, i.e. the area process describing the reference phase of the points, and the intensity of the points do not appear to be independent, our particular data are not adaequately described as an interrupted point process. However, this finding precludes in no way that the model of an interrupted point process might be useful in other histopathological applications. Tessellation methods for quantitative texture analysis of tumour tissue In the present study the tumour texture was not only quantified in terms of an explorative study of the point process of the centroids of tumour cell nucleus profiles, but in addition a virtual tessellation was generated which led to a partition of the epithelial image phase into planar cells, each related to a nucleus profile. In principle this yields additional information to the aforementioned second-order statistics K (r ) and g(r ), because it provides data on the size distribution of the planar tessellation cells. This distribution itself can be estimated and plotted, or one may extract simple distribution characteristics such as mean, S D , skewness, kurtosis, etc. This approach was applied for the first time to squamous cell carcinomas. It turned out that the mean boundary length of the cells in the nodal-positive group was significantly larger than that of the nodal-negative group. Note that the application of our tessellation technique is by no means restricted to squamous cell

carcinomas, but can be applied also to other frequent solid tumours such as squamous cell carcinomas in other organs (e.g. the skin, the lung, or the uterine cervix), or in squamous cell carcinomas and urothelial carcinomas of the urinary bladder. In such applications, it may well come out that the characteristics estimated by the tessellation method bring additional information for the distinction of closely similar looking tumour tissues. Tessellation methods have been applied for the study of tumour tissue only very scarcely (Sudbø et al. 2000a, 2000b, 2000c); for applications in electron microscopy see (Beil et al., 2005a, 2005b, 2006). It must be admitted that the amount of overhead for the present methodology was high, but it seems not utopic that in the future the whole process may become automatic. It is only necessary to segment the tumour cells and the nuclei by appropriate stains, for this aim automatic immunohistochemical stains could be of value. Once the segmentation is performed, the whole algorithm runs automatically. Results on tumour cell nucleus profiles A main result of our study was the finding that the two groups showed significant differences w.r.t. quantitative characteristics of the tumour cell nucleus profiles. The mean area and the SD of the areas was increased in the group of the nodal positive cases. Further distributional characteristics such as the skewness and the kurtosis showed no significant difference. Interestingly, only the mean nuclear area and tumour stage proved to be significant in a multivariate approach, when multiple input variables were tested for significance in a stepwise logistic regression analysis on nodal state. All other characteristics (mean, standard deviation, skewness, and kurtosis of the area and boundary length of the cells; mean, standard deviation, skewness and kurtosis of nuclear profile boundary length; standard deviation, skewness and kurtosis of nuclear profile area) were eliminated as insignificant in the logistic model. Moreover, the mean nuclear area was a better predictor of nodal state than conventional grading. The finding that the mean nucleus profile area was higher in the nodal positive group fits well to the observations of other authors. In general, nuclear morphometry has suggested that more aggressive squamous cell carcinomas have larger nuclei. However, it has been debated whether the nuclear morphometry may be considered as a significant predictor of nodal state. Some have provided data corroborating this view (Sekine et al., 2003; Lavie et al., 2006; Natarajan et al., 2010), whereas other authors could not confirm the prognostic value of 2D nuclear morphometry (Dobros et al., 1999). Also it has been questioned whether the stereological parameter of nuclear size, the volume weighted mean nuclear volume, is of prognostic value w.r.t. nodal state (Bjerregaard et al., 1993). The conflicting results may in part be due to the fact that squamous cell carcinomas of such different locations such as the oro- and hypopharynx, the oral mucosa, the skin, the esophagus, the  C 2014 The Authors C 2014 Royal Microscopical Society, 256, 46–60 Journal of Microscopy 

SPATIAL STATISTICS OF SQUAMOUS CELL CARCINOMAS

vulva and the uterine cervix were considered, which may well have different biological properties. The present paper appears to confirm the affirmative position w.r.t. the prognostic value of nuclear morphometry and also indicates that the virtual cells of pN+ cases differ from pN0 cases. Acknowledgements The skillful technical assistance of Rolf Kunft and Michael Held is gratefully acknowledged.

References Barnes, L., Eveson, J.W., Reichart, P. & Sidransky, D. (Eds.) (2005) World Health Organization classification of tumours. Pathology and Genetics of Head and Neck Tumours. IARC Press, Lyon. Baddeley, A. & Turner, R. (2005) Spatstat: an R package for analyzing spatial point patterns. J. Stat. Software 12, 1–42. Beil, M., Braxmeier, H., Fleischer, F., Schmidt, V. & Walther, P. (2005a) Quantitative analysis of keratin filament networks in scanning electron microscopy images of cancer cells. J. Microsc. 220, 84–95. Beil, M., Eckel, S., Fleischer, F., Schmidt, H., Schmidt, V. & Walther, P. (2006) Fitting of random tessellation models to cytoskeleton network data. J. Theoret. Biol. 241, 62–72. Beil, M., Fleischer, F., Paschke, S. & Schmidt, V. (2005b) Statistical analysis of 3D centromeric heterochromatin structure in interphase nuclei. J. Microsc. 217, 60–68. Bjerregaard, B., Andreasson, B., Visfeldt, J. & Bock, J.E. (1993) The significance of histology and morphometry in predicting lymph node metastases in patients with squamous cell carcinoma of the vulva. Gynecol. Oncol. 50, 323–329. Bundgaard, T., Sorensen, F.B., Gaihede, M., Sogaard, H. & Overgaard, J. (1992) Stereologic, histopathologic, flow cytometric, and clinical parameters in the prognostic evaluation of 74 patients with intraoral squamous cell carcinomas. Cancer 70, 1–13. Chiu, S.N., Stoyan, D., Kendall, W.S. & Mecke, J. (2013) Stochastic Geometry and its Applications. 3rd edn. Wiley, Chichester. Dobros, W., Gil, K., Chlap, Z. & Olszewski, E. (1999) The use of nuclear morphometry for the prediction of survival in patients with advanced cancer of the larynx. Eur. Arch. Oto-Rhino-Laryngol. 256, 257–261. Dreyer, T., Knoblauch, I., Doudkine, A., MacAulay, C.E., Garner, D., Palcic, B. & Popella, C. (2001) Nuclear texture features for classifying benign vs. dysplastic or malignant squamous epithelium of the larynx. Anal. Quant. Cytol. Histol. 23, 193–200. Fleischer, F., Gloaguen, C., Schmidt, H., Schmidt, V. & Schweiggert, F. (2008) Simulation algorithm of typical modulated Poisson-Voronoi cells and application to telecommunication network modelling. Japan. J. Indust. Appl. Math. 25, 305–330. Glazer, E.S., Bartels, P.H., Prasad, A.R., Yozwiak, M.L., Bartels, H.G., Einspahr, J.G., Alberts, D.S. & Krouse, R.S. (2011) Nuclear morphometry identifies a distinct aggressive cellular phenotype in cutaneous squamous cell carcinoma. Cancer Prev. Res. 4, 1770–1777. Gloaguen, C., Fleischer, F., Schmidt, H. & Schmidt, V. (2006) Fitting of stochastic telecommunications network models via distance measures and Monte-Carlo tests. Telecommun. Syst. 31, 353– 377.

 C 2014 The Authors C 2014 Royal Microscopical Society, 256, 46–60 Journal of Microscopy 

59

Huang, T., Chen, M., Wu, M. & Wu, X. (2008) Image analysis of DNA content and nuclear morphometry for predicting radiosensitivity of nasopharyngeal carcinoma. Anal. Quant. Cytol. Histol. 30, 169–174. Illian, J., Penttinen, A., Stoyan, H. & Stoyan, D. (2008) Statistical Analysis and Modelling of Spatial Point Patterns. Chichester, Wiley. J¨ahne, B. (2001) Digital Image Processing. Springer, Berlin. Janot, F., Klijanienko, J., Russo, A., et al. (1996) Prognostic value of clinicopathological parameters in head and neck squamous cell carcinoma: a prospective analysis. Br. J. Cancer 73, 531–538. Kautz, M., Berger, U., Stoyan, D., et al. (2011) Desynchronizing effects of lightning strike disturbances on cyclic forest dynamics in mangrove plantations. Aquat. Bot. 95, 173–181. Kleinbaum, D.G. (1994) Logistic Regression: A Self-learning Text. Springer, New York. Klenk, S., Schmidt, V. & Spodarev, E. (2006) A new algorithmic approach to the computation of Minkowski functionals of polyconvex sets. Comp. Geom. Th. Appl. 34, 127–148. Lavie, O., Maini, I., Pilip, A., et al. (2006) Computerized nuclear morphometry for the prediction of inguinal lymph nodes metastases in squamous cell carcinoma of the vulva. Int. J. Gyn. Cancer 16, 556– 561. Liu, O., Zhang, H., Wang, Y., et al. (2012) Stereology study of oral verrucous carcinoma. J. BUON 17, 343–349. Maiolino, P., Restucci, B., Papparella, S. & De Vico, G. (2002) Nuclear morphometry in squamous cell carcinomas of canine skin. J. Comp. Pathol. 127, 114–117. Martano, M., Damiano, S., Restucci, B., Paciello, O., Russo, V. & Maiolino, P. (2006) Nuclear morphometry in canine acanthomatous ameloblastomas and squamous cell carcinomas. Eur. J. Histochem. 50, 125– 130. Mattfeldt, T. (2005) Explorative statistical analysis of planar point processes in microscopy. J. Microsc. 220, 131–139. Mattfeldt, T., Eckel, S., Fleischer, F. & Schmidt, V. (2006) Statistical analysis of reduced pair correlation functions of capillaries in the prostate gland. J. Microsc. 223, 107–119. Mattfeldt, T., Eckel, S., Fleischer, F. & Schmidt, V. (2007) Statistical modelling of the geometry of prostatic capillaries on the basis of stationary Strauss hard-core processes. J. Microsc. 228, 272–281. Mayer, J. (2003) On quality improvement of scientific software: theory, methods, and application in the GeoStoch Development. Doctoral Dissertation, University of Ulm. Mayer, J., Schmidt, V. & Schweiggert, F. (2004) A unified simulation framework for spatial stochastic models. Simul. Modelling Pract. Theory 12, 307–326. Narvaez, D., Kanitakis, J., Euvrard, S., Schmitt, D., Faure, M. & Claudy, A. (1997) Comparative nuclear morphometric analysis of aggressive and non-aggressive squamous cell carcinomas of the skin. Acta Derm.Venereol. 77, 115–117. Natarajan, S., Mahajan, S., Boaz, K. & George, T. (2010) Prediction of lymph node metastases by preoperative nuclear morphometry in oral squamous cell carcinoma: a comparative image analysis study. Ind. J. Cancer 47, 406–411. Okon, K., Basta, A. & Stachura, J. (1998) Morphometry separates squamous cell carcinoma of the vulva into two distinct groups. Polish J. Pathol. 49, 293–295. Pektas, Z.O., Keskin, A., Gunhan, O. & Karslioglu, Y. (2006) Evaluation of nuclear morphometry and DNA ploidy status for detection of malignant and premalignant oral lesions: quantitative cytologic assessment

60

T. MATTFELDT AND F. FLEISCHER

and review of methods for cytomorphometric measurements. J. Oral Maxillofac. Surg. 64, 628–635. Ripley, B.D. (1981) Spatial Statistics. Wiley, Chichester. Sachs, L. (2004) Angewandte Statistik. 11th edn. Springer, Berlin. SAS Institute (2000) SAS/STAT User’s Guide, Version 8. Cary, NC. Schmidt, V. & Spodarev, E. (2005) Joint estimators for the specific intrinsic volumes of stationary random sets. Stoch. Proc. Appl. 115, 959– 981. Sekine, J., Uehara, M., Hideshima, K., Irie, A. & Inokuchi, T. (2003) Predictability of lymph node metastases by preoperative nuclear morphometry in squamous cell carcinoma of the tongue. Cancer Det. Prev. 27, 427–433. Serra, J. (1988) Image Analysis and Mathematical Morphology. Academic Press, London. Soille, P. (2003) Morphological Image Analysis. Springer, Berlin. Sorensen, F.B., Bichel, P. & Jakobsen, A. (1992) DNA level and stereologic estimates of nuclear volume in squamous cell carcinomas of the uterine cervix. A comparative study with analysis of prognostic impact. Cancer 69, 187–199.

Stoyan, D. (1979) Interrupted point processes. Biom. J. 7, 607–610. Stoyan, D. & Stoyan, H. (1994) Fractals, Random Shapes and Point Fields. Methods of Geometrical Statistics. Wiley, Chichester. Sudbø, J., Bankfalvi, A., Bryne, M., et al. (2000b) Prognostic value of graph theory-based tissue architecture analysis in carcinomas of the tongue. Lab. Invest. 80, 1881–1889. Sudbø, J., Marcelpoil, R. & Reith, A. (2000) Caveats: numerical requirements in graph theory based on quantitation of tissue architecture. Anal. Cell. Pathol. 21, 59–69. Sudbø, J., Marcelpoil, R. & Reith, A. (2000a) New algorithms based on the Voronoi diagram applied in a pilot study on normal mucosa and carcinomas. Anal. Cell. Pathol. 21, 71–86. Wahi, P.N., Cohen, B., Luthra, U.K. & Torloni, H. (1971) Histological typing of oral and oropharyngeal tumours. World Health Organization, Geneva. Weyn, B., Tjalma, W., Van de Wouwer, G., Van Daele, A., Scheunders, P., Jacob, W. & Van Marck, E. (2000) Validation of nuclear texture, density, morphometry and tissue syntactic structure analysis as prognosticators of cervical carcinoma. Anal. Quant. Cytol. Histol. 22, 373–382.

 C 2014 The Authors C 2014 Royal Microscopical Society, 256, 46–60 Journal of Microscopy