Label-Free Classification of Apoptosis, Ferroptosis and Necroptosis Using Digital Holographic Cytometry

Featured Application: Digital holographic cytometry for classiﬁcation of programmed cell death. Abstract: Apoptosis, ferroptosis and necroptosis are three distinct forms of programmed cell death. Each of these pathways can be exploited to terminate cancer cells. One promising therapeutic strategy is to activate alternative programmed cell death pathways subsequent to cancer cells evolving mechanisms to evade apoptosis. However, the interplay between distinct programmed cell death pathways and cancer progression is complex and can paradoxically promote the disease. There is a need for high-throughput assays for real-time classiﬁcation of programmed cell death, both to further investigate these important biologic processes and to assess the case-by-case e ﬃ cacy of targeting each pathway in patient-derived tumor cells. Here, we sought to develop a label-free, live-imaging-based assay for classifying forms of programmed cell death with single cell resolution. We used digital holographic cytometry (DHC) to monitor human melanoma cells undergoing apoptosis, ferroptosis, and necroptosis. We developed and validated models that used DHC-derived features to classify each form of cell death with 91–93% accuracy in the test sets. We conclude that high-accuracy, high-throughput, label-free classiﬁcation of apoptosis, ferroptosis and necroptosis can be achieved with DHC.


Introduction
The ability to evade programmed cell death is one of the hallmarks of cancer, leading to therapeutic resistance or failure [1]. Although apoptosis is perhaps the most famous, other forms of programed cell death have been well-characterized [2]. Previous studies have demonstrated that cancer cells that are resistant to one type of cell death may be more vulnerable to another. Ferroptosis and necroptosis, for example, can be exploited to induce death in cancer cells that have successfully evaded apoptosis [3][4][5][6][7][8]. However, the interplay between each of these programed cell death pathways and their involvement in cancer progression is complex. For example, pan-caspase inhibitors have been shown to inhibit apoptosis while simultaneously inducing necroptosis [9]. Induction of necroptosis can be antitumorigenic in some settings or promote tumor progression in others [10][11][12][13][14]. In another example, vulnerability of melanoma cells to induction of ferroptosis is dependent on the phenotypic subtype [15]. Such observations highlight the need for further investigation and suggest that the therapeutic potential of targeting distinct programmed cell death pathways might require case-by-case assessment. Methods that efficiently and simultaneously quantify apoptosis, ferroptosis, and necroptosis would be valuable for further exploring the clinical utility of agents that regulate cell death pathways.
Apoptosis, ferroptosis and necroptosis are all active forms of cell death, meaning they comprise defined molecular pathways that drive organized termination. The key regulatory molecules that activate, inhibit, or orchestrate each process permit molecular classification of cell death pathways using a variety of methods [2]. For example, changes in expression, post-translational processing, or subcellular localization of key proteins can be monitored using immunological methods such as fluorescence microscopy, immunoblotting, flow analysis, or enzyme-linked immunoassays (ELISA). Other useful methodologies detect outer membrane permeabilization or mitochondria function and health-both of which change during the process of apoptosis [16]. While each of these approaches has distinct advantages and limitations, none permit both high-throughput classification of all three forms of cell death in real time and with single cell resolution.
Digital Holographic Cytometry (DHC) is a label-free technique that allows for live cell imaging and real time morphologic assessment [17]. DHC uses a laser that is split into two beams-a reference beam and an object beam. The reference beam is undisturbed, while the object beam is directed through the specimen. When the two beams are merged, an image can be generated through interpretation of the resulting interference pattern, with pixel intensity representing optical thickness-a function of both physical thickness and optical density. Cells are identified using standard image segmentation techniques, and individual cell features (optical volume, optical thickness, perimeter length, etc.) are derived. Quantitative phase images can be captured over a period of time, permitting tracking of DHC-derived features as cells undergo change, for example, due to exposure to a therapeutic agent. DHC has previously been used to distinguish different cell states based upon morphology, including identification of pre-apoptotic and dying cells [18][19][20][21][22].
In addition to molecular changes, a distinct set of morphological changes accompanies each form of programmed cell death. Apoptosis is associated with an increase in cell thickness, plasma membrane blebbing, decrease in cell volume, and chromatin condensation [23,24]. Ferroptosis is associated with iron-dependent membrane lipid oxidation and damage to the plasma membrane [24,25]. Necroptosis is associated with cell swelling, an increase in cell thickness, decrease in cell area, and eventual rupture of the plasma membrane [4,23,24]. We therefore reasoned that analysis of cellular morphologic change with DHC might permit label-free, live, and high-throughput classification of apoptosis, ferroptosis and necroptosis with single cell resolution. The goal of this study was to train and test a classifier to separate three types of cell death using DHC-derived cell features.

Cell Culture
The BRAF V600E -mutated, human-derived, metastatic melanoma cell line 501mel (CVCL_4633) was obtained from the University of California, San Francisco (a generous gift from Dr. Boris Bastian). It was obtained under ethical guidelines and informed consent was signed [26]. As all information regarding the donating patient has been de-identified, the use of the line for this study is not considered Human Subjects Research by the National Institutes of Health. Cells were cultured in growth media consisting of 10% Fetal Bovine Serum (VWR, Radnor, PA, USA), 1% L-Glutamine (Gibco Thermo Fisher Scientific, Gaithersburg, MD, United States), 1% Non-essential Amino Acid supplement (Gibco Thermo Fisher Scientific, Gaithersburg, MD, United States), and 1% Penicillin-Streptomycin (Gibco Thermo Fisher Scientific, Gaithersburg, MD, United States) in RPMI-1640 medium (Gibco Thermo Fisher Scientific, Gaithersburg, MD, United States). The cells were passaged twice a week at a ratio of 1 to 10. The cells were kept in a standard 37 • C humidified mammalian cell incubator with 5% CO 2 (Thermo Fisher Scientific, Waltham, MA, United States).
Preliminary dose-response experiments were conducted to determine the IC50 of each compound. Briefly, the cells were seeded in 6-well plates (Sarstedt, Nümbrecht, Germany) at a density of 60k For classification experiments, cells were cultured and treated as above, except that each plate contained a single well of each compound at the calculated IC50. The 6-well plate was then loaded onto the DHC platform and imaged via AppSuite. The parameters were set to 48 h, with images taken every 1 h from 20 random fields. The training and test sets were generated using independent passages, plates of cells, and freshly generated working solutions of each compound, and were conducted at different times. In two test set experiments, the erastin-treated well was lost due to technical reasons resulting in a sample size of two instead of four for this condition.

Digital Holographic Imaging and Cytometry
The M4 HoloMonitor (Phase Holographic Imaging, Lund, Sweden) uses a form of quantitative phase imaging called digital holographic microscopy to create a digitally reconstructed image of cells [17]. Images were segmented and analyzed using the Hstudio software (v2.7.5, Phase Holographic Imaging, Lund, Sweden). HStudio was used to derive 32 quantitative cell features for individual cells. For each condition, image stacks were monitored for cell death and time points selected at which cells were undergoing morphological changes prior to cell death. In control conditions, M-phase cells were excluded from the analyses.

Data Analysis and Statistics
Pearson's correlation coefficient was used to assess the association between all features. Blocks of features with significant correlation (r > 0.96) were refined to single features. One-way ANOVAs with Tukey's multiple comparisons test were used to determine significant differences in unique feature means between each of the treatments. The Tukey multiple comparison procedure was chosen as the goal was to adjust for all pairwise comparisons. Features that were significantly associated (adjusted p value < 0.0001) with at least two treatments were selected for further analysis. Decision trees were constructed by first identifying individual features that best separated each pair of conditions. Thresholds were identified by optimizing for accuracy of classification. Finally, the ordering of nodes was conducted manually through a rational consideration of the feature distributions. Sensitivity, specificity, positive predictive value (PPV), negative predictive value (NPV), and accuracy were calculated using standard methods [27].

Identification of Informative DHC-Derived Features
To investigate whether DHC-derived features could be used to distinguish cells undergoing apoptosis, ferroptosis and necroptosis, we generated four replicate datasets ( Figure 1). Human melanoma cells (501mel) were treated with compounds previously established as inducing apoptosis [28][29][30][31] (staurosporine, 0.1 µm), ferroptosis [8,25,[32][33][34][35] (erastin, 2 µm), or necroptosis [36][37][38][39][40][41] (shikonin, 1.5 µm). Cells treated with only DMSO were included as controls (not treated, NT). Cultures were live-imaged and DHC-derived features were generated for each cell. Two of the experiments contained all four conditions, whereas two were performed without erastin. One of the four-condition experiments was designated the Training Set and the rest as Test Sets. To first determine whether data generated from independent experiments were comparable, we evaluated the optical thickness (average pixel intensity in object) of all cells within the NT conditions ( Figure 2). Optical thickness is a fundamental measurement of DHC, from which other morphological features are derived, so we reasoned that comparable datasets should yield comparable distributions of this feature in NT conditions, where biological variability should be minimal. Indeed, optical thickness distributions were comparable across datasets, indicative of a stable platform and reproducible features. We next compared images from each condition ( Figure 3). We observed notable and distinct morphological changes upon treatment with each of the three compounds. The distribution of average optical thickness demonstrated somewhat greater variability in the treated conditions as compared to the NT conditions, indicative of the expected increase in biological variability ( Figure 2). We therefore investigated whether changes in morphological features were sufficiently consistent across experiments to classify each form of programed cell death despite biological variability.
Many of the morphological features generated by HStudio are highly correlated. In order to avoid over-training, we sought to limit the features considered in classification models. We calculated the correlation coefficient for each combination of features in the Training Set ( Figure 4). For features with perfect or near-perfect (r > 0.96) correlations, only one was selected for further analysis. The choice of feature was either arbitrary or, where appropriate, driven by ease of biological interpretation (for example, optical thickness was selected over phaseshift). We compared the distribution of all selected features across the four conditions in the Training Set. Nine features demonstrated a significant difference in distribution (adj. P < 0.0001) between at least two conditions ( Figure 5 and Supplemental Tables 1-4). These nine features were selected as the basis for two approaches to classifier generation ( Figure 1).    Many of the morphological features generated by HStudio are highly correlated. In order to avoid over-training, we sought to limit the features considered in classification models. We calculated the correlation coefficient for each combination of features in the Training Set (Figure 4). For features with perfect or near-perfect (r > 0.96) correlations, only one was selected for further analysis. The choice of feature was either arbitrary or, where appropriate, driven by ease of biological interpretation (for example, optical thickness was selected over phaseshift). We compared the distribution of all selected features across the four conditions in the Training Set. Nine features demonstrated a significant difference in distribution (adj. P < 0.0001) between at least two conditions ( Figure 5 and Supplemental Tables 1-4). These nine features were selected as the basis for two approaches to classifier generation (Figure 1).  Many of the morphological features generated by HStudio are highly correlated. In order to avoid over-training, we sought to limit the features considered in classification models. We calculated the correlation coefficient for each combination of features in the Training Set (Figure 4). For features with perfect or near-perfect (r > 0.96) correlations, only one was selected for further analysis. The choice of feature was either arbitrary or, where appropriate, driven by ease of biological interpretation (for example, optical thickness was selected over phaseshift). We compared the distribution of all selected features across the four conditions in the Training Set. Nine features demonstrated a significant difference in distribution (adj. P < 0.0001) between at least two conditions ( Figure 5 and Supplemental Tables 1-4). These nine features were selected as the basis for two approaches to classifier generation (Figure 1).

Optical Volume and Optical Thickness Distinguish Form of Programmed Cell Death
In our first approach, we sought to train a classifier that, when presented with a dying cell, could predict the form of programmed cell death. For this reason, healthy cells were excluded from this analysis (Figure 1). Cells undergoing necroptosis presented a significant increase in thickness, as previously observed using other systems, whereas cells undergoing ferroptosis presented substantially decreased volume ( Figure 6). Consequently, the form of programmed cell death could be classified with high accuracy using just these two features (Figure 7). We observed an overall accuracy of 89% of the Training Set. When the trained model was then presented with unmodified data (i.e., not normalized for batch effects) from each of the three Test Sets, accuracies remained high, ranging from 87 to 97% (Supplemental Table 5). Accuracy were high for the datasets lacking the erastin condition, indicating the performance of the model was not dependent on receiving equally-sized groups as input. When all Test Set data were combined, the model performed with an overall accuracy of 91% (Figure 7).

Optical Volume and Optical Thickness Distinguish Form of Programmed Cell Death
In our first approach, we sought to train a classifier that, when presented with a dying cell, could predict the form of programmed cell death. For this reason, healthy cells were excluded from this analysis (Figure 1). Cells undergoing necroptosis presented a significant increase in thickness, as previously observed using other systems, whereas cells undergoing ferroptosis presented substantially decreased volume ( Figure 6). Consequently, the form of programmed cell death could be classified with high accuracy using just these two features (Figure 7). We observed an overall accuracy of 89% of the Training Set. When the trained model was then presented with unmodified data (i.e., not normalized for batch effects) from each of the three Test Sets, accuracies remained high, ranging from 87 to 97% (Supplemental Tables 5). Accuracy were high for the datasets lacking the erastin condition, indicating the performance of the model was not dependent on receiving equallysized groups as input. When all Test Set data were combined, the model performed with an overall accuracy of 91% (Figure 7).

3.3.A four-Node Classifier Distinguishes Form of Programmed Cell Death from Healthy Cells
Although encouraging, our first approach is only useful when classifying a cell pre-identified as undergoing cell death. A classifier that could identify different forms of programmed cell death from the diverse morphologies adopted by healthy cells could be a more versatile tool. This is a more challenging demand due, in large part, to the morphological heterogeneity of both healthy cells and

A four-Node Classifier Distinguishes Form of Programmed Cell Death from Healthy Cells
Although encouraging, our first approach is only useful when classifying a cell pre-identified as undergoing cell death. A classifier that could identify different forms of programmed cell death from the diverse morphologies adopted by healthy cells could be a more versatile tool. This is a more challenging demand due, in large part, to the morphological heterogeneity of both healthy cells and staurosporine-treated cells (Figures 3 and 5). Although we noted several features with distributions that were significantly different across each of the four conditions (for example, Figure 5, optical thickness average (ave)), the range of values in each condition were sufficiently overlapping to prohibit accurate classification of individual cells based upon any one or two features. However, when considering four features, we succeeded in constructing a decision tree that generated overall accuracies of 92% and 93% when presented with the Training and Test Sets, respectively (Figure 8). Necroptosis and ferroptosis remained identifiable with 97% and 95% accuracies based upon optical thickness and optical volume alone. Of the remaining unclassified cells, apoptotic cells were identified with an accuracy of 88% by a combination of large fluctuation in phaseshift (potentially due to the established process of membrane blebbing) and a decrease in length. Similar to the two-node classifier, presentation of each Test Set individually resulted in a tight range of high overall accuracies (91-98%), even in the absence of erastin, indicating the performance of the model was not dependent on a single condition or Test Set (Supplemental Table 6).
Appl. Sci. 2020, 10, x FOR PEER REVIEW 8 of 13 due to the established process of membrane blebbing) and a decrease in length. Similar to the twonode classifier, presentation of each Test Set individually resulted in a tight range of high overall accuracies (91-98%), even in the absence of erastin, indicating the performance of the model was not dependent on a single condition or Test Set (Supplemental Tables 6).

Discussion
Here we present the application of DHC to the classification of human melanoma cells undergoing three active forms of programmed cell death-apoptosis, ferroptosis, and necroptosis. This approach to monitoring and distinguishing forms of programmed cell death has several advantages over other methods. First, digital holographic imaging is label-free and therefore does not require the manipulation of the cells through addition of dyes or introduction of exogenous genetic material. Second, the assay is compatible with live-imaging, permitting real-time monitoring of cell viability and assessment of the kinetics of each form of programmed cell death. Third, the

Discussion
Here we present the application of DHC to the classification of human melanoma cells undergoing three active forms of programmed cell death-apoptosis, ferroptosis, and necroptosis. This approach to monitoring and distinguishing forms of programmed cell death has several advantages over other methods. First, digital holographic imaging is label-free and therefore does not require the manipulation of the cells through addition of dyes or introduction of exogenous genetic material. Second, the assay is compatible with live-imaging, permitting real-time monitoring of cell viability and assessment of the kinetics of each form of programmed cell death. Third, the approach has single-cell resolution, and can therefore provide information on the percentages of a population that undergo different forms of programmed cell death. Finally, DHC is a non-invasive and non-endpoint assay, such that once the analysis is complete, the cells remain unadulterated and can be used for further molecular analyses. We expect this assay will be a valuable addition to the toolbox for investigating programed cell death, and that it will be particularly useful for assessing the dynamics and interplay of different forms of cell death in the same culture.
The strengths of this assay make it particularly well suited for real-time assessment of the therapeutic response of cells directly isolated from human tumors. The potential clinical utility of assessing tumor-specific sensitivity to a collection of potential single agent or combination therapies using in vitro approaches has long been appreciated [42]. However, several demanding requirements impede the successful implementation of this approach as a widespread clinically valuable tool for personalized medicine. First, in order to provide efficient turnaround time and to most accurately capture the in situ cell state, tumor-derived cells must be minimally expanded in vitro. As a label-free approach that requires few cells per condition, the assay presented here requires little in vitro manipulation and selection. Second, an ideal assay would capture the phenotypic heterogeneity of both the naïve tumor cells and their response to therapeutic agents. Capturing heterogeneity of response is crucial for assessing the adoption of "persister states"-cells that survive an initial toxic treatment, often leading to therapeutic resistance, but also increasing susceptibility to other therapeutic agents [2]. The assay presented here could, for example, simultaneously compare not only the efficacy of two different agents, but also determine whether phenotypic subsets of tumor cells were predominately affected, and whether different forms of programmed cell death were triggered. Such information would be instructive when selecting candidate combination therapies. For these reasons, DHC platforms have already been implemented in the high-throughput real-time screening of tumor-derived cells [43,44]. The approach presented here offers additional granularity to the information regarding cell death provided by these systems.
We present here the classification of three forms of programmed cell death. However, other major forms of programmed cell death, most notably, pyroptosis, have also been described 2 . We limited our analyses to the forms of cell death that we could confidently and uniformly induce based upon prior work. It is possible that similar approaches could yield models that also identify pyroptosis. One limitation of this study is the inclusion of only single agents to induce each form of cell death. Each agent used in this study (staurosporin, erastin, shinkonin) is an established and classic inducer of the associated form of cell death [25,[29][30][31][32][33][34][35][37][38][39][40][41]. However, it remains possible that the morphological changes measured here are specific to the agent, not the form of cell death. Similarly, an additional limitation to the work presented here is the analysis of a single human melanoma line. As we and others have applied DHC to the classification of different cellular phenotypes [18,[45][46][47][48][49][50][51], one emerging concern is that single models cannot be trained that accurately classify diverse cell lines or conditions. Indeed, as the approach is based exclusively on morphology, it should be expected that cells with distinct morphologies cannot be assessed by the same model with any degree of accuracy. For these reasons, we emphasize that, as for any assay, the optimal parameters must be determined for each cellular system, experimental condition, reagent, and biological question. We do not encourage the blind application of the exact decision nodes and thresholds defined here to cell lines other than human melanoma cells. Rather, it is our intent, that by visualizing each step of model development, we have provided a workflow for generating high-accuracy classifiers using DHC-derived features that is accessible to all cell biologists interested in applying this method to their cellular systems.

Conclusions
We conclude that three forms of programmed cell death-apoptosis, ferroptosis and necroptosis-can be accurately classified using DHC-derived features, label-free, in real-time, and with single-cell resolution.
Supplementary Materials: The following are available online at http://www.mdpi.com/2076-3417/10/13/4439/s1, Table S1: Nine DHC-derived features from Training Set. Table S2: Nine DHC-derived features from Test Set 1. Table S3: Nine DHC-derived features from Test Set 2. Table S4: Nine DHC-derived features from Test Set 3. Table  S5: Performance of two-node decision tree with each individual Test Set. Table S6: Performance of four-node decision tree with each individual Test Set.