Intratumoral Heterogeneity and Immune Response Indicators to Predict Overall Survival in a Retrospective Study of HER2-Borderline (IHC 2+) Breast Cancer Patients

Breast cancer (BC) categorized as human epidermal growth factor receptor 2 (HER2) borderline [2+ by immunohistochemistry (IHC 2+)] presents challenges for the testing, frequently obscured by intratumoral heterogeneity (ITH). This leads to difficulties in therapy decisions. We aimed to establish prognostic models of overall survival (OS) of these patients, which take into account spatial aspects of ITH and tumor microenvironment by using hexagonal tiling analytics of digital image analysis (DIA). In particular, we assessed the prognostic value of Immunogradient indicators at the tumor–stroma interface zone (IZ) as a feature of antitumor immune response. Surgical excision samples stained for estrogen receptor (ER), progesterone receptor (PR), Ki67, HER2, and CD8 from 275 patients with HER2 IHC 2+ invasive ductal BC were used in the study. DIA outputs were subsampled by HexT for ITH quantification and tumor microenvironment extraction for Immunogradient indicators. Multiple Cox regression revealed HER2 membrane completeness (HER2 MC) (HR: 0.18, p = 0.0007), its spatial entropy (HR: 0.37, p = 0.0341), and ER contrast (HR: 0.21, p = 0.0449) as independent predictors of better OS, with worse OS predicted by pT status (HR: 6.04, p = 0.0014) in the HER2 non-amplified patients. In the HER2-amplified patients, HER2 MC contrast (HR: 0.35, p = 0.0367) and CEP17 copy number (HR: 0.19, p = 0.0035) were independent predictors of better OS along with worse OS predicted by pN status (HR: 4.75, p = 0.0018). In the non-amplified tumors, three Immunogradient indicators provided the independent prognostic value: CD8 density in the tumor aspect of the IZ and CD8 center of mass were associated with better OS (HR: 0.23, p = 0.0079 and 0.14, p = 0.0014, respectively), and CD8 density variance along the tumor edge predicted worse OS (HR: 9.45, p = 0.0002). Combining these three computational indicators of the CD8 cell spatial distribution within the tumor microenvironment augmented prognostic stratification of the patients. In the HER2-amplified group, CD8 cell density in the tumor aspect of the IZ was the only independent immune response feature to predict better OS (HR: 0.22, p = 0.0047). In conclusion, we present novel prognostic models, based on computational ITH and Immunogradient indicators of the IHC biomarkers, in HER2 IHC 2+ BC patients.


INTRODUCTION
Breast cancer (BC) is a complex and diverse disease with distinct clinical, pathological, and molecular characteristics. The multifaceted nature of the disease leads to diverse clinical outcomes and therapeutic responses. BC has been classified into several biologically distinct subtypes: luminal A, luminal B, human epidermal growth factor receptor 2 (HER2)-enriched (HER2), basal-like, and normal-like by gene expression profiling analysis (1,2), requiring different treatment strategies. This categorization of the BC subtypes has been adapted for clinical practice and is mainly based on immunohistochemistry (IHC) assessment of estrogen receptor (ER), progesterone receptor (PR), HER2, and Ki67 expression.
Routinely used predictive features, including clinicopathological parameters (age, tumor size, lymph node status, and histological grade) and biomarkers (ER, PR, and HER2) are insufficient for personalized clinical decisions in BC patients (3). Novel prognostic BC biomarkers have been intensively investigated as recently reviewed by Wu et al. (4). In particular, robust biomarkers are in demand for HER2-positive disease to improve selection of patients for current and emerging therapies of HER2-positive metastatic BC (5) as well as for prediction of resistance for anti-HER2 therapies, recurrence (6,7), and particular consequences of the disease (8). Novel approaches based on pathology image analytics and machine learning methods open new perspectives for predictive modeling and clinical decision support (9,10). Importantly, both molecular and image-based biomarkers can be explored and validated using The Cancer Genome Atlas (TCGA) Data Portal (11).
HER2 amplification and overexpression occur in approximately 15%-20% of invasive BC cases and are associated with worse patient survival as compared with nonamplified HER2 BC (12)(13)(14)(15). A positive HER2 status predicts better effect of HER2-targeted therapies, and therefore, its accurate detection is essential for treatment decisions (16,17).
While the majority of tumors can be categorized as either HER2-positive or HER2-negative by IHC and in situ hybridization (ISH) techniques, which are regarded as the standard methods to assess HER2 status in BC, borderline tumors do account for up to 18% of BCs (18,19) and present challenges for patient assessment and therapy choices. In 2018, the American Society of Clinical Oncology (ASCO) and College of American Pathologists (CAP) updated the guidelines for HER2 testing with revised criteria for HER2 IHC borderline (IHC 2+) classification. This mainly focused on less common fluorescence ISH (FISH) patterns (ASCO/CAP groups 2, 3, and 4) and recommended to integrate them with a concomitant IHC review for a final HER2 result determination (20). The ambiguous FISH equivocal group (18), which poses therapeutic dilemmas, was removed, which resulted in an increased frequency of HER2-negative cases (21)(22)(23)(24). Nevertheless, some studies report that HER2 equivocal tumors present similar clinical behaviors to the HER2-negative BC (25,26), while others find differences in clinicopathological and prognostic aspects between these two categories (23). This suggests that the equivocal category represents an intermediate state between HER2-positive and HER2-negative tumors (27,28).
Approximately 15%-30% of IHC 2+ cases are HER2amplified (29), while the remaining IHC 2+ and IHC 1+ HER2 non-amplified tumors were recently designated as a relatively common "HER2-low" category, accounting for approximately 40%-55% of BC (30)(31)(32). This concept becomes important with the advent of a new generation of anti-HER2 agents. Specifically, ongoing clinical trials have demonstrated high efficacy of antibody-drug conjugates that are designed to target and deliver chemotherapy inside cancer cells in this particular subset of BC patients (33)(34)(35). The HER2-low BC group is not formally defined at present, but if treatment options will become available, the current dichotomous HER2 guidelines will have to be revised further to distinguish truly HER2-negative from HER2-low breast cancer (31).
The intratumoral heterogeneity (ITH) of HER2, at both protein expression and gene amplification levels, is a common feature of HER2-borderline tumors, which further complicates the assessment of HER2 status (36)(37)(38)(39)(40). In addition to the heterogeneous HER2 expression, the variable expression of hormone receptors (HRs) also contributes to the ITH and may further affect clinical outcomes and responses to treatment of BC (41,42). Potential interactions between HR and HER2 signaling pathways, which could impact development of resistance to endocrine and anti-HER2 therapies, have been highlighted by several preclinical and clinical studies (43)(44)(45)(46)(47)(48).
Current pathological IHC methods are based on the assessment of a proportion of HER2-positive tumor cells; however, ITH of HER2 expression may present a challenge in some tumors to be categorized with a single value (0, 1+, 2+, and 3+). Digital image analysis (DIA) has opened new opportunities in HER2 IHC assessment by providing biomarker quantification with increased accuracy, precision, reproducibility, and capacity (49)(50)(51)(52)(53)(54)(55). Studies have demonstrated that DIA can reliably distinguish HER2 IHC negative (0-1+) and positive (3+) cases and reduce the proportion of IHC 2+ cases (51)(52)(53)(54). Importantly, continuous data and spatial aspects of IHC biomarker distribution can be revealed by DIA (56)(57)(58). Several diversity metrics (the Shannon entropy (59), the Simpson index (60), and Rao's quadratic entropy (61) have been adapted for molecular, genetic, and microenvironmental heterogeneity assessments in BC) (62)(63)(64)(65). Potts et al. examined HER2 expression ITH in BC by combining semiquantitative analysis with ecology diversity statistics (64). Both cell-level and tumor-level heterogeneities were evaluated, but the authors had doubts about the insufficient number of regions sampled to make an assessment of heterogeneity at a tumor level. Several recent studies (56)(57)(58) showed a successful assessment of ITH of IHC biomarkers in whole slide images (WSIs) based on hexagonal grid subsampling of DIA data; importantly, this methodology enabled retrieval of prognostically informative spatial heterogeneity indicators of tissue biomarker expression.
Although ITH may challenge the efficacy of therapy, it may be also associated with favorable prognostic effects, since a greater mutational load could lead to an increased tumor neo-antigen generation that attracts immune cells and stimulate antitumor immunity (66,67). However, immunogenicity is different among BC subtypes, with generally higher mutational load, higher numbers of tumor-infiltrating lymphocytes (TILs), and higher programmed death-ligand 1 (PD-L1) expression in triplenegative and lower in HR-positive subtypes (68)(69)(70)(71). These differences may impact the efficacy of therapy with immune checkpoint inhibitors with significant responses achieved only in patients with triple-negative BC so far (72)(73)(74).
TILs have been recognized as a potential biomarker of survival on BC patients (75,76); however, their prognostic significance varies in BC types (70). A positive prognostic role of CD8 cells has been demonstrated in ER-negative and triplenegative BC (77)(78)(79), but its prognostic value in HR-positive BC remains unclear (70,80). Recent studies have shown that the distance between immune cells and cancer cells is clinically and prognostically important in BC (81)(82)(83)(84). The methods for assessing TILs and their spatial distributions have been the focus of many studies. Recently, Krijgsman et al. (85) first applied an automated deep learning approach that identifies the tumor boundary and detects CD8-positive cells in IHC images, and then they analyzed the spatial distribution of CD8 lymphocytes in ER-positive invasive BC. They found that only the SD of the CD8 density (but not the mean of CD8 density) distribution was significantly associated with better survival, hypothesizing that it reflects the contribution from local highdensity areas. In another study (84), the immune scores of cell abundance and spatial heterogeneity were quantified using a combination of fully automated H&E-stained image analysis and spatial statistics. High immune spatial scores, but not the abundance scores, were associated with poor prognosis in ER positive BC. Rasmusson et al. proposed a hexagonal grid-based methodology to automatically detect the tumor-host interface zone (IZ) and compute the immune cell density profile across the interface. The computed Immunogradient indicators provided the independent prognostic value in HR-positive breast and colorectal cancer patients (86).
In our study, we investigated ITH and immune response properties of HER2 IHC 2+ borderline BC patients with regard to their prognostic value. We utilized image DIA of IHC for ER, PR, Ki67, HER2, and CD8, with subsequent hexagonal grid analytics to extract combined prognostic overall survival (OS) models in HER2 IHC 2+ FISH-negative and FISHpositive patients.
HER2 expression was scored as 0 (no staining, or incomplete membrane staining that is faint or barely perceptible and within ≤10% of the invasive tumor cells); 1+ (incomplete membrane staining that is faint or barely perceptible and within >10% of the invasive tumor cells); 2+ (weak-to-moderate complete membrane staining observed in >10% of tumor cells); or 3+ (circumferential membrane staining that is complete, is intense, and in >10% of tumor cells) according to the 2018 ASCO/CAP guidelines (20). IHC 0 and IHC 1+ were defined as HER2 negative, IHC 2+ was categorized as HER2 borderline, and IHC 3 + was categorized as HER2 positive.

Fluorescence In Situ Hybridization
HER2 FISH was performed on FFPE sections using the PathVysion HER2 DNA probe kit and Paraffin pretreatment kit (Abbott-Vysis, Inc., Downers Grove, IL, USA) as described in detail previously (87). Briefly, 4 µm thick sections were mounted on positively charged slides and dried overnight at 56°C. Subsequently, deparaffinization, dehydration, and pretreatment procedures were performed. After the digestion with protease, the hybridization mixture containing two fluorescently labeled DNA probes recognizing the HER2 locus (17q11.2-q12) and the centromeric region of CEP17 (17p11.1-q11.1) was applied to the target tissue. Denaturation and hybridization were performed in a hybridizer (Dako Diagnostics, Glostrup, Denmark). Then slides were washed, counterstained with DAPI, and coverslipped (Invitrogen Corporation, Carlsbad, USA). The samples were analyzed using a fluorescence microscope (Zeiss, Axio Imager.Z2, Gottingen, Germany) equipped with single-pass filters for DAPI, HER2, and CEP17, under a 63× oil immersion objective. All tumors were tested routinely by dual-probe FISH assay for final HER2 classification according to the ASCO/CAP guidelines (20).

Digital Image Acquisition, Analysis, and Calculation of Indicators
For the analysis of ER, PR, Ki67, and HER2, sections were scanned using a ScanScope XT Slide Scanner (Leica Aperio Technologies, Vista, CA, USA) at ×20 objective magnification (0.5 mm per pixel); CD8 IHC slides were scanned using an Aperio AT2 Slide Scanner (Leica Biosystems, Buffalo Grove, IL, USA) at ×20 objective magnification (0.5 mm per pixel). The DIA was performed on the WSIs with HALO ™ software (version 3.0311.174; Indica Labs, Corrales, NM, USA) by three operators (RG, RA, and GR). Initially, the tissue was classified into the tumor, stroma, and background (consisting of glass, necrosis, and artifacts) by HALO AI ™ classifier. Subsequently, the HALO Multiplex IHC and Membrane algorithms (versions 1.2 and 1.4, respectively) were applied to obtain coordinates of the cells in the IHC WSI. For quality assurance, all image analysis results were approved by the breast pathologist (RG).
The automated extraction of the IZ and Immunogradient indicators is described in detail in (86). In our study, an IZ width of seven hexagon ranks (hexagon side length 65 mm) was used. CD8 cell density was calculated in both 1) the WSI stroma and tumor areas and 2) within the tumor-stroma IZ, which consists of stroma (S), tumor (T), and tumor edge (TE) aspects. Subsequently, Immunogradient indicators (center of mass (CM) and immunodrop) representing CD8 cell density profiles across the IZ were computed. The CM indicator reflects CD8 cell density increase towards the tumor within the IZ, while the immunodrop indicator reflects an abrupt decrease of CD8 cell density across the TE (IZ rank 0) from stroma (IZ rank −1) to tumor (IZ rank 1), represented by the CD8 cell density ratio between rank −1 and rank 1.

Statistical Analysis
All continuous variables were tested for normal distribution by Kolmogorov-Smirnov test and compared by two-tiled Student's t-test (for normally distributed variables) or the Mann-Whitney U test (for non-normally distributed variables). A logtransformation was applied to normalize the asymmetric distributions of immune response variables and to meet the assumptions of parametric statistical tests; they were used in oneway ANOVA followed by Bonferroni's post-hoc test for pairwise comparisons and a two-sided Welch's t-test for homogeneity of variances. Fisher's exact test was used to assess the differences in clinicopathological variables among the analyzed groups.
A factor analysis was performed using the factoring method based on principal component analysis; factors were retained based on the threshold of an eigenvalue of 1; lastly, a general orthogonal varimax rotation of the initial factors was applied.
The optimal cutoff value for each indicator was determined using Cutoff Finder (89) to test the predictions of OS. The Kaplan-Meier method was applied to estimate the OS distributions with the log-rank test to compare survival differences between the stratified groups. To assess the prognostic factors, univariate and multivariate analyses were performed using the Cox proportional-hazards models. The "best" subset of variables to be included in the multivariate Cox proportional-hazards models was identified by leave-oneout cross-validation (90). All p-values were considered significant at the <0.05 level. Statistical analyses were performed with SAS software (version 9.4; SAS Institute Inc., Cary, NC, USA); plots were generated by R (version 4.1.0).

Clinicopathological and Follow-Up Characteristics
Clinicopathological and follow-up characteristics of the HER2 non-amplified and HER2-amplified groups are summarized in Table 1. The median follow-up period was 64 (range 2-102) and 52 months (range 0.7-100) in the non-amplified and amplified HER2 cohorts, respectively. Forty-two patients died during the follow-up, including 22 (13.7%) and 20 (17.1%) in the nonamplified and amplified tumor subsets, respectively.
The HER2-amplified group revealed significantly higher histological grade (p < 0.001) and higher frequency of increased CEP17 copy number (p = 0.0002) as compared with the HER2 non-amplified group ( Table 1). Of note, 55 (34.8%) and 67 (57.3%) cases with CEP17 copy number ≥3 were detected in the HER2 non-amplified and HER2-amplified groups, respectively. No significant differences between the groups regarding the patient age, tumor stage, and node involvement were found.

Summary Statistics of Explored Indicators
Summary statistics of the variables in the HER2 non-amplified and HER2-amplified groups are presented in Supplementary  Table 1; the variance plots of the significant differences are presented in Supplementary Figure 1.
In general, expression rates of ER and PR were higher, while Ki67 was lower in the HER2 non-amplified group. No significant difference in CD8 cell density distribution between tumor and stroma areas was observed in both the HER2 non-amplified (t = 1.72, p = 0.0867) and HER2-amplified (t = 1.07, p = 0.2841) groups. Also, the mean of CD8 density within the IZ was significantly higher in the S aspect than in the T aspect in both the HER2 non-amplified (t = 6.56, p < 0.001) and HER2amplified (t = 6.17, p < 0.001) groups. The variance of CD8 cells was the highest in the S aspect, less in the TE aspect, and lowest in the T aspect of the IZ in both the HER2 non-amplified and HER2-amplified groups (p < 0.0001) (data not shown). No significant differences of CD8 cell densities neither in tumor nor stroma areas nor inside the IZ (T, TE, and S aspects) were found between the groups. ITH (higher contrast, dissimilarity, and entropy but lower energy and homogeneity) was higher only for Ki67 in the HER2-amplified group.
Factor Analysis of Immunohistochemistry, Fluorescence In Situ Hybridization, Immune Response, and Intratumoral Heterogeneity Indicators in HER2 Non-Amplified and HER2-Amplified Groups A factor analysis was performed on the combined set of DIA IHC, FISH, immune response, and ITH data and six orthogonally independent factors in each patient group were extracted. The patterns of the factors are plotted in Supplementary Figures 2, 3, factor loadings obtained after varimax rotation are presented in Supplementary Tables 2, 3 for the HER2 non-amplified and HER2-amplified groups, respectively.
In the HER2 non-amplified BC cases, Factor 1 was characterized by positive loadings of the variables indicative of CD8 density within the IZ T, TE, and S aspects and was named CD8 density factor. Factor 2 showed positive loadings of HER2 FISH variables (HER2 copy number, HER2/CEP17 ratio, percentage of amplified cells calculated from HER2/CEP17 ratio, and percentage of amplified cells calculated by HER2 signals only) and was named the HER2 amplification factor. Factor 3 was characterized by increasing CD8 densities towards the T aspect of the IZ (strong positive loadings of the CD8 CM and its SD) and by moderate loading of CD8 density in the T Similarly, in HER2-amplified tumors, Factor 1 was the HER2 amplification factor, Factor 2 was the CD8 density factor, and Factor 3 (CD8 density gradient factor) was the main sources of variance. Factor 4 was characterized by strong positive loadings of Ki67% and Ki67 entropy indicators and by moderate negative loading of ER entropy. Factor 5 was represented by the percentage of both HRs along with the PR entropy. Factor 6 was characterized by strong positive loading of a single HER2 MC variable.
Prognostic Significance of Clinicopathological Parameters, Immunohistochemistry, Fluorescence In Situ Hybridization, Immune Response, and Intratumoral Heterogeneity Indicators in HER2 Non-Amplified and HER2-Amplified Patients We explored the potential of the clinicopathological parameters, IHC, FISH, immune response, and ITH indicators for predicting OS of the patients by univariate survival analysis. Statistically significant indicators and their hazard ratios are presented in Table 2. For the HER2 non-amplified group, higher T stage, lymph node status (pN), CD8 density in the S aspect, SD of CD8 density in the S and TE aspects, immunodrop of CD8 density, and Haralick's texture indicators reflecting homogeneity of HER2 and HER2 MC (energy, homogeneity) were associated with shorter OS. Meanwhile, higher HER2 expression, CD8 densities in the tumor area and T aspect within IZ along with its variance, CM for CD8 density and its variance, Haralick's texture indicators reflecting heterogeneity of HER2 and HER2 MC (contrast, dissimilarity, and entropy), and ER contrast were associated with longer OS. In the HER2-amplified patients, worse OS was associated with higher T stage, pN, immunodrop of CD8 density, HER2 MC homogeneity, Ki67 entropy, and PR AshD (bimodality), while in the presence of higher CEP17 copy number, the remaining Immunogradient indicators, HER2 entropy, HER2 MC contrast, and dissimilarity were associated with better OS.
All the variables significantly associated with outcome at a univariate analysis (p < 0.05, Table 2) were assessed for their independent prognostic value in the multivariate Cox regression models.
To investigate any added prognostic value of the indicators, three models in each group were generated from different variable sets (Table 3). Models 1 and 4 were obtained from the pathology and IHC data, including the ITH indicators; FISH  indicators were additionally used in the HER2-amplified group. In models 2 and 5, the IHC CD8 density and Immunogradient indicators were added to the variables tested in models 1 and 4. Models 3 and 6 were obtained from the pathological and CD8 indicators, without inclusion of any ER, PR, Ki67, and HER2 variables.
In the HER2 non-amplified group, higher values of HER2 MC, HER2 MC entropy, and ER contrast indicators were independent features of better OS, while higher tumor stage was associated with worse OS (model 1). Model 2 revealed a marked increase of prognostic power contributed by the immune response indicators in the data set (model likelihood ratio 56.1 achieved in model 2 compared with that of 27.1 in model 1); better OS was associated with higher CD8_CM and CD8_d_T cell densities, and worse OS with higher CD8_d_TE_sd. Remarkably, models 2 and 3 included the same three immune response indicators as independent prognostic factors, reflecting different properties of the local CD8 densities within tumor microenvironment: CD8_d_T (absolute density in the tumor aspect of IZ) and CD8_CM (positive IZ density gradient towards the tumor) were both associated with longer OS, while CD8_d_TE_sd (variance of the CD8 cell density along the IZ) was a feature of worse prognosis.
For the HER2-amplified group, no significant prognostic IHC (ER, PR, HER2, and Ki67 global expression levels) indicators were found by the univariate analyses; therefore, only a set of ITH and FISH indicators along with the pathological variables were used in model 4. Models 5 and 6 were built with the same sets of variables as in the HER2 non-amplified group. In model 4, higher values of HER2 MC contrast and CEP17 copy number indicators predicted better OS, while pN was associated with worse OS. The prognostic power of model 5 was increased by adding immune response indicators (likelihood ratio 29.03 of model 5 compared with 17.64 of model 4), where higher CD8 density in the tumor aspect of IZ predicted better OS. The latter indicator was also an independent factor of better OS in the context of worse OS predicted by pN status in model 6.
The Kaplan-Meier survival probability plots demonstrating an association between the prognostic factors and OS are presented for the HER2 non-amplified and HER2-amplified groups in Figures 1, 2, respectively.
Combined CD8 Immunogradient Prognostic Score in HER2 Non-Amplified Patient Group To further assess the added prognostic value of the independent immune response features revealed by the multivariate regression analysis in the HER2 non-amplified group, a combined CD8 Immunogradient prognostic score was calculated by summing corresponding scores (0/1) for each factor (CD8_CM, CD8_d_T, and CD8_d_TE_sd), assigning the score 1 for good or 0 for poor prognosis. The combined CD8 Immunogradient prognostic score allowed stratification of patients into three prognostic groups with 5-year OS probability of 98%, 80%, and 49% for the score of 3, 2, and 1, respectively ( Figure 3). Of note, there were no patients with all three indicators assigned a score of 0.

DISCUSSION
In this study, we present prognostic models for patients with HER2 IHC borderline (2+) BC patients, based on the expression levels of ER, PR, Ki67, HER2, and CD8 densities in the tumor tissue assessed by DIA. These biomarkers were augmented by a set of computational indicators that quantify spatial aspects of ITH and tumor microenvironment. Importantly, the CD8 (immune response) indicators markedly strengthened the  models in both the HER2 non-amplified and HER2-amplified groups. Furthermore, these latter indicators outperformed pathological variables and enabled independent prognostic stratification of the HER2 non-amplified BC patients. DIA enables extraction of IHC data and spatial aspects from WSI with high capacity, not available by conventional IHC scoring. Additional processing of the DIA-generated data by hexagonal tiling enabled extraction of Haralick's texture measures of ITH of the biomarkers of the prognostic value, as reported previously. In this study, we found that, of the IHC variables explored in the HER2 non-amplified group, only HER2 expression percentage and HER2 MC were significantly  Table 3, models 1 and 2). Similar findings of beneficial prognostic impact of higher HER2 MC were reported recently in early HR-positive BC patients, where better prognosis of higher HER2 expression was found in a univariate analysis (58). In another study, a trend of more favorable prognosis with respect to relapse-free survival has been shown for the ER-positive, HER2 non-amplified tumors with higher levels of HER2 RNA (91). HER2 MC status shows the status of HER2 expression, as HER2 protein is localized on the cell membrane. This means that HER2 MC entropy, which is indicative of MC spatial heterogeneity, also reveals information about the ITH of the HER2 protein expression. We observed a non-linear relationship between the HER2 MC and its spatial heterogeneity in our study (Supplementary Figure 4), represented by high ITH values in    Supplementary Figure 2 and Supplementary Table 2). This HER2&ER heterogeneity factor reflects the higher ITH of both HER2 and ER proteins in the tumors with decreased ER expression. A majority of the patients with HER2-amplified tumors received adjuvant trastuzumab treatment (87,74.4%). The OS of these patients is likely to have been impacted by the targeted therapy; therefore, the prognostic models obtained in this subgroup should be taken with caution. One can speculate that any potential effect of the targeted therapy could be related to our finding of ITH of HER2 expression, represented by the HER2 MC contrast indicator as independent predictor of better OS ( Table 3, HR = 0.35, p = 0.0367, Model 4 and HR = 0.243, p = 0.0077, Model 5) but not by the HER2 MC indicator. The effect on better OS caused by a higher CEP17 copy number in this group is not clear, and it may be related to various treatment modalities applied in HER2-amplified BC patients. Several studies have reported an association between CEP17 copy number gain and responsiveness to anthracycline-based chemotherapy (92)(93)(94). Also, in addition to HER2, chromosome 17 includes other genes involved in BC pathogenesis and DNA repair, such as BRCA1, TOP2A, TP53, and RAD51C (95, 96); therefore, various abnormalities of chromosome 17 may affect prognosis and treatment response.
In this study, we tested the prognostic value of CD8 cell densities quantified by DIA in the tumor and stroma compartments and applied a recently proposed method, based on hexagonal grid analytics of the DIA data to compute CD8 local density profiles (Immunogradient) across automatically detected tumor-stroma IZ (86). This method actually tests if the immune cells reveal increasing densities towards the tumor at the tumor/host interface and therefore is expected to be more sensitive to capture "spatial behavior" of TILs. Indeed, higher CD8 cell densities in the tumor compartment were associated with better OS in univariate analyses in both patient subgroups ( Table 2, HR = 0.37, p = 0.017 and HR = 0.38, p = 0.024 in the HER2 non-amplified and HER2-amplified groups, respectively); however, they did not provide the independent prognostic value in our models. In contrast, three Immunogradient indicators provided the independent prognostic value in the non-amplified tumors: CD8 density in the tumor aspect of IZ (CD8_d_T) and positive IZ CD8 density gradient towards the tumor (CD8_CM) were associated with better OS, while the variance (SD) of CD8 density (CD8_d_TE_sd) along the TE predicted worse OS. A strong prognostic stratification was achieved by aggregating these three independent spatial properties of the CD8 cell distribution in the tumor microenvironment into a combined CD8 Immunogradient prognostic score; this represents an instance of computational augmentation of a single IHC biomarker ( Figure 3). Remarkably, these three indicators were sufficient to predict OS independently of any other variables ( Table 3, model 3) with statistical power obtained from pathology and IHC data supplemented with ITH indicators ( Table 3, model 1). Finally, the prognostic power was doubled by adding the immune response indicators to the model ( Table 3, model 2). Our findings are similar to the results presented in the study of Rasmusson et al. (86), where both CD8 density in the tumor aspect of IZ and CM for CD8 cell density within the IZ indicators were independent predictors of better OS in early HR-positive BC. Although several studies have reported a higher density of CD8 cells to be associated with a favorable prognosis in node-negative BC (97), or in combination with CD163 (98), other studies have shown an adverse prognostic effect of increased CD8 lymphocytes in patients with HR-positive/HER2-negative tumors (77,99,100) or reported no significant association between CD8 cells and patient outcome (79). These contradictory results in HRpositive BC may be related to different methodologies applied, lacking precision in the assessment of spatial aspects of TIL distributions within the tissues (100,101). Recently, Dieci et al. (102) highlighted the need of deeper insight into the mechanisms on which the interaction between HR-positive/HER2-negative BC tumor and immune cells relies, as various factors such as menopausal status, estrogen levels, and endocrine treatments may be involved in the modulation of the tumor microenvironment (102). Therefore, methods with appropriate discriminatory spatial precision are needed to expose the prognostic role of TILs in luminal-like BC.
We did not find significant differences in CD8 cell densities between the HER2-amplified and non-amplified groups, which could be explained by the fact that the HER2-amplified group was composed of both molecular subgroups showing HER2 positivity, namely, luminal B and HER2-enriched. Previous studies reported that HER2-enriched subtype is more immunogenic than the luminal B (103). In our study, the only immune response indicator-density of CD8 in the T aspect of IZ-provided an independent association with better OS in HER2-amplified BC patients ( Table 3, models 5 and 6). Extensive TIL infiltration has been associated with better outcomes (pathological complete response, event-free survival, and disease-free survival) in HER2-positive BC (70,104,105). However, studies evaluating the prognostic significance of CD8 TILs reported conflicting results (77,79,100,106,107), suggesting that the association between CD8 cells and prognosis depends on lymphocyte types, their tissue location, analysis methods, and other factors. Indeed, the interaction between immune system and tumor as well as prognostic effects of TILs in HER2-positive BC is impacted by various combined therapy modalities, including anti-HER2 therapy, chemotherapy, and hormonal therapy (108). Trastuzumab therapy effect depends on immune response (109), it has both cytotoxic and immunological effects on tumor cells (110)(111)(112), and better therapeutic efficacy is achieved in tumors with high TILs (113)(114)(115). However, this was not confirmed by other studies (116,117).
Our study has limitations, related to its retrospective and monocentric design and lack of well-structured information about applied therapies and responses achieved. In particular, it is relevant to the prognostic modeling in the HER2-positive patient cohort.
In conclusion, we present prognostic OS models based on computational ITH, tumor microenvironment, and immune response indicators of the IHC biomarkers in HER2 IHC 2+ borderline BC patients. The ITH indicators (HER2 MC entropy and ER contrast in FISH-negative and HER2 MC contrast in FISH-positive tumors) provided an independent contribution to predict better OS. In FISH-negative tumors, antitumor immune response, assessed by the CD8 IZ Immunogradient indicators, provided prognostic stratification independent and superior to other pathology and IHC variables.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

ETHICS STATEMENT
The studies involving human participants were reviewed and approved by Lithuanian Bioethics Committee. The patients/ participants provided their written informed consent to participate in this study.

AUTHOR CONTRIBUTIONS
GR, ArL, and AR designed the study. GR in collaboration with AiL collected the samples and participated in IHC slide staining. GR, RG, and RA performed digital image analyses. GR and RA performed statistical analyses. DZ participated in the statistical analysis of CD8. GR in collaboration with AR and ArL drafted essential parts of the manuscript. All authors reviewed the analysis results and read and approved the final manuscript.