Breast Carcinoma, Intratumour Heterogeneity and Histological Grading, Using Geostatistics

Tumour progression is currently believed to result from genetic instability. Chromosomal patterns specific of a type of cancer are frequent even though phenotypic spatial heterogeneity is omnipresent. The latter is the usual cause of histological grading imprecision, a well documented problem, without any fully satisfactory solution up to now. The present article addresses this problem in breast carcinoma. The assessment of a genetic marker for human tumours requires quantifiable measures of intratumoral heterogeneity. If any invariance paradigm representing a stochastic or geostatistic function could be discovered, this might help in solving the grading problem. A novel methodological approach using geostatistics to measure heterogeneity is used. Twenty tumours from the three usual (Scarff‐Bloom and Richardson) grades were obtained and paraffin sections stained by MIB‐1 (Ki‐67) and peroxidase staining. Whole two‐dimensional sections were sampled. Morphometric grids of variable sizes allowed a simple and fast recording of positions of epithelial nuclei, marked or not by MIB‐1. The geostatistical method is based here upon the asymptotic behaviour of dispersion variance. Measure of asymptotic exponent of dispersion variance shows an increase from grade 1 to grade 3. Preliminary results are encouraging: grades 1 and 3 on one hand and 2 and 3 on the other hand are totally separated. The final proof of an improved grading using this measure will of course require a confrontation with the results of survival studies.

Tumour progression is currently believed to result from genetic instability. Chromosomal patterns specific of a type of cancer are frequent even though phenotypic spatial heterogeneity is omnipresent. The latter is the usual cause of histological grading imprecision, a well documented problem, without any fully satisfactory solution up to now. The present article addresses this problem in breast carcinoma. The assessment of a genetic marker for human tumours requires quantifiable measures of intratumoral heterogeneity. If any invariance paradigm representing a stochastic or geostatistic function could be discovered, this might help in solving the grading problem. A novel methodological approach using geostatistics to measure heterogeneity is used. Twenty tumours from the three usual (Scarff-Bloom and Richardson) grades were obtained and paraffin sections stained by MIB-1 (Ki-67) and peroxidase staining. Whole two-dimensional sections were sampled. Morphometric grids of variable sizes allowed a simple and fast recording of positions of epithelial nuclei, marked or not by MIB-1. The geostatistical method is based here upon the asymptotic behaviour of dispersion variance. Measure of asymptotic exponent of dispersion variance shows an increase from grade 1 to grade 3. Preliminary results are encouraging: grades 1 and 3 on one hand and 2 and 3 on the other hand are totally separated. The final proof of an improved grading using this measure will of course require a confrontation with the results of survival studies.

Introduction
Determining the accuracy of the histological prognostic grading of a cancer is a crucial matter in terms of therapeutic decisions. In breast carcinoma, various morphological, cytonuclear and differentiation characteristics and number of mitoses, are used to construct histological grading systems. The best known is the Scarff-Bloom and Richardson (SBR) [6]. It has been shown that this grading system (as many others) is far from being perfect [2,14]. One of the concerns with the histopathological grading of breast cancer is that tumour grading is a subjective evaluation that may affect reproducibility [10,34].
Even if one could capture the whole histological pattern of a tumour (an impossible endeavour in routine pathological practice), the problem would remain of selecting the worst sections (but they might be all more or less bad depending upon the exact place) and then finding a model to characterise the intensity of the phenomenon and its variance. We had previously attempted this by determining a quantitative morphometric index from a few tumours and showing the feasibility of such an approach. Although the work was too tedious to study important series of patients, we showed that heterogeneity in breast cancer was everywhere, at all levels (blocks, sections and slides) [11]. A previous study had proposed a measure of heterogeneity based upon the Minimal Spanning Tree [13].
A study of tumoral heterogeneity implies a systematic sampling. Recently, several experiments were carried out using several samples from the same tumour, in breast [11,16,27] and prostate [26]. Heterogeneity was found at all levels of sampling. Therefore there were hardly any feasible grading recommendations.
A recent meeting of the Kananaskis Tumour Heterogeneity working group resulted in the definition of a general framework to assess the use of genetic markers of tumour progression in the context of intratumour heterogeneity [8]. This meeting emphasized the need for a reproducible and quantifiable measure of heterogeneity.
If an invariance paradigm representing a stochastic or geostatistic function could be established to measure heterogeneity, a novel grading method could be proposed to assess heterogeneity using geostatistics.
In this article, a new measure of heterogeneity is introduced. It is based on the asymptotic behaviour of dispersion variance in a spatial point process defined by the existence of a positive cell for a given marker at each point of the domain. The relevancy of this measure to histopathological grading is assessed by a nonparametric discriminant analysis.

Samples
Twenty carcinomas of the breast were surgically removed.
Blocks were routinely fixed in 10% neutral buffered formalin for 24 hr at room temperature and then paraffin embedded according to standard pathology labora-  II  IDC  10  1  II  IDC   11  1  II  IDC  12  1  II  ILC   13  1  II  IDC   14  1  II  IDC  15  1  II  IDC   16  1  III  IDC   17  1  III  IDC  18  2  III  IDC   19  1  III  IDC  20 1 III IDC * IDC: infiltrating ductal carcinoma, ILC: infiltrating lobular carcinoma. tory practice. This will allow archival material to be used in future retrospective studies. The largest diametral sections only were studied.
All tumours were classified according to the grading method of Scarff-Bloom and Richardson [6]. Information about these specimens is to be found in Table 1.

Microwave treatment
After removing paraffin with toluen, sections were gradually rehydrated in ethanol and distilled water. They were heated in a microwave oven (450 W) for 15 min, in a 10 mM sodium citrate buffer, pH 7.3. Sections were then cooled using water.

Labelling
Paraffin sections (3-5 µm thick) were cut and labelled by monoclonal anti-MIB-1 (equivalent of Ki-67; Immunotech, France) because of its cell proliferation properties [7], using a dilution of 1/20. This marker was detected by a mouse antibody labelled with perox-idase and revealed by DAB (brown) with amplification (Ventana Medical Systems, USA).
Controls were performed to check the specificity of the MIB-1 marker, by omitting the first incubation. All controls were negative and are, therefore, not mentioned further. Nuclei were stained with Harris-hematoxyline solution.

Image acquisition
This was achieved by using a computer in semiautomatic mode. In each field, the coordinates of MIB-1 labelled epithelial nuclei were obtained (software "IPS Image", Alcatel TITN, France) and stored. Full automation was difficult because of segmentation problems between epithelial and stromal cells and frequently overlapping epithelial cells. Labelled cell marks were identified visually, avoiding any artefact. Clumped cells were also separated by visual inspection.
Stromal cells could be rejected visually as they were clearly more elongated.

Morphometric model
The position of markers (epithelial cell nuclei and MIB-1) allowed their sorting into a morphometric grid, obtaining the dispersion variance data.
At each step, a tissue portion of dimension 200 µm × 200 µm, at microscope magnification 25 × 1.6, was numerized into 550×550 pixels (Fig. 1). This was done in semi-automatic mode (picking cell coordinates one by one and saving them in the computer). At the next step, the whole image window was moved by 200 µm one way until the edges of the tumour were reached, after which the scanning passed above to the next field, and so on until the whole section had been explored using staggered row scanning. All fields must be contiguous. The data series obtained for each tumour block consisted of the 2D coodinates of all labelled cells in 40 to 50 image fields (depending upon the size of the tumour).

Geostatistics
The theory of geostatistics was invented by Matheron [21,22]. It is a corpus of mathematical/statistical methods based upon random models on limited spatial fields. It has seldom been used in biology. The most common applications of geostatistics are in mining (especially for precious ore and diamonds), for reserve estimation [9,18].
We demonstrate here that geostatistics may be used to determine grading classes by bypassing heterogeneity problems in breast cancer. The distribution of cells positive for a given marker in a breast tumour specimen can be considered as a unique, partially unknown, realisation z(x) of a spatial random process Z(x). Grading is in fact a geostatistical inference based upon a unique (single) realisation [21] of the histological spatial process. In this case, the spatial process Z(x) is the point wise process defined by the existence at the point x of a cell positive for a given marker. The method used here is based on an approach through spatial and point  statistics [12,29] and dispersion variance, which stems from geostatistics.
It is important to understand that each of our specimens can be viewed as a unique realisation z(x) of the stochastic process Z(x), and that this forbids most classical statistical inferences. Stationarity (see Appendix) helps bypassing in part uniqueness difficulties.
Only two-dimensional (2D) space will be considered here, so as to simplify our work, although a generalisation to 3D would be straightforward, as the most plausible hypothesis in large fields is isotropy [19].
Let a complete section V , of arbitrary shape, be partitioned into disjoint sub-domains v 1 , . . . , v k , with equal-form and size |v| (Fig. 2).
Let the dispersion variance over the partition V [9, 28] be whereẑ(v i ) is the spatial average of the process z(x) which is a non-biased estimator of the point process mean µ = E{Z(x)}. In this study,ẑ(v i ) represents the density of (positive) marked cells over sub-domain v i . The dispersion variance reflects the dispersion of the density of marked cells over sub-domains of a same size |v| around the overall density (ẑ(V )) and over the whole section V . As |v| and |V | become larger, the dispersion variance (s 2 (v|V |)) should converge to zero. The speed of convergence toward zero can be used as an indicator of heterogeneity: the more the process is homogeneous, the faster the dispersion variance decreases. If the point process Z(x) is homogeneous at all observation scales, as in the case of a Poisson point process, s 2 (v|V ) converges to zero as K/|v|, where K is real positive valued [12]. An ergodic process shows a similar asymptotic behaviour [20], homogeneous for observation scales greater than the integral range (see Appendix) [20]. On the other hand, if no convergence towards zero occurs, the spatial average loses its meaning and can no longer describe the process.
Following the conclusions of the Kananaskis Tumour Heterogeneity working group [8], the asymptotic behaviour of the dispersion variance is used here to assess heterogeneity. A Poisson point process, with an asymptotic convergence to zero as K/|v|, can be viewed here as a null model, that represents homogeneity at all observation scales.
A linear model for the asymptotic behaviour of the dispersion variance can be obtained. Let a log-log graph of s 2 (v|V ) versus |v| show the convergence of the dispersion variance towards a linear alignment with slope −α (0 < α 1) and intercept β, i.e., log meaning that the dispersion variance behaves asymptotically as whereσ 2 is the estimation of the point process Z(x) and variance (Var{Z(x)}) is obtained using the following dispersion variance property: The positive exponent α in Eq. (4) measures the speed of convergence dispersion variance toward 0 and K = e β represents the integral range [20] when the point process Z(x) is ergodic (see Appendix).
The asymptotic exponent of the dispersion variance, α, reflects the speed of convergence of the spatial averageẑ(v) to the statistical mean µ = E{Z(x)} of the spatial point process Z(x). The more the process is heterogeneous, the slower the convergence of the dispersion variance to zero and the smaller the value of α. The asymptotic exponent α is used to assess heterogeneity.

Practical implementation
The following procedure was applied to estimate the asymptotic exponent α and the offset β. Domain V may be tiled in any shape by design (Fig. 2); however, the use of rectangular sections allowed partitioning and recombining in an easier way. The size of the domain V varied between a minimum of 4000 µm × 800 µm and maximum of 6000 µm × 1200 µm depending on the size of the tissue block. Progressive partitioning of V was used with 2 i , i = 2, . . . , 14 sub-domains. This led to a maximum of 16384 sub-domains. This maximal partition depth was chosen such that at least a couple of cells not necessarily labelled subsist in the smallest sub-domain.
The density of marked cells inẑ(v) in each subdomain v was obtained from the stored coordinates.
Graphs of s 2 (v|V ) from all specimens were shown and the dispersion variance was normalised so as to allow all curves to begin at the same point and the size of the sub-domain |v| to be scaled by the size of the initial partition |v| min = |V |/2 14 .
The values of α and β were estimated by linear regression over the last part of the log-log s 2 (v|V ) versus |v| curve, where linear alignment occurs. Clearly, for a homogeneous Poisson process, α converges rapidly toward 1. The log-log s 2 (v|V ) versus |v| graph is then compared to a standard line with slope −1, representing the behaviour of a homogeneous Poisson process.
The robustness of the estimation of α versus the block location of the tumour has been analysed by applying the procedure above to two blocks of from the same tumour. The values in (Table 2) show that the estimated α is not dependent of location of blocks into the tumour. However, a more complete analysis of robustness will require a confidence interval estimation. The estimation of the integral range using K = e β is not robust. A slight estimation error in β can degenerate into a large estimation error for K. More generally, the values of E{e X } and e E{X} might be far from one another.

Results
The stationarity condition was verified by comparing the curve of s 2 (v|V ) obtained over different subdomains of the same sample V .
All graphs show a similar behaviour for small |v| (the slope of log-log dispersion variance graph is equal to −1) which is compatible with homogeneity at finest scale of observation. For coarser scales, the behaviour of the dispersion variance graphs follows two possible paths: grade 1 behave in a linear log-log way, but higher grades yield slightly concave curves (Fig. 3).
For all specimens, linear alignment occurred over the five last points of the dispersion variance graph (which represents 2 2 to 2 7 subdomains). These points were used to estimate α and β by linear regression.
The estimated asymptotic exponents of dispersion variance (α) for grade 3 samples were all lower than 0.5 which means a considerable level of heterogeneity. On the other hand, grade 1 and some grade 2 specimens showed a high level exponent (near 1), which is compatible with a relative homogeneity.
Parameter (α) is a measure of heterogeneity; it was used to discriminate histopathological grades (Fig. 4).
The measure of heterogeneity α separates SBR grades 1 and 2 from grade 3 (Mann-Whitney U test, p 0.05) ( Table 3). However, grade 1 and grade 2 are not significantly separated. Two kinds of behaviour of the dispersion variance can be detected for grades 1 and 2: (i) a completely linear curve with high slope (α > 0.75), which is compatible with a homogeneous process at the scale of observation; (ii) a semi-linear The failure of the Mann-Whitney test to separate between grades 1 and 2 might mean that some tumours of grades 2 (4 cases out of 10) assessed as grade 2 by the pathologist are in fact of grade 1, and that one assessed as grade 1 (1 out of 4) was in fact of grade 2. More precise analyses will require a larger set of data.

Discussion and conclusions
The samples were used without any possibility of bias (systematic morphometry) and the procedure is easily reproducible; further image processing should allow full automation.
The most important result is that there exists a geostatistical parameter (the asymptotic exponent of the dispersion variance α) that can be used as a measure of spatial heterogeneity and grading. Parameter α shows a progressive, statistically significant increase from grades 1 to 3 passing by grade 2 (Fig. 4).
An objective automatic histo-pathological grading procedure using the asymptotic exponent of the dispersion variance (α) should eventually be applicable in a clinical environment.
The dispersion variance α clearly bears a relationship with heterogeneity.  Our conjecture about grading is that grading and heterogeneity are akin. The full confirmation will require a confrontation between our grades and the patients' outcome after at least five years of follow up.

Acknowledgements
We wish to thank Dr. Christian Lantuéjoul (Ecole des Mines, Centre de Géostatistiques) for his help and remarkable expertise in geostatistics and the nice discussions we had together. We thank the referees, notably one for pointing out a reference [13] and Dr. Angela Downs-Rigaut for her help in revealing some unclear points.

A.1. Stationarity
Stationarity is defined by the invariance by translation of the distribution law of the spatial stochastic process Z(x). This property validates the assumption that the observed value z(x) at each point is a sample of a single random variable. Therefore it is sufficient to analyse only a subpart of the whole natural phenomenon which is most often too large to be studied fully. A sampling model must therefore be used to make sure of the representativity of the sample analysed. Stationarity is a strong property here, as it concerns directly the distribution law of any tuples of values of Z.
A weaker stationarity condition requires only that the expectation of punctual and doublets of values of the function exist and are invariant by translation: In practice, we always found this sufficient to ensure stationarity [20].

A.2. Ergodicity
For making statistical inferences, the stationarity condition is necessary but not sufficient as the observed values of Z(x) are not statistically independent. Ergodicity ensures that an estimation of statistical parameter of the stochastic process based on spatial averages is consistent, i.e., that spatial averages converge to the statistical expectation if the sample size is sufficient.
Let the statistical mean of a stochastic process µ = E{Z(x)} be estimated by the spatial average on a sample |V | (z (V )): This estimator is unbiased. For an ergodic function, the variance of the estimator can be approximated for large sample (|V | large) by: where A is the integral range [20].

A.3. Integral range
When |V |/A = N , the variance estimator will behave as the variance of an estimator based upon N independent observation of variance σ 2 (var{z(V )} ≈ σ 2 /N ). The integral range can be interpreted as the scale of correlation of an ergodic process Z(x).