Multidimensional profiling of drug‐treated cells by Imaging Mass Cytometry

In pharmaceutical research, high‐content screening is an integral part of lead candidate development. Measuring drug response in vitro by examining over 40 parameters, including biomarkers, signaling molecules, cell morphological changes, proliferation indices, and toxicity in a single sample, could significantly enhance discovery of new therapeutics. As a proof of concept, we present here a workflow for multidimensional Imaging Mass Cytometry™ (IMC™) and data processing with open source computational tools. CellProfiler was used to identify single cells through establishing cellular boundaries, followed by histoCAT™ (histology topography cytometry analysis toolbox) for extracting single‐cell quantitative information visualized as t‐SNE plots and heatmaps. Human breast cancer‐derived cell lines SKBR3, HCC1143, and MCF‐7 were screened for expression of cellular markers to generate digital images with a resolution comparable to conventional fluorescence microscopy. Predicted pharmacodynamic effects were measured in MCF‐7 cells dosed with three target‐specific compounds: growth stimulatory EGF, microtubule depolymerization agent nocodazole, and genotoxic chemotherapeutic drug etoposide. We show strong pairwise correlation between nuclear markers pHistone3S28, Ki‐67, and p4E‐BP1T37/T46 in classified mitotic cells and anticorrelation with cell surface markers. Our study demonstrates that IMC data expand the number of measured parameters in single cells and brings higher‐dimension analysis to the field of cell‐based screening in early lead compound discovery.

In pharmaceutical research, high-content screening is an integral part of lead candidate development. Measuring drug response in vitro by examining over 40 parameters, including biomarkers, signaling molecules, cell morphological changes, proliferation indices, and toxicity in a single sample, could significantly enhance discovery of new therapeutics. As a proof of concept, we present here a workflow for multidimensional Imaging Mass Cytometry TM (IMC TM ) and data processing with open source computational tools. CellProfiler was used to identify single cells through establishing cellular boundaries, followed by histoCAT TM (histology topography cytometry analysis toolbox) for extracting single-cell quantitative information visualized as t-SNE plots and heatmaps. Human breast cancer-derived cell lines SKBR3, HCC1143, and MCF-7 were screened for expression of cellular markers to generate digital images with a resolution comparable to conventional fluorescence microscopy. Predicted pharmacodynamic effects were measured in MCF-7 cells dosed with three target-specific compounds: growth stimulatory EGF, microtubule depolymerization agent nocodazole, and genotoxic chemotherapeutic drug etoposide. We show strong pairwise correlation between nuclear markers pHistone3 S28 , Ki-67, and p4E-BP1 T37/T46 in classified mitotic cells and anticorrelation with cell surface markers. Our study demonstrates that IMC data expand the number of measured parameters in single cells and brings higher-dimension analysis to the field of cellbased screening in early lead compound discovery.
Discovery of new treatments in oncology research relies extensively on the use of human-derived cell culture models [1]. High-content cell-based screens are widely applied in pharmaceutical drug development to prioritize lead molecules for animal testing [2]. These assays rely on the use of primary and cancer cell lines and mostly monitor cytotoxicity and proliferation. With the advent of sophisticated genomics and proteomics technologies for soluble proteins and mass cytometry for single cells, it is becoming possible to increase the multidimensionality of in vitro screens. Databases of genotypic and phenotypic profiles [3][4][5] across cancer drug panels have provided a comparative analysis between in vitro studies and clinical therapeutic responses in vivo [6]. One of these datasets developed and shared by the National Cancer Institute 60 (NCI-60) in the late 1980s was the first in vitro discovery screening tool for pharmacologic compounds inducing growth inhibition across 60 cancer cell lines [7]. Over the years, several concerns were raised regarding the use of in vitro established human tumorderived cell lines for drug testing. Concerns included genetic instability and dedifferentiation. However, the study of genetic mutations arising in immortalized cells remained a field of interest in drug discovery, providing knowledge on the dysregulation of cellular signaling pathways and the effects small-molecule inhibitors have on human tumor cell lines [8].
High-content imaging of cancer cell lines in response to drug treatment is a standard assay applied in preclinical studies for identification of different mechanisms of drug action [9][10][11][12]. The main criterion in every morphophenotypic screen for early drug discovery is the selection of biomarkers and detection modalities (antibodies, chemical and enzymatic probes, reporter detection tags). In the field of fluorescent cellular imaging microscopy, biomarker analysis in single cells is limited to 4-6 due to the overlapping spectra of fluorescent dyes [13]. Repetitive rounds of staining using the same biological sample are required to achieve multiplexing [14]. New imaging technologies are needed to significantly increase multiplicity in a single sample preparation experiment and perform replicate analyses [15]. High-content imaging assays applied to screening drug perturbations in heterogeneous cancer cell models could uncover additional modes of drug action and biomarkers for clinical trials [16].
Imaging Mass Cytometry is an emerging and transformative technique in the field of digital histopathology applied to complex tissue sections [17][18][19]. The Hyperion TM Imaging System (Fluidigm Ò , Toronto, Canada) can measure up to 40 parameters simultaneously in formalin-fixed, paraffin-embedded (FFPE), and frozen human tissue sections with subcellular resolution [20]. Sample preparation is very similar to standard immunohistochemical protocols [17], where the tissue section on a slide is first deparaffinized and treated with an antigen retrieval buffer for antigen epitope exposure in the case of FFPE, then stained once with a mixture of isotopically pure metal-labeled antibodies specific to structural and cell type biomarkers. After drying and inserting the slide into the ablation chamber, regions of interest (ROIs) are chosen and recorded by a camera integrated with the Hyperion Tissue Imager (Fluidigm Ò ). These ROIs are then ablated in 1-lm steps as the slide moves under the laser. Each laser pulse vaporizes a 1-lm 2 area of the tissue and generates a plume of particles. The resulting particles are carried by a stream of helium/argon mixture into a CyTOF Ò instrument, an inductively coupled plasma time-of-flight mass cytometer, described elsewhere [18]. Data are stored as raw binary files of ion counts for each mass channel (parameter) per pixel. These files can be converted into a stack of single-parameter matrices (grayscale images), which are used to create merged multiparametric pseudocolor images and perform quantitative pixel-or segmentation-based analysis. Bodenmiller et al. are developing multivariate computational tools to visualize and analyze multiplexed images of human tissue sections generated by IMC [21]. However, there is no single software package or analysis workflow that could currently be applied to answer specific biological questions.
In this proof-of-principle study, we set out to develop a comprehensive workflow for IMC data analysis based on recent advances in imaging algorithmic methods to visualize and measure multiple biomarkers in model cell lines cultured in chamber slides (Fig. 1). Several human breast cancer-derived cell lines, SKBR3, HCC1143, and MCF-7, were screened for expression of surface and intracellular markers. Predicted pharmacodynamic effects were studied in MCF-7 cells dosed with three target-specific compounds: growth stimulatory epidermal growth factor (EGF), microtubule depolymerization agent nocodazole, and genotoxic chemotherapeutic drug etoposide [22,23]. Our method analysis workflow demonstrates the high multiplexing capability of IMC for in vitro research and offers quantitation of multiple biomarkers at subcellular resolution. Future improvements of the technology toward a higher acquisition speed of multiple samples will expand IMC applications as an imaging platform for in vitro cell-based drug profiling. ) were used to seed cells at a concentration of 0.1e6ÁmL À1 in 0.5 mL per chamber and grown in 37°C, 5% CO 2 , 100% humidity for 48 h prior to drug compound addition. Each breast cancer cell line was seeded using individual chamber slides to avoid any cell cross-contamination, and each IMC experiment was carried out independently. In chronological order, seeding of HCC1143 was the first sample preparation, then we run the second experiment using SKBR3, and last we used MCF-7 for the drug-treated cells in vitro study.

In vitro drug treatment
All chemical compounds were dissolved and aliquoted at the supplier's recommended concentrations for long-term storage in 100% DMSO (Sigma-Aldrich Cat. No. 276855). Final concentrations for each compound were selected from previously published data [24][25][26] by premixing the required initial drug stock volume with full growth media and then transferring 0.5 mL volume to corresponding cell chambers (Table S1): 10 lM of etoposide (Cell Signaling Technology Ò , Danvers, MA, USA; Cat. No. 2200S); 0.5 lM of nocodazole (Sigma-Aldrich Cat. No. M1404); 10 ngÁmL À1 of human EGF (Sigma-Aldrich Cat. No. E9644); 2% DMSO in growth media. One chamber was left with fresh growth media as nontreated control. Cells were exposed to compounds for 48 h at 37°C, 5% CO 2 , and 100% humidity. Preparation of additional eight chamber sample slides to replicate our experiment for IMC analysis was not required due to the known in vitro predicted effect of the drugs we tested on this cell line.

Immunocytochemistry and antibody validation
Immunostaining of live cells was performed in chamber slides. Following drug treatment, fresh media in all chambers was supplemented with Cell-ID TM Intercalator-103Rh (Fluidigm Cat. No. 201103A) at 1 : 500 for 15 min at 37°C, 5% CO 2 for dead cell identification. The chamber slides were washed with DPBS (Thermo Fisher Scientific Cat. No. 14190-144) at room temperature (RT) for 5 min. Surface metal-labeled antibody cocktail (Table S2) was made in 0.5% BSA/DPBS from antibody stocks at a dilution of 1 : 100 for each. Cell monolayers were stained with 250 lL antibody mix for 120 min, RT. Following a washing step, cells were fixed with 250 lL freshly made 1.6% formaldehyde (Thermo Fisher Scientific Cat. No. 28906) in DPBS for 15 min at RT. After fixative removal, cells were permeabilized with fresh 0.3% Saponin/DPBS (Sigma-Aldrich Cat. No. S7900-25G) for 30 min at RT and then blocked with 1% BSA/DPBS for 1 h at 37°C. After blocking, cells were incubated overnight at 4°C in 250 lL freshly prepared cocktail of metal-labeled antibodies targeting intracellular markers (Table S3) diluted 1 : 50 in 0.3% Saponin/DBPS. The following day cells were washed once with 0.5% BSA and stained with 250 lL of nuclear Cell-ID Intercalator-Ir for 30 min (1 : 1000), RT. After the last wash, each chamber was filled with 250 lL deionized water for 5 min for salt removal. Chambers were detached from slides and left to dry at RT until IMC analysis. The commercial antibodies of our panels have been previously validated by immunocytochemistry applied to immunofluorescence analysis of adherent cells. A similar sample preparation method was applied to the validation of these antibodies labeled with metals regarding their antigen specificity detected by IMC. The selection of titers tested for surface markers (1 : 100) and intracellular markers (1 : 50) refers to rigorous quality control antibody validation functional assays using single-cell suspension Mass Cytometry, which has the same single isotopic detection modalities of antigens as IMC. The list of antibodies used for each breast cancer cell line experiment is cross-referenced in Tables S2 and S3. A total number of 14 antibodies were tested on HCC1143 cells, 20 markers for SKBR3 and 25 for MCF-7. The markers from each panel were selected based on their specificity to surface membrane, cytoplasmic and nuclear components of single cells, their previous applications in the context of breast cancer cell line studies, and their commercial availability. Expanding a panel to 37 markers requires additional custom labeling of target-specific monoclonal antibodies with nonused metal isotopes.

Data acquisition by IMC
Samples were analyzed with the Hyperion Imaging System (Fluidigm). The dried slide was loaded into the imaging module, where an optical previewing of the ROIs was recorded for laser ablation. Areas of dimension 1000 9 1000 lm were acquired for HCC1143 and 1400 9 1400 lm for SKBR3 cell lines experiments. Replicate ROIs of 1500 9 1500 lm size were collected for MCF-7 exposed to chemical compounds, and a single ROI was ablated for nontreated control and DMSO-treated cells. Each ROI sample acquisition took 4.5 h at an ablation frequency of 200 Hz. The resulting data files were stored in MCD binary format. Multicolor images were generated with open source IMAGEJ 1.51 software [27]. Zoom-in regions for each multicolor image were made by selecting and cropping areas of dimension 400 9 400 lm from the original picture without any change of the pixel size. The resolution of each IMC image shown is 1 lm for 1 pixel.

Image analysis software tools
For each recorded ROI, stacks of 16-bit and 32-bit singlechannel TIFF files were exported from MCD binary files using MCD TM Viewer 1.0 (Fluidigm). Cell-based morphological segmentation was carried out with two image processing pipelines for 16-bit or 32-bit TIFF files and CELLPROFILER 2.2.0 (CP, Broad Institute, Cambridge, MA, USA), a widely adopted software in the open source image analysis community which has been continually improved since its availability in 2005 to read and analyze cell images using advanced algorithms [28]. The 16-bit TIFF unstacked image format was used as input data in CELLPROFILER to run a segmentation pipeline with a set of sequential modules to generate and save an unsigned 16-bit integer singlecell mask TIFF image. Inputs of 16-bit TIFF images with their corresponding segmentation mask were uploaded in histoCAT to open a session data analysis. Due to the high number of 29 single isotopic channel TIFF images for each replicate and drug treatment condition, histoCAT analysis was generated faster with a 16-bit format instead of a 32bit image format which is bigger in size. The 32-bit TIFF unstacked images were loaded in CELLPROFILER to run a segmentation and mean intensity multiparametric measurement analysis of individual cells. The resulting outputs were exported and saved as SQLite single-cell object measurements database files directly compatible with CELLPRO-FILER ANALYST for downstream computational analysis [29]. Spatial distribution maps, dimensionality reduction, and unsupervised clustering for 16-bit single images were performed using the histoCAT 1.73 open source toolbox [21]. Supervised analysis by support vector machine-learning classification on 32-bit datasets was processed with CELLPRO-FILER ANALYST 2.2.1 [30]. Hierarchical clustering and correlation heat maps of classified cell populations were prepared using the web-based tool Morpheus (GenePattern, Broad Institute) [31]. Clustering and visualization of different types of markers in classified populations of cells across all controls and drug treatments were performed with the free software environment R 3.5.3, the packages pheatmap 1.0.12, igraph 1.2.4 and edgebundleR 0.1.4 [32].

Statistical analysis
Data from nontreated control, DMSO-only, and drug-treated cells were analyzed using a nonparametric Mann-Whitney test with a two-tailed P value for every channel independently. Statistical significance was defined as P < 0.05. All statistical analysis was performed with GRAPH-PAD TM PRISM Ò 7.04 software [33].

Multiparametric characterization of breast cancer cell lines by IMC
Imaging Mass Cytometry analysis revealed heterogeneity in the expression of surface, intracellular, and nuclear markers in human breast cancer cell lines cultured in chamber slides. As an example, we selected the SKBR3 cell line, which is characterized by an invasive phenotype and increased proliferation in vitro and is used as a model to delineate mechanisms of resistance to ErbB2-targeted clinical therapies [34,35]. As seen in Fig. 2, the culture consists of mostly smalland medium-size rounded epithelial cells and occasional large multinucleated cells. Most cells show high levels of the cytoskeletal marker pan-keratin. Human epidermal growth factor receptor 2 (HER2) and EGFR are seen at the surface membrane, while tumor suppressor and transcription factor p53 are localized in the nuclei. Presence of lysosomal organelles is outlined in the cytoplasm with an antibody targeting the lysosome-associated membrane protein 2 (LAMP-2), known as CD107b. Image resolution is similar to fluorescence microscopy. Ductal breast carcinoma-derived triple-negative HCC1143 cell line was used to demonstrate mesenchymal morphology ( Fig. 3) with high vimentin detection in the cytoplasm, nuclear localization of p53, and cytoplasmic distribution of the basallike breast cancer marker cytokeratin 5 [36,37]. The proliferation marker Ki-67 identifies cells in the active state of the cell cycle. Combination of multiple markers in an in vitro cell profiling assay allows more precise analysis of cellular states of differentiation and invasiveness characteristic for breast cancer tumors.

In vitro drug effect profiling by IMC
To evaluate how IMC may be translated into future applications in the field of cell-based drug profiling, we initiated a limited drug treatment study using one model breast cancer cell line. The phenotypic and functional responses of wild-type p53-expressing MCF-7 cells to etoposide, nocodazole, or EGF compounds were investigated. In preclinical in vitro drug discovery, this model is considered to be reproducible, fast, and inexpensive despite the lack of clinical correlation to drug response in vivo due to breast tumor heterogeneity [38]. To identify how these chemical compounds induce cytostatic and cytotoxic effects on this noninvasive luminal breast cancer cell line, we generated multiplexed pseudocolor images of affected cellular compartments (Fig. 4). For example, nuclear size and proliferative activity of each cell were analyzed using the nuclear intercalator-Ir stain (Cell-ID Intercalator-Ir), cell cycle-S phase marker (Cell-ID IdU), proliferation marker Ki-67, cell cycle regulatory proteins cyclin B1 and D3, and estrogen receptor alpha (ERa). Antibody against pH2A.X S139 was used as a nuclear marker for DNA double-strand breaks and pHistone3 S28 for mitotic cells. Cytoskeletal markers pan-keratin and cytokeratin 19 were used to follow changes in morphology and cell size. Drug response of MCF-7 morphology, adhesion, and cell-to-cell interaction were visualized with surface membrane proteins, such as the membrane-bound mucin marker MUC-1 overexpressed in breast carcinoma [39], cell adhesion molecule EpCAM, integrin CD29, tetraspanin CD81, the low-expressed epidermal growth factor receptor EGFR, regulators of integrin-mediated cell adhesion CD98 and CD47 [40], and E-cadherin as epithelial interaction protein (Table S2). Long-term exposure of MCF-7 cells to EGF induced significant morphological changes and an increase in cell size. Lineage shift from a classic epithelial cobblestone monolayer to a more invasive mesenchymal type was observed. The presence of large spindle-shaped cells that seemed to lose contact with other cells and migrate into empty spaces may be related to an increase in cell motility (Fig. 4A). Activation of the EGFR signal transduction pathway by EGF is known to play a role in cellular motility and size increase in the tumor aggressiveness [41].
Screening of MCF-7 cells treated with etoposide, nocodazole, or EGF shows multivariate phenotypic response profiles compared to nontreated and DMSOtreated controls (Fig. S1).
Treatment of MCF-7 cells with the antineoplastic agent nocodazole resulted in an increase in the population of pHistone3 S28 -positive cells with a mitotic index threefold higher compared to DMSO-treated or nontreated control (Table S4). Another cytostatic phenotype was observed after prolonged drug exposure: subpopulations of cyclin D3-positive and tetraploid cells exiting mitosis without cytokinesis as a result of microtubule disruption by nocodazole (Fig. 4B). Presence of cyclin D3-positive and tetraploid cells in the G1/S interphase confirms nocodazole effect on MCF-7 cell cycle [24].
Incubation of MCF-7 with the anticancer agent etoposide shows inhibition of cell proliferation, elevated numbers of cells positive for nuclear pH2A.X S139 protein, and compaction of the monolayer compared to EGF-or nocodazole-treated samples (Fig. 4C). Exposure of tumor-derived cell lines to this inhibitor induces the formation of a stable covalent topoisomerase II-cleaved DNA complex, which causes multiple breaks in double-stranded DNA. Etoposide treatment of MCF-7 activates phosphorylation of H2A.X, which is required for checkpoint-mediated cell cycle arrest, promotes DNA repair, and maintains genomic stability. Thus, pH2A.X S139 is a sensitive biomarker for measuring levels of drug-induced double-stranded DNA breaks in cancer cells [42].

High-dimensional phenotype analysis of in vitro drug response
Multivariate phenotypic responses of MCF-7 cells following chemical perturbations were quantified with a set of available software tools. Unfortunately, there is no single software package for this type of analysis. Therefore, we employed a workflow that includes a wide variety of tools designed for cellular imaging. Hopefully, this work will underscore the need for software specifically developed for analysis of multiparametric data acquired by IMC. First, we used CELLPROFILER 2.2.0 (CP) to identify single cells through establishing nuclear and cellular boundaries [43]. Accuracy and robustness of cell segmentation were assessed by using Cell-ID Intercalator-Ir to identify cell nuclei and pan-keratin for cytoplasm (Fig. 5A, Figs S2A and S3). The segmentation masks were generated for each ROI and then exported for further downstream analysis of changes in MCF-7 cellular phenotypes in response to drug treatment.
Second, to visualize and analyze data we used the open source platform histoCAT, which enables highly multiplexed, quantitative, and detailed analysis of cell phenotypes, microenvironment interactions, and tissue architecture [21]. The workflow of histoCAT consists of overlaying the segmentation masks previously generated by CP to extract single-cell level information.
To evaluate the diversity of phenotypic profiles collected from drug-treated MCF-7 cells, we then selected and generated Z-score normalized mean intensities of pHistone3 S28 and pH2AX S139 as drug-sensitive protein markers. Spatial heat maps of both channels were displayed for various drug-treated MCF-7 ROIs using a 99th percentile cutoff (Fig. 5B,C, Fig. S2B). This quantitative and visual approach allows us to compare at the same scale the spatial distribution of these parameters between ROIs. For example, pH2A.X S139 nuclear expression level in etoposide-treated MCF-7 cells is higher than in nocodazole or EGF-treated cells on average. Mitotic cell marker pHistone3 S28 expression level in nocodazole-treated MCF-7 is twofold higher than that identified for other drugs. Related single-cell Z-score mean intensity values for nuclear markers (pHistone3 S28 , pH2A.X S139 , Ki-67, p4EBP1 T37/T46 , p53, cyclin D3), surface markers (CD98, CD81, CD29, CD49e, CD47), and cell size parameters (area, perimeter, major and minor axis lengths) were exported from each drug treatment dataset to assess statistically significant changes of protein expression levels and cellular morphology ( Fig. 5D and Fig. S4). We can see differences between both controls and drug-treated samples in the expression of nuclear markers, with a coincrease of pH2A.X S139 and Ki-67 levels, and phosphorylation of the translation repression protein 4E-BP1 under etoposide treatment. This suggests the presence of a nonquiescent population of MCF-7 undergoing DNA repair and maintaining high proliferative activity. For nocodazole and EGF treatments, low Ki-67 levels may be related to the disruption of cell cycle and lineage shift induced in response to respective compound [44]. EGF treatment shows increased cell size and expression of adhesion surface markers compared to other conditions. These quantitative results are in agreement with multicolor images of treated cells. Finally, the heterogeneity of MCF-7 cell phenotypes was analyzed using the t-distributed stochastic neighbor embedding (t-SNE) algorithm, a data dimensionality reduction method based on the similarity of selected markers or cell types [45]. We processed replicates of ROI for each compound using Z-scorenormalized protein markers and cell size parameters to generate two-dimensional t-SNE plots. This visualization tool shows signal distribution of phosphoproteins, size parameters, and nuclear and cell surface membrane markers over the different ROIs ( Fig. 6A and Fig. S5). For example, the population of mitotic cells was grouped in a single region on the map with high expression levels of p4EBP1 T37/T46 , pHistone3 S28 , and Ki-67. The high level of pH2A.X S139 was detected in the region of etoposide-treated cells and correlates with the strong Ir-intercalator signal due to compaction of cellular DNA in response to the drug. Area plot reveals the presence of larger size cells in EGF and nocodazole treatments, in line with EGF-induced lineage shift and cytokinesis inhibition by nocodazole (Fig. 6B).

Classification of drug-treated MCF-7 phenotypic changes by support vector machine learning
Here, we applied the support vector machine-learning classification ruler termed Fast Gentle Boosting for accurate IMC dataset multiparametric analysis of nuclear states of compound-treated MCF-7 cells [29]. This classifier was integrated into CELLPROFILER ANALYST (CPA) and trained [30] to identify five different nuclear classes of MCF-7 cells: pHistone3 S28 -positive mitotic cells, pH2A.X S139 for cells undergoing DNA repair, p53 as tumor suppressor activation, cyclin D3 for G1/S transition in interphase, and Ki-67 for proliferating cells. A composite image of all nuclear markers is shown in Fig. 7A. These single-cellular objects were then scored and classified for a supervised analysis of every protein marker and cell size parameter per specified nuclear subpopulation. Accuracy of image classification on each dataset of MCF-7 drug compound was evaluated by cross-validation between the true labeled cells as classification training set and the image full scoring as predicted label set for each nuclear class. Classification reports were displayed as confusion matrices where each value corresponds to average cross-validation metrics between the trained (true label) and scored (predicted label) single-cell objects for each nuclear class on a performance scale from 0 to 1 (Fig. 7B and Fig. S6). As shown, classification of each nuclear class per drug compound dataset shows relatively good performance between 0.5 and 1.
Raw mean intensity values of protein expression for each of the five nuclear markers were exported per classified MCF-7 population and compared to evaluate statistical significance (Fig. 7C and Fig. S7). Bar charts show various levels of protein markers among the five predicted classes of nuclear states, with higher protein detection in the corresponding class, except for the Ki-67 proliferative marker, which shows stronger detection in pHistone3 S28 mitotic cells.
Multiparametric MCF-7 drug compound data prelabeled by nuclear classification were reprocessed using hierarchical clustering similarity heat maps. First, we computed pairwise similarities for each class of cells per drug exposure with Pearson correlation coefficient between each measured parameter. This method is applicable to non-normalized data and can be used when data measurements vary between samples. Similarity Pearson correlation heat maps were computed using the browser-based tool Morpheus [31]. The core interface in Morpheus is a graphical color heat map representing protein marker and size parameter measurements applied to each single-cell object. Determination of the Pearson correlation coefficient between each measured parameter generates a reorganized heat map of multiparametric similarity relationships for each classified MCF-7 drug dataset [46]. Using these hierarchical similarity clustering heat maps, we compared the grouping pattern of mitotic pHistone3 S28positive cells across the multiple ROIs of compoundtreated MCF-7 ( Fig. 8 and Fig. S8). Naturally occurring isotopes of iridium, 191 Ir and 193 Ir, are present in the nuclear staining reagent Cell-ID Intercalator-Ir and are grouped consistently across all heat maps. The two cytoskeleton biomarkers pan-keratin and cytokeratin 19 strongly correlate as well, since anti-pan-keratin and anti-CK19 antibodies identify cytokeratin isoforms highly expressed in MCF-7. Size parameters including perimeter, area, and axis length of each mitotic cell are grouped together after hierarchical clustering across all drug treatments.
Cell surface protein markers (CD29, CD98, CD81, CD47, EpCAM, and CD49e) generated a cluster with high correlation and redundancy in the pHistone3 S28 mitotic cell subclass. However, EGF treatment influenced expression of these markers, and some level of dissimilarity occurred between the surface membrane cluster and two integrin-mediated regulators, CD98 and CD47. This is consistent with the fact that MCF-7 cells switch to a more mesenchymal phenotype after EGF treatment and higher motility, with a concomitant increase in the expression level of some adhesion proteins.
The intracellular and nuclear biomarkers pHis-tone3 S28 , p4EBP1 T37/T46 , and Ki-67 for mitotic cells  cluster together in all conditions, showing a specific phenotypic profile. By looking at the correlation between the mitotic marker pHistone3 S28 and other markers, we see that p53, cyclin B1, and cleaved poly (ADP-ribose) (cPARP) show variations in their correlation similarities depending on the drug treatment. p53 and cPARP exhibit close relationships with other nuclear markers, since both bind to nuclear DNA to be functional. Cyclin B1 shows strong correlation with pHistone3 S28 in cells from nontreated control and EGF-treated ROIs. This correlation is less pronounced in the case of etoposide and nocodazole. The clustering of different types of markers in the class of pHistone3 S28 mitotic cells was measured by comparing the similarity level of each target-specific cluster across all controls and drug treatments. First, we reprocessed cell surface membrane CD markers for each treatment as a new hierarchical similarity heatmap by using Pearson correlation as distance metric and complete linkage method (Fig. 9A). As expected, each individual cluster is aligned and identified around the matrix diagonal component. Both nontreated and DMSO controls group of protein markers show close similarities to each other, and low or negative correlation with drug-treated markers. We used a hierarchical visualization edge bundle circle plot diagram to highlight the positive intra-and intersimilarities across surface membrane protein parameters (Fig. 9B). Visualization of the high pairwise correlations within each cluster is shown as edge connecting nodes, and bundled lines correspond to close inter-relationships between both control conditions. Drug-treated clusters of proteins show poor or no interdependency. A similar visual workflow was applied to the multivariate analysis by pairwise Pearson correlation of specific nuclear markers and cytoplasmic markers types (Fig. S9).

Discussion
The goal of this study was to demonstrate a multidimensional IMC workflow and quantitative data processing using open source computational tools for in vitro drug response research. The multiplexed measurements acquired from a single image may provide a more informative and reliable screening of drug leads for further clinical development [47]. In this work, phenotypic features of tumor-derived cell lines such as SKBR3, HCC1143, and MCF-7, which recapitulate different subtypes found in breast tumors, were characterized. The wild-type p53-expressing MCF-7 epithelial cell line, a model of luminal breast adenocarcinoma, was selected for in vitro drug treatment and exposed for 48 h to three classes of bioactive compounds: EGF, nocodazole, and etoposide. The cells were stained with a mixture of metal-tagged antibodies against cell membrane, cytoplasmic, and nuclei markers and subjected to IMC. The panel of metal tag antibodies used in our study was designed to minimize any potential spill over cross-interference between neighbor isotopic mass channels [48,49]. When analyzing our IMC data, we did not encounter strong and significant channel contamination from one metal-labeled antibody to another. The collected data allowed us to develop a workflow for cellular imaging analysis, which relies on established computational methods [50]. Multiplexed images of cell compartments give a detailed and broad visualization of biomarker localization and cell size parameter changes. Cell-based morphological segmentation of generated digital images is a prerequisite to identify and export multiple features of interest at single-cell resolution. Exploratory downstream analysis was performed with the histoCAT platform to generate, normalize, measure, and compare multiple images according to their multidimensional content in the form of spatial distribution heat maps, dimensionality reduction, and unsupervised clustering. In this report, we also describe cell classification by machine learning to prelabel cells before measuring their multivariate similarities. Prior knowledge of different cellular states is important for the calculation of distance and similarity metrics between each marker in the evaluation of the heterogeneity of cells and the effects of different compounds on cancer cells [22]. In the example of pHistone3 S28 -positive mitotic cells, we generated a multivariate high-dimensional heat map that allows a comparative study of mechanisms of action of multiple drugs. Additional nuclear classes such as pH2A.X S139 , p53, cyclin D3, and Ki-67 were analyzed using similar multiparametric clustering approaches (data not shown). The strategy could enable new ways to perform deep proteomic drug screening on different models for translational medicine and preclinical trials [51,52].
Other in vitro assays using stem cells, or cocultures of different cell types, could be studied by IMC as an additional approach to designing high-content analysis to investigate cell-to-cell interactions and signaling. A recent targeted therapy against cancer oncogenes described a cytotoxic small-molecule compound that activated the steroid receptor signaling pathways and led to cell stress and death [53]. IMC may assist researchers to better understand perturbations in this and other pathways through high-dimensional analysis. Three-dimensional tumor spheroids grown in vitro are more representative of the complexity of cancer tissue than two-dimensional cell cultures and can be analyzed by IMC through the preparation of sections of FFPE pelleted spheroids [54]. We summarized in Table 1 several advantages and challenges during preparation on slide, staining procedure, acquisition by IMC, and software data analysis for three different models of sample. IMC was useful in the biodistribution study of the platinum-containing anticancer drug cisplatin in pancreatic adenocarcinoma xenograft model and could be applied to the emerging metalcontaining chemotherapies and photodynamic therapies, which are considered as alternative approaches to classic compounds [55,56]. Therefore, the goal for IMC in preclinical profiling of metallodrugs will be to efficiently and effectively identify cellular responses of in vitro tumor-derived cell lines for lead optimization in target-based drug discovery [57]. A comparison of IMC performance metrics to other high-content imaging system technologies in the field of single-cell analysis is shown in Table 2.
In conclusion, we demonstrate a comprehensive image analysis workflow for cell-based IMC data analysis of surface and intracellular markers in drug-treated cells. The correlation coefficient for pairs of multiple parameters and the similar distances among a collection of treatment profiles facilitate downstream analysis and allow for direct data visualization. Image-based cell profiling studies compute statistical estimates of the  likelihood of equivalence between two drug-induced profiles. Similarity measurements quantify proximity between profiles, since they detect deviations from one sample to another regardless of the absolute magnitude. This procedure is useful in finding relations and groups of samples that share common properties in high-content screening. Predicted pharmacodynamic effects were visualized and quantified in MCF-7 cells dosed with three target-specific compounds. Strong pairwise correlation between nuclear markers pHis-tone3 S28 , Ki-67, and p4E-BP1 T37/T46 in mitotic cells and anticorrelation with cell surface markers CD29, CD98, CD81, CD47, and EpCAM was demonstrated.

Supporting information
Additional supporting information may be found online in the Supporting Information section at the end of the article. .X S139 and cell area (scale bar = 10 lm). Fig. S3. Segmentation accuracy of nuclei and cytoplasmic contours of identified MCF-7 cells drug-treated using binary image overlap. The identified nucleus and cytoplasm from each cell as object is converted to a binary format called test 'image' and overlaid to the binary converted tiff image used for segmentation (DNA-Ir for nuclei identification and pan-keratin for cell identification) as 'ground truth'. Both binary images are overlapped to calculate statistics of the closeness from the test image to its ground truth. Performance segmentation statistics parameters is ranging on a scale from 0 to 1, where 1 means perfect overlap and 0 no overlap. Nuclei segmentation shows high precision (number of true positive pixels/(number of true positive pixels + number of false positive pixels) close to 1 over the different images. For cell segmentation the precision is slightly lower but still higher than 0.5 value. Fig. S4. Comparison of average Z-score expression level of p53, cyclin D3, CD47, CD49e and Cell-ID Intercalator-Ir between each control and drug treatment condition. Data are presented as mean AE SEM. ****P < 0.0001 (unpaired t-test).  Fig. S7. Comparative expression levels of protein markers Ki-67, cyclin D3, pH2A.X S139 , p53 and pHis-tone3 S28 (rows) per nuclear class and drug treatment condition (column) of MCF-7 cells. Bar charts represent non-normalized mean intensity (+SD) of each nuclear marker to its related class. P < 0.05 and P < 0.0001 (unpaired t-test) are included. Fig. S8. Heat maps of Pearson correlation coefficients of multiple parameters for non-treated and DMSOtreated controls ROIs on mitotic pHistone3 S28 cells.
Cluster highlighting of cell size parameters (black), surface markers (blue), Cell-ID iridium isotopes (green), pan-keratin and CK19 (orange) and nucleus markers (yellow). Fig. S9. Graphical representations of Pearson correlation coefficients for nuclei markers (A) and cytoplasmic markers (B) across all controls and drug treatments in the classified population of mitotic pHis-tone3 S28 MCF-7 cells. Hierarchical similarities heatmaps show all pairwise correlations values between each protein, with diagonal components highlighted as multiparametric clusters identified for each condition (Non-treated, blue; DMSO, red; EGF, green; etoposide, black; nocodazole, purple). Intra and inter-relationships between these clusters are shown as Hierarchical edgebundle visual graphs. Bundled lines connecting parametric protein nodes correspond to positive correlation values higher than 0.3. Table S1. List of chemical compounds used for MCF-7 cell treatment, 48 h exposure. Table S2. Metal-labeled antibody panel for cell surface markers. Table S3. Metal-labeled antibody panel for intracellular markers. Table S4. Mitotic index of compound-treated MCF-7 generated by support vector machine classification with Fast Gentle Boosting ruler. Population of pHistone3 S28 of all ROIs replicates per drug treatment and control were counted and divided by the total number of cells detected to determine the mitotic index of each condition.