Assessment and Visualization of Ki67 Heterogeneity in Breast Cancers through Digital Image Analysis ()
1. Introduction
Breast cancer is the most common malignancy and leading cause of cancer-related fatalities among women worldwide [1] . In addition to histological grading and staging, immunohistochemical markers, including hormone receptors HER2 and Ki67, are essential in breast cancer diagnosis [2] . These biomarkers not only serve as treatment indicators but also offer valuable prognostic insights. Ki67, a nuclear protein expressed in proliferative cells, serves as the foundation for determining the Ki67 index (KI), denoting the percentage of immunoreactive tumor cells among all tumor cells, and thereby serving as a measure of tumor proliferation [3] . The KI has been used for subtyping hormone-positive breast carcinomas into luminal type A (low KI) and luminal type B (high KI), which predict prognosis and indicate suitability for chemoradiotherapy [4] [5] [6] . Typically, a threshold of 14% is adopted to demarcate low and high proliferation categories. Therefore, it is critical to precisely quantify the KI. However, in daily clinical practice, KI assessment predominantly relies on the subjective estimation of pathologists. Despite meticulous adherence to the standard international protocol during counting, substantial intra- and inter-observer variations persist [6] [7] .
Many factors contribute to observer variation [6] , including pre-analytical considerations (specimen type and cold ischemic time), analytical aspects (antibodies and autostainers), and interpretation factors (scoring method and area chosen). Pathologists often choose different areas for KI counting due to intratumoral heterogeneity. Ki67 intratumoral heterogeneity presents a challenge for assessment and quantification and is a daily phenomenon encountered by pathologists. In most breast cancers, Ki67-positive cells exhibit non-uniform distribution throughout the tumor, with peripheral tumor cells often displaying greater activity and higher KI. Some studies reported that hotspots tend to occur most frequently at tumor edges [8] . Conversely, to investigate the predictive value of Ki67 heterogeneity on prognosis, researchers have measured the KI gaps between hotspots and the average from sampled fields [9] [10] [11] . However, evaluating this heterogeneity has proven challenging due to the inability to manually obtain comprehensive KI informationacross the entire tumor face.
In recent decades, due to advancements in deep learning and convolutional neural networks, digital image analysis (DIA) has found applications in automatic KI assessment. Both commercial and open-source image analysis software solutions have become readily accessible [12] [13] . Compared to human evaluation, these algorithms are more efficient and reliable in quantification [14] [15] . Numerous pathologists and researchers have used these tools to study KI in breast cancer patients. However, most studies choose to evaluate KI in limited regions (tissue microarray or hotspots), serving as proxies for the whole tumor [16] [17] . This approach is favored because a full-face tumor whole-slide image (WSI) typically includes thousands of high-power fields (HPF). The use of WSIs would consume substantial computational power and time for analysis. In 2016, researchers applied DIA to breast cancer Ki67 WSIs to obtain a comprehensive Ki67 index. They also introduced sophisticated indicators to describe Ki67 heterogeneity [18] [19] . Retrieving complete Ki67 profiles within WSIs, as opposed to obtaining a single Ki67 index number, is the first step in decoding heterogeneity.
In order to survey how Ki67 heterogeneity affects breast cancer subtyping and how many areas are adequate for Ki67 index estimation, in this study, we used DIA to assess comprehensive Ki67 profiles in breast cancer at the WSI level and to visualize intratumoral heterogeneity. In addition, we utilized this data to simulate the impact of sample size and heterogeneity on KI estimation.
2. Materials and Methods
2.1. Case Collection
We collected data from 50 excised breast cancer cases that were diagnosed in the Department of Anatomical Pathology at Far Eastern Memorial Hospital, New Taipei City, Taiwan area, between January 2020 and December 2020. Inclusion criteria encompassed only cases of breast invasive carcinoma of no special type. Mucinous carcinoma, metaplastic carcinoma and lobular carcinoma were not included. All patients are over 35-year-old and include histological grading 1 to 3. Two cases were excluded due to the invasive area being less than 10 HPF. Consequently, our study analyzed a total of 42 cases. The clinical and pathological details of the patients are shown in Table 1.
![]()
Figure 1. Ki67 index analysis and sampling simulation. (A) The tumor area in Case 24 was annotated, and the software cropped the region into tiles. The neural network identified positive tumor cells (PTC, red), negative tumor cells (NTC, blue), and stromal cells (yellow) in tile images and calculated the Ki67 index. (B) Tiles containing more than 100 tumor cells were included, and the Ki67 data was summarized to generate a boxplot, histogram, and heat map. (C) The KI data was also used in sample simulation to generate a hit rate probability histogram for different sample sizes. The red line represents global Ki67, with the target range falling between the two green lines.
2.2. Ethical Approval
This study received approval from the Research Ethics Review Committee of the Far Eastern Memorial Hospital (No. 110012-E). All slides used in the clinical diagnosis setting were retrieved from the repository. These slides and images were kept anonymous and did not contain any personal information. Therefore, the requirement for informed consent was waived.
2.3. Specimen Preparation and Ki67 Staining
Ki67-immunostained slides, utilized in a routine diagnostic context, were sourced from the repository of the Anatomical Pathology Department at the Far Eastern Memorial Hospital. The specimens were fixed in 10% neutral-buffered formalin 6 for 72 h and subsequently embedded in paraffin. The tissues were cut into 5 μm-thick sections and mounted on hydrophilic slides. The Ki67 immunostaining procedure was performed using a Benchmark Ultra Automated Staining System (Ventana Medical Systems, Tucson, AZ, USA). The slides were heated to 96˚C for 34 min to facilitate antigen retrieval and then subjected to incubation with the Ki67 antibody (SP6, Biocare Medical, Pacheco, CA, USA) at a 1:200 dilution. To visualize the immunoreaction, the ultraView Universal DAB Detection Kit (Ventana Medical Systems) was used, with hematoxylin serving as the counterstain. The Ki67 immunostained slides were scanned using a Hamamatsu S210 microscope (Hamamatsu Photonics, Hamamatsu, Japan) at 40× magnification (resolution of 0.23 μm per pixel), and the WSIs were saved in NDPI format.
2.4. Ki67 Index Analysis
We employed SmartPathology software (ver. 1.2.0, Quanta Computer, Taipei, Taiwan area) for the WSI KI analysis (Figure 1). First, a licensed pathologist (CMH) identified the tumor and captured a rectangular area containing the invasive component. The software then automatically segmented the field into image tiles (1936 × 1216 pixels) and saved each tile in PNG format with sequential numbering as filenames. The real size of each tile is 0.124 mm2, equivalent to one HPF [20] . Each case comprised 207 - 4215 tiles (1924 tiles on average). These tiles were analyzed using a Breast Cancer Model that identified cells and classified them into three categories: positive tumor cells, negative tumor cells, and stromal cells. The Breast Cancer Model was a mask region-based convolutional neural networks algorithm. For evaluating Ki67 index in tumor area, the interclass coefficient between algorithm and manual counting could reach 0.99. However, for non-tumor or lymphocyte-rich area, the algorithm might misrecognize stomal cells as tumor cells. The KI was defined as the percentage of positive tumor cells among the total tumor cells. After all the tiles were inferred, the data from these tiles were compiled into a CSV file for each case. When computing the global KI (the percentage of positive tumor cells across the entire tumor), we excluded tiles with fewer than 100 tumor cells. Due to the presence of outlier data (tiles exhibiting extremely high KI), we defined hotspots as tiles featuring KI values at the 90th percentile (Pareto hotspots), rather than focusing only on those with the highest KI [18] . For subtyping, we adapted the Saint Gallen International Expert Consensus on the primary therapy of early breast cancer in 2011, they suggested 14% of KI as threshold to categorize tumors into two groups: low proliferative (<14%) and high proliferative (≥14%).
2.5. Sampling Simulator
To investigate how intratumoral heterogeneity impacts the accuracy of KI sampling, that is, the number of tiles (or HPF) that pathologists should select to approach the global KI, we used a Python-based sampling simulator. In this simulator, varying numbers of tiles (n = 1, 2, 3…) were chosen to estimate the global KI, with a predefined target range of global KI ± 2.5%. If the average KI of the sampled tiles, KI (n), fell within the target range, it was classified as a “Hit”. The program conducted 10,000 iterations of the sampling experiment for each sample tile number (n), allowing us to compute the probability of achieving a “Hit” for each n, which we termed the “Hit rate (n)”. Using Case 24 as an example (Figure 1), when we randomly selected one tile (n = 1) to estimate the global KI, we found that in 1465 out of 10 000 tests, KI (1) fell within the target range, yielding a hit rate (1) of 0.1465. When the sample size was increased to 128 (n = 128), the hit rate reached 0.9519, which was higher than 0.95. We repeated the simulation until the hit rate exceeded 0.95, at which point we recorded the sample size as n = 95. Theoretically, the higher the KI heterogeneity, the larger the number of sampling tiles required to achieve a 0.95 hit rate, leading to a higher value for n = 95.
2.6. Statistical Analysis
Statistical analyses were performed using SPSS version 21 for Windows (Chicago, Armonk, NY, USA). To statistically and spatially represent the KI distribution, we employed various visualization techniques, including boxplots, histograms, and heat maps, which were generated using the Seaborn plugin (v.0.12.1) in Python (ver. 3.6.10).
3. Results
3.1. Cases Information
The details of the 42 cases summarized in the KI data are presented in Table 2 and Figure 2. The global KI ranged from 2.50% to 42.79%, while the hotspot KI ranged from 5.92% to 88.46%. Concerning the global KI, 23 cases were categorized as low-proliferative tumors (<14%), whereas only 10 cases fell into the low-proliferative group based on hotspot KI. Therefore, 30% (13/42) of the cases showed discrepancies between the different classification methods. In addition, we calculated cellularity (tumor cells per tile), which ranged from 172 to 490 per tile.
* means Cases with different KI tiers by global KI and hotspot KI.
![]()
Figure 2. Global KI and Hotspot KI of all cases. Using the global KI, 23 cases were identified as low proliferative tumors (<14%), and of those, 13 cases were reclassified as high proliferative when assessed by the hotspot KI.
![]()
Figure 3. Ki67 profiles and histograms. (A) The boxplot displays the Ki67 index distribution for all cases. (B) The KI histogram for Case 21 exhibits a bimodal distribution and heterogeneity. (C) The histogram for Case 1 shows a homogenous distribution.
![]()
Figure 4. Hotspots and Ki67 distribution patterns. (A) Hotspot in Case 12. (B) Hotspot in Case 18. Some red blood cells and lymphocytes (arrows) were mistakenly identified as positive tumor cells (red) and negative tumor cells (blue). (C) The heat map for Case 22 shows a homogenous randomly-distributed pattern. (D) The heat map for Case 35 shows a heterogenous randomly-distributed pattern. (E) The heat map for Case 10 shows a peripheral pattern. (F) Case 33 shows central scar with low activity. (G)-(H) Case 21 and Case 16 feature central hotspots.
![]()
Figure 5. Correlation between n (95) (A) and global KI and (B) hotspot KI.
3.2. Statistically Heterogenicity
In Figure 3(A), a boxplot shows the distribution of intratumoral KI values in 42 cases. Some patients exhibited a wide interquartile range (IQR), as seen in Cases 13, 21, and 35, indicating a higher degree of heterogeneity in KI levels. In general, patients with a higher global KI also had a higher interquartile range. The histogram representing Case 21 (Figure 3(B)) exhibits a heterogeneous, bimodal distribution characterized by significant variation in KI values. In contrast, the histogram for Case 1 (Figure 3(C)) displays a homogenous distribution, with most tiles exhibiting similar KI values with narrow variations. Beyond the IQR, we observed outlier data points with extremely high KI values, most of which were attributed to the software incorrectly identifying red blood cells or lymphocytes as positive tumor cells (Figure 4(A) & Figure 4(B)). Therefore, we opted to use the 90th percentile KI (Pareto hotspot) instead of the highest KI to represent the hotspot KI.
3.3. Spatial Heterogenicity
We used a heat map to elucidate the spatial distribution of KI, which can be categorized into three patterns based on the location of hotspots: random, peripheral, and centered (Figures 4(C)-(H)). Among our cohort of patients, nearly half (52%, 22/42) had a randomly distributed Ki67 density. The peripheral type (43%, 18/42) showed increased proliferation activity at the tumor edges. Only two patients (5%, 2/42) displayed hot central zones.
3.4. Sampling Simulation
The results of the sampling simulation are presented in Table 2, with n (95) ranging from 4 to 373 (67 on average). Therefore, in most cases, pathologists must sample more than 10 tiles to estimate global KI. Generally, in accordance with our assumption, cases with greater heterogeneity (wider IQR) had a higher n (95). One exception to this trend was seen in Case 25, which displayed a low global KI and IQR but exhibited the highest n (95) among the 373 tiles. This result may have been caused by the presence of a large number of outlier data points. In addition, n (95) demonstrated a positive correlation with global KI and hotspot KI (Figure 5(A), Figure 5(B)).
4. Discussion
The Ki67 index has been used in the field of breast pathology for several decades, and Ki67 intratumoral heterogeneity is a well-known issue. Nevertheless, numerous unresolved questions stem from this heterogeneity, including whether global KI or hotspot KI should be reported. Where should we define hotspots, and how do we locate them? Which threshold should be used to divide tumors into low- or high-proliferative groups? Does the degree of statistical or spatial heterogeneity matter, and how should we measure it? How many HPF or tumor cells should pathologists count? These questions pose a considerable challenge for resolution, primarily due to the limitations faced by human pathologists in quantifying millions of cells or conducting a comprehensive evaluation of a single Ki67 slide. Due to the advancements made in deep learning and computational capabilities, researchers are now able to apply DIA to WSI to address these questions. In this study, we used image analysis software to infer full-face tumor Ki67 slides and retrieved Ki67 profiles to generate boxplots and histograms that statistically visualize heterogeneity. Heat maps were employed to reveal the spatial information and different Ki67 distribution patterns. Based on sampling simulations, our findings suggest that more than 10 HPFs (approximately 1.24 mm2) are required for global KI estimation.
Due to the intratumoral heterogeneity of Ki67 in breast cancers, global KI and hotspot KI are correlated but significantly different. As noted in this study, setting the threshold at 14% resulted in a reclassification of one-third of cases, upgrading them to a highly proliferative category based on hotspot KI instead of global KI. For most pathologists, the Ki67 index has a similar mitotic count; therefore, we tended to evaluate focal hotspots rather than the whole tumor. Some researchers posit that the most proliferative areas represent tumor behavior and predict prognosis more accurately [8] . However, most early studies on KI in breast cancer primarily relied on tissue microarrays (random samples) or genetic profiling rather than a comprehensive examination of full-face tumor sections [16] [21] [22] . Therefore, the use of global KI may align more consistently with the results of these studies. Regardless of whether global KI or hotspot KI is employed, the primary challenge lies in the manual impossibility of conducting a thorough, full-face tumor Ki67 evaluation; pathologists can only sample some areas and count a limited number of cells (up to 1000). These comprehensive Ki67 profiles extracted from WSI can provide more information than a single percentage value, whether global a hotspot-derived.
Hotspot identification is another often overlooked problem in pathology due to the absence of a well-defined practical framework. Intuitively, pathologists rely on subjective judgments to identify areas with the highest Ki67 positivity when observing low-power fields, designating these regions as hotspots [11] [23] . However, there is no consensus regarding the optimal size of such areas and the minimum number of cells required to qualify as a hotspot. In our study, we established an image tile, similar to an HPF, as the fundamental unit for dividing the WSI and chose 100 tumor cells as the minimal requirement. Statistically, we could identify the tiles with the highest KI from the Ki67 profile, but these might be a few outliers that do not represent tumor behavior. Therefore, the concept of “Pareto hotspots” was adopted [18] , and we considered that these areas might offer a more accurate reflection of tumor proliferation compared to sporadic, extreme hotspots. Spatially, we easily identified hotspots within the heat map, which clearly visualizes the Ki67 density.
The next question pertained to determining the most appropriate threshold for clinical use. Many studies have proposed the separation of low- and high-proliferation groups as follows: 10%, 14%, or 20% [4] [5] [24] . Recently, the Gallen/Vienna Consensus Conference on Early Breast Cancer Treatment Standards has recommended three categories: low (<5%), intermediate (5% - 30%), and high (30%) [25] . In addition to pre-analytical factors, inter-laboratory, inter-observer, and even intra-observer variations hinder the threshold setup. Therefore, a more precise and reliable quantification method is required. Using DIA, we can not only obtain a comprehensive Ki67 profile, but the data also exhibit greater accuracy and reproducibility compared to human estimates [8] [14] . The precise quantification of KI can solidify the Ki67 study, enabling the establishment of a threshold for patient classification and even considering KI as a continuous variable [2] .
Some researchers have suggested that Ki67 intratumoral heterogeneity could be a prognostic factor, and various parameters have been used to indicate the degree of heterogeneity. The most commonly used parameter is known as the “heterogenicity score”, typically defined as the difference between the hotspot KI and the global KI (or lowest KI). Two studies found that heterogeneity scores were independently associated with prognosis and lymph node metastasis [10] [23] . However, the KI data in these two studies were obtained through manual counting in the selected fields. Plancoulaine et al. [18] used the comprehensive Ki67 profile obtained through DIA and proposed parameters to quantify heterogeneity, including entropy and bimodality. Their findings revealed that Ashman’s D, an indicator of bimodality, was an independent predictor of overall survival [19] . Furthermore, Ki67 expression was quantified using the hexagonal tiling approach, and heat maps were generated to depict spatial heterogeneity. Taking inspiration from these trailblazing studies, we employed heat maps to visualize Ki67 density and initially classified the spatial distribution into three patterns: randomly distributed, peripheral, and central. The peripheral pattern conforms to the general assumption held by pathologists: the most proliferative cells are situated at the tumor edges, where they actively search for mitotic activity. A classic example of this pattern can be observed in tumors with central necrosis, as they are so large that only peripheral tumor cells have sufficient blood supply to grow. However, nearly half of the tumors exhibited a randomly distributed pattern, indicating that the hotspots were evenly distributed within the tumor. Only a minority of patients displayed proliferative centers. These different patterns could be the result of differences in factors such as blood supply, immune responses, or the tumor microenvironment (hormones or growth factors). This spatial information provides a novel perspective for interpreting Ki67 expression and understanding tumor behaviors.
In the assessment of KI, the matter of sample size is another long-standing question. Early studies suggested a requirement of 500 - 1000 cells [21] . According to the protocol established by the International Ki67 Breast Cancer Working Group, a minimum of three HPFs (including more than 100 cells) from different Ki67 density areas is required for estimation. In our study, we conducted simulations with various sample sizes to estimate the global KI, and most cases required more than 10 HPFs. Reasonably, a higher global KI value corresponds to increased heterogeneity and a greater number of tiles (HPF) to be examined [26] . In clinical practice, it could be imprecise to use KI from core biopsy to represent KI in whole tumors, especially in cases of tumors with a high global KI. Therefore, it is necessary to repeat the Ki67 evaluation when excised specimens are available [27] . Nevertheless, it is imperative to establish a comprehensive Ki67 profile, given that sampled KI values from limited areas cannot represent tumors with high heterogeneity.
There are several limitations to this study, the first of which is the absence of a correlation test between Ki67 and genetic profiles. The question of which parameters within the Ki67 profile play a more decisive role in predicting genetic subtypes remains unanswered. A technological challenge we encountered was outlier data owing to an imperfect algorithm. The software incorrectly identified lymphocytes or red blood cells as tumor cells. However, our team is working on enhancing and accelerating KI assessment by applying tumor masks [12] [28] . Another concern is the coexistence of ductal carcinoma in situ within the invasive components. In theory, KI should only be counted in the invasive regions; however, we did not exclude in situ lesions. This issue could potentially be addressed through area-segmentation algorithms capable of distinguishing invasive lesions from in situ lesions.
In next step, we are planning to solidify the result of study in advance. Firstly, epithelial masks will be applied to KI quantification for more accurate evaluation. Also, the statistic analysis and heat map generation will be performed automatically and packaged into the software. Therefore, data generation will be more efficient. Secondly, more breast cancer cases are required and subtypes, such as luminal type or triple-negative type, should be separated into different group. This is because clinical utility of KI might work differently in subtypes. For example, in hormone-positive early breast cancer, KI is an indicator for adjuvant chemotherapy, but for triple-negative breast cancer, it could be a prognostic factor. Lastly, overall or disease-free survival rate should be correlated with their comprehensive KI profiles and KI heat map patterns. Therefore, it is possible to figure out which parameters and thresholds have the most important clinical utility.
In conclusion, using DIA, we are now capable of addressing Ki67 heterogeneity and revealing distribution patterns that were previously deemed implausible. Future studies should focus on more robust software to consistently generate spatial information in a wider range of cases, enabling the determination of parameters that accurately represent tumor behavior. The comprehensive Ki67 profile holds the potential to revolutionize the approach to conducting proliferative activity surveys in breast cancers.
Data Availability
The datasets used and analyzed in the current study are available from the corresponding author upon reasonable request.
Acknowledgements
The authors thank Dr. Chu Hsiu-Yi for her support and assistance in writing.
Author Contributions
Chien-Hui Wu: Methodology, Investigation, Writing—Original Draft. Min-Hsiang Chang: Methodology, Investigation, Software, Formal analysis, Writing—Original Draft. Hsin-Hsiu Tsai: Software. Yi-Ting Peng: Software.
Abbreviations
DIA, digital image analysis;
HPF, high-power fields;
IQR, interquartile range;
KI, Ki67 index;
WSI, whole-slide image.