The Quantitative-Phase Dynamics of Apoptosis and Lytic Cell Death

Cell viability and cytotoxicity assays are highly important for drug screening and cytotoxicity tests of antineoplastic or other therapeutic drugs. Even though biochemical-based tests are very helpful to obtain preliminary preview, their results should be confirmed by methods based on direct cell death assessment. In this study, time-dependent changes in quantitative phase-based parameters during cell death were determined and methodology useable for rapid and label-free assessment of direct cell death was introduced. The goal of our study was distinction between apoptosis and primary lytic cell death based on morphologic features. We have distinguished the lytic and non-lytic type of cell death according to their end-point features (Dance of Death typical for apoptosis versus swelling and membrane rupture typical for all kinds of necrosis common for necroptosis, pyroptosis, ferroptosis and accidental cell death). Our method utilizes Quantitative Phase Imaging (QPI) which enables the time-lapse observation of subtle changes in cell mass distribution. According to our results, morphological and dynamical features extracted from QPI micrographs are suitable for cell death detection (76% accuracy in comparison with manual annotation). Furthermore, based on QPI data alone and machine learning, we were able to classify typical dynamical changes of cell morphology during both caspase 3,7-dependent and -independent cell death subroutines. The main parameters used for label-free detection of these cell death modalities were cell density (pg/pixel) and average intensity change of cell pixels further designated as Cell Dynamic Score (CDS). To the best of our knowledge, this is the first study introducing CDS and cell density as a parameter typical for individual cell death subroutines with prediction accuracy 75.4% for caspase 3,7-dependent and -independent cell death.

colourimetric endpoint visualization of the analyzed parameter and belong among indirect assays. Such indirect endpoint analyses are prone to the misleading results as we have shown in our previous work 5 . Although morphological aspects of cell death are not generally recommended to determine cell death subroutines 1 , it would be a mistake to completely ignore them. Recent progress in Quantitative Phase Imaging techniques (QPI) has enabled the observation of time-dependent subtle changes unrecognizable to the naked eye (such as cell mass distribution) on micrographs. These changes in cell mass distribution, cell density, micro-blebbing of the cell membrane, nuclear shape, homogeneity of cell content, and many other parameters, which can be typical for individual cell death subroutines, can be observed without fixation, labelling or cell harvesting.
In this article, we demonstrate a methodology useable for rapid assessment of direct cell death that is based on NCCD recommendations (Fig. 1). The goals of our study were (a) estimation of time point of cell demise and (b) distinction between apoptosis and primary lytic cell death based on morphologic features. We were able to distinguish the lytic and non-lytic type of cell death according to their end-point features (Dance of Death typical for apoptosis versus swelling and membrane rupture typical for all kinds of necrosis common for necroptosis, pyroptosis, ferroptosis and accidental cell death). Using advanced quantitative label-free phase imaging (QPI), we were also able to observe typical dynamical changes of cell morphology during both caspase 3, 7-dependent and -independent cell death subroutines and to determine the moment of cell death. Dynamical, time-dependent  (2) prediction of timepoint of death for each cell. Long-short term memory (LSTM) neuronal network is used for this step and all features mentioned in the step 2 were included for this prediction. (3) Cell death type classification is based on two parameters: density and cell dynamic score (CDS). Based on them, Support Vector Machine classifier (SVM) is used for cell death type classification as "apoptotic" or "lytic". The input parameters for this analysis are: (i) interval minus 10 h to timepoint of death (predicted in the previous step), and (ii) manually annotated death type (apoptosis vs lytic cell death). changes in cell mass distribution maps were used for the label-free distinction between typical morphological patterns appropriate for caspase-dependent and caspase-independent cell death subroutines.
As model cells, prostate cancer cell lines DU145, LNCaP, and benign cell line PNT1A established by immortalization of normal adult prostatic epithelial cells were used. These cell lines were not selected because of the clinical relevance, but rather because of their radically different size and accessibility for the automatic segmentation algorithm. We expected cellular size to be the factor influencing the rate of morphological changes during cell death. Cell death was triggered by doxorubicin and staurosporine, two commonly known inducers of apoptosis with a distinct induction of caspase cleavage and consequently a distinct morphological manifestation 6,7 . Cell death induced by black phosphorus (BP) was detected as an example of difficult detection conditions 8 .
Cell tracking is an essential step in cell image analysis. Even though many methods exist 9 , there is no universal and sufficiently robust method applicable to QPI, especially for touching cells and if correct tracking through whole image sequence is needed. Simple available tools like TrackMate 10 have shown to be insufficient, thus we have developed a new tracking method tailored for our dataset.
The number of existing methods for the detection of cell death in label-free time-lapse images is very limited. In 11 authors detect cell death event in Phase Contrast Microscopy using features time-series (size, roundness, speed etc.) and classify each time point as alive or death with transductive support vector machine. In comparison, our technique based on quantitative microscopy data uses more features and more advanced Long Short-Term Memory (LSTM) neural network.
There are several techniques for static cell image classification including extraction of features and application of classifier 12,13 or application of convolutional neural network 14 . These techniques can be also possibly applied for a distinction of different types of cell death, however, we decided for using only two features for this distinction in order to make our results easily interpretable, while also incorporating new features describing cell dynamics. cell culture and cultured cell conditions. LNCaP cell line was established from a lymph node metastase of the hormone-refractory patient and contains a mutation in the androgen receptor (AR) gene. This mutation creates a promiscuous AR that can bind to different types of steroids. LNCaP cells are AR-positive, PSA-positive, PTEN-negative and harbor wild-type p53 15,16 . PNT1A is an immortalized non-tumorigenic epithelial cell line. PNT1A cells harbour wild-type p53. However, SV40 induced T-antigen expression inhibits the activity of p53. This cell line had lost the expression of AR and prostate-specific antigen (PSA) 17 . DU-145 cell line is derived from the metastatic site in the brain and contains P223L and V274F mutations in p53. This cell line is PSA and AR-negative and androgen independent 18 . All cell lines used in this study were purchased from HPA Culture Collections (Salisbury, UK) and were cultured in RPMI-1640 medium with 10% FBS. The medium was supplemented with antibiotics (penicillin 100 U/ml and streptomycin 0.1 mg/ml). Cells were maintained at www.nature.com/scientificreports www.nature.com/scientificreports/ 37 °C in a humidified (60%) incubator with 5% CO 2 (Sanyo, Japan). Cell lines were not tested on mycoplasma contamination.
Correlative time-lapse quantitative phase-fluorescence imaging. QPI and fluorescence imaging were performed by using multimodal holographic microscope Q-PHASE (Telight a.s., Brno, Czech Republic). www.nature.com/scientificreports www.nature.com/scientificreports/ To determine the amount of caspase-3/7 product accumulation, cells were loaded with 2 µM CellEvent TM Caspase-3/7 Green Detection Reagent (Life Technologies, Carlsbad, CA, USA) according to the manufacturer's protocol and visualized using FITC 488 nm filter. To detect the cells with a loss of plasma membrane integrity, cells were stained with 1 ug/ml propidium iodide (Sigma Aldrich Co., St. Louis, MO, USA) and visualized using TRITC 542 nm filter. Nuclear morphology and chromatin condensation were analyzed using Hoechst 33342 nuclear staining (ENZO, Lausen, Switzerland) and visualized using DAPI 461 nm filter. Cells were cultivated in Flow chambers μ-Slide I Lauer Family (Ibidi, Martinsried, Germany). To maintain standard cultivation conditions (37 °C, humidified air (60%) with 5% CO 2 ) during time-lapse experiments, cells were placed in the gas chamber H201 -for Mad City Labs Z100/Z500 piezo Z-stages (Okolab, Ottaviano NA, Italy). To image enough cells in one field of view, Nikon Plan 10/0.30 was chosen. For each of three cell lines and each of three treatments, seven fields of view were observed with the frame rate 3 mins/frame for 24 or 48 h respectively.
Holograms were captured by CCD camera (XIMEA MR4021 MC-VELETA), fluorescence images were captured using ANDOR Zyla 5.5 sCMOS camera. Complete quantitative phase image reconstruction and image processing were performed in Q-PHASE control software. Cell dry mass values were derived according to 19 and 20 from the phase (Eq. (1)), where m is cell dry mass density (in pg/μm 2 ), ϕ is detected phase (in rad), λ is wavelength in μm (0.65 μm in Q-PHASE), and α is specific refraction increment (≈0.18 μm 3 /pg). All values in the formula except the ϕ are constant. The value of ϕ (Phase) is measured directly by the microscope. Integrated phase shift through a cell is proportional to its dry mass, which enables studying changes in cell mass distribution 20 . cell dry mass tracking. The custom method for automatic cell tracking and measuring of selected features was developed and implemented in MATLAB. Although tracking is not the main contribution of this article, it is a necessary step for our analysis and there is no available method suitable for QPI. The proposed tracking method consists of two main steps -nuclei tracking followed by expansion of each nucleus region to the whole cell.
Nuclei tracking is done by Gaussian Mixture Model (GMM) fitting to nuclei image (Hoechst 33342) in each frame by Expectation-Maximization (EM) algorithm (a similar method is used for segmentation in 21 ). For each frame, Gaussians are provided from the previous frame and their parameters are updated by several steps of EM. Before the application of GMM, background pixels are eliminated by segmentation using Maximally Stable Extremal Region (MSER) 22 . Gaussians parameters are optimized for each nuclei cluster separately. In the first frame, positions of Gaussians are initialized manually and the covariance matrices are initialized as an identity matrix. The whole algorithm is summarized in Fig. 2 and for more details see 23 . This method does not recognise division of cells. For the next analysis, we use only cells, which occur in the whole image sequence without division or joining with other cell (which appears mainly due to tracking and segmentation errors).
The previous step provides tracked nuclei, but for extraction of most of the cell features, we need cell segmentation. However, nuclei can be used as seeds for segmentation in each frame. We can use simple thresholding for segmentation of foreground in QPI cell image, where the threshold value can be the same in all frames, which provides sufficient results for a distinction between foreground and background, but a separation of single cells is very problematic. We use tracked nuclei as seeds for proper division of the foreground binary mask F t (obtained by tresholding) to the single cells binary mask S t . Seeded watershed 24 is used on the negative of the original image O t with nuclei tracking result as seedpoints. Watershed results (boundary lines) are then used for division of the foreground mask.
However, the simple application of the watershed algorithm leads to a fast mass exchange between cells due to segmentation errors, which is very undesirable for precise feature extraction. For this reason, we introduced a simple modification named Movement Regularized Watershed (MRW), in order not to allow dramatic contour www.nature.com/scientificreports www.nature.com/scientificreports/ changes between frames. This can be achieved by incorporating the mask from the previous frame to the actual frame watershed calculation. This can be done by modifying O t and F t as where − S x y ( , )   www.nature.com/scientificreports www.nature.com/scientificreports/ morphological changes, but it is not much dependent on the segmentation quality, because the same mask is used in both frames. All these features were evaluated in all frames, where the result is a set of signals describing the cell behaviour in time.
Label-free algorithm for cell death detection. Besides the significant mass decrease of dead cells, various feature evolvement types were observed during detected cell deaths, which complicate the expert specification of the detected phenomenon, thus machine learning approach was chosen for this task. Bidirectional Long Short-term Memory (BiLSTM) networks 25 have shown to be very successful for signal classification and regression tasks, thus it is suitable for the purpose of detection of cell death in analysed cell signals. Manual annotation of cell death timepoint is converted to a signal with a Gaussian curve at a death timepoint. The proposed approach is inspired by the detection of correlation filters 26 , where the aim is to regress the Gaussian curve on the desired www.nature.com/scientificreports www.nature.com/scientificreports/ position, where Gaussian represents uncertainty in position. Long-short term memory (LSTM) network has been trained for the regression of Gaussian curves created on the time of cell death (with sigma 50), where the whole method is summarized in Fig. 1.
According to 27 , the network was set to two BiLSTM layers (with 100 units) and 3 fully connected layers (with 100, 50 and 100 neurons) with ReLu and dropout with probability 0.5. The whole network was optimized using Mean Square Error loss and ADAM optimizer 28 with learning rate 10 −3 , β 1 = 0.9, β 2 = 0.999, gradient clipping to norm 1 29 , weight decay 10 −3 , batch size 256 and 40 epochs. We augment our training dataset with random clipping (shortening each signal by at most 1/3 of length). One network was trained for all 9 cell lines/experiments evaluated by cross-validation, where one FOV from all experiments was used for testing and the rest for training. Cell death time was identified as a maximum in the network response with a value higher than the chosen threshold 0.4. Quality of cell detection is evaluated in terms of accuracy, where detection is considered as correct if it detects a death in a ±5 h (±100 frames) window from ground truth cell death or if death is not detected and the cell is labelled as alive (at the end of the experiment).
Cell death type identification. We observed two distinct dynamical patterns of cell mass manipulation during cellular death, where each cell was manually labelled as one of these types. To quantify these morphological types we extracted average values of extracted features in 200 frames (10 h) before cell death (shorter time window is used if more timepoints are not available). We trained linear Support Vector Machine (SVM) classifier for automatic classification of these cell death types. SVM classifier was trained for each cell line (one for all three treatments) because cell lines are morphologically different. All possible subsets of features were tested and classification accuracy was evaluated. Only density and CDS were finally used because accuracy does not increase significantly with the addition of other features (accuracy 75.4% and 76.2% for two and three features, respectively), where two features can be comprehensibly visualized. These two features have also a distinct biological meaning and its difference in an average aligned signal can be confirmed in Fig. 1.
Furthermore, the correctness of the classification was confirmed by the analysis of nucleus shrinkage and by Casp 3,7-PI signal onset delay in fluorescence data. Nucleus shrinkage was measured as an average Hoechst 33342 brightness in 200 frames (10 h) before cell death (based on nuclei mask produced by GMM tracking). Casp 3,7-PI signal onset delay was measured as a time difference (delay) between the time when Casp 3,7 signals reach 1/3 of its average value after cell death (average in widnow of 200 frames) and the cell death time (which corresponds to a steep PI signal increase). The output of this analysis is shown in Fig. 3d,e.

Results
Automatic cell death detection. After induction of cell death by using 0.5 µM staurosporine, 0.1 µM doxorubicin, or 400 µg/mL of BP, respectively, 11 morphological parameters detectable by QPI (see Fig. 1) and their real-time values were collected. Effector caspase 3/7 activity, nuclear morphology and membrane integrity (propidium iodide signal; PI) were verified using time-lapse fluorescent microscopy to estimate ground truth for cell death detection based on cell mass parameters. The decrease in cell dry mass (pg), cell density (pg/pixel), and the Cell Dynamic Score (CDS; see the methodological part for equation) measured by QPI was in clear relationship with the onset of propidium iodide and caspase signalling. Furthermore, dead cells showed no longer any significant changes in CDS. (see Fig. 4). Based on QPI features, automatic and label-free detection of IT 50 (time of half-maximal inhibition effect for given treatment concentration) can be performed (see Fig. 5). Its evaluation is based on automatic cell detection of cell death using the LSTM network, where we can easily accumulate the number of dead cells assuming that all cells are alive at the start of the image sequence. This parameter is important because IC 50 values can be significantly different in different time points (such as 24 h versus 48 h treatment 30 ). For IT 50 of tested compounds and cell lines see Fig. 5. Automatic detection of cell death showed 76% accuracy compared to the manual detection of cell death based on the fluorescent PI signal and morphological criteria visible to the naked eye.
Cell density and cell dynamic score reflect subroutine of cell death. Based on QPI data, we were able to distinguish two different subroutines of cell death: a) cells having high cell density and intensive blebbing of the plasma membrane (high CDS) followed by plasma membrane rupture and b) cells having low cell density and low CDS before plasma membrane rupture (see Fig. 6 and Supplementary Videos 1-2). These variants occupied extreme positions on the CDS versus cell density plot (see Fig. 3b and Supplementary Fig. 1). Time-lapse images and QPI features of these cells are shown in Fig. 3a and Supplementary Figs. 3 and 8. Based on fluorescent data referring caspase 3, 7 activity, nuclear morphology and membrane integrity (PI signal); see Fig. 6 and Supplementary Figs. 2-12, the a) type of cell death with high CDS and high cell density is bona fide apoptosis because PI signal was delayed over caspase signals and shrinkage of the nucleus was apparent (Fig. 3e). Moreover, cells with high CDS and high cell density were only rarely presented in z-VAD-FMK treated population, see Fig. 3c; Z-VAD-FMK is a cell-permeable pan-caspase inhibitor that irreversibly binds to the catalytic site of caspase proteases and inhibits apoptosis. On the other hand, the b) type of cell death is bona fide lytic cell death as PI and caspase signals were displayed simultaneously and no significant shrinkage of the nucleus was apparent. Cells near the dividing line showed a)-like and b)-like features to varying degrees. Casp 3,7-PI signal onset delay and nucleus shrinkage (described in Methods) were quantified for both automatic and manual labels of cell death types (see Fig. 3d and Supplementary Table 1), which confirmed the existence of these two cell death subroutines. Furthermore, apoptosis and lytic cell death were visually detected by the expert according to the fluorescent signals and eye-visible cell morphology. The automatic label-free distinction of cell death subroutine showed 75.4% accuracy (average of all cells and treatments) compared manual distinction of cell death type based on fluorescent signals and eye-visible morphological criteria; (see Fig. 3f). Interestingly, no cells showed a gradual cell rounding (2020) 10:1566 | https://doi.org/10.1038/s41598-020-58474-w www.nature.com/scientificreports www.nature.com/scientificreports/ and loss of surface contact during staurosporine treatment. Other treatments caused such phenomena relatively often. For illustration see Fig. 7.

Discussion
The basic requirement for a good method of cell death detection is the ability to detect a robust parameter in a highly reproducible and if possible inexpensive manner under changing conditions caused by the cell heterogeneity. In this study, we suggest that morphological and dynamical features extracted from the QPI cell image are suitable for cell death detection. QPI enables the time-lapse observation of subtle changes in the cell mass distribution unrecognizable by the naked eye, and therefore provides detailed information on cell morphology and cell mass topography during cell death. Cell dry mass can be calculated directly from phase values detected in each pixel [31][32][33] . Cell dry mass is the direct result of biosynthetic and degradative processes within a cell and is, therefore, a promising indicator of cell fate including cell death 34,35 . The distribution of cell mass during the processes of cell death is significantly changing in time with a typical steep decrease of the cell mass due to the rupture of the plasma membrane or cell fragmentation. This phenomenon is bona fide universal for all cell types and was also observed in target cells after contact with a cytotoxic T-cells 34 . Plasma membrane rupture can be otherwise measured by quantification of the release of intracellular enzymes from the cell into the cell culture medium. Nevertheless, the enzymatic activity of these enzymes can be seriously affected by differences in the pH between intracellular environment and culture medium and by time spent in the extracellular space 3 .
Based on QPI data, we were also able to distinguish two specific subroutines of cell death (Fig. 3a-e). While the cell death subroutine a) (bona fide apoptosis) showed a high cell density and strong fluctuation in CDS before cell death; during the subroutine b) (bona fide lytic cell death) the cell density was low and CDS was relatively stable after the initial decline. These variants occupied extreme positions on the CDS versus cell density plot. Cells near the dividing line (see Fig. 3b) probably succumb to another subroutine of cell death such as pyroptosis or necroptosis, as neither staurosporine nor doxorubicin is a specific apoptosis inductor 36,37 . The a) type of cell death was characteristic by a decrease in the cell area, a gradual cell rounding and a loss of surface contact (this phenomenon was not observed in the case of staurosporine treatment) followed by membrane blebbing (known as "dance of death"); see Fig. 3 and Supplementary Fig. 2 to 12 38 . Since membrane blebbing during apoptosis results from caspase-mediated activation of ROCK I, it can be assumed that apoptosis was involved in these cases 39 . At the final stages of apoptosis, the actin cytoskeleton is degraded 40,41 . Due to the hydrophobic nature of the apoptotic bodies, they undergo the plasma membrane fusion in the culture medium. The membrane of this post-apoptotic body subsequently cracks and culminates in secondary necrosis (in our case depicted by a steep decrease of cell mass) in the absence of phagocytes 42,43 . The b) type of observed cell death presented persisting large cytoplasmic membrane blebs or multiple small blebs and cell swelling leading to the final cell membrane rupture. This type of dying cells was adherent during the whole process. It can be assumed that lytic cell death was involved in these cases 3 . Although several studies based on QPI detection of cell death have been published earlier, they do not include the distinguishing of specific subroutines of cell death and have failed to capture the entire process of cell death, including the early stages due to short intervals of QPI capturing [44][45][46] . From an image analysis point of view, we introduce a new powerful technique for cell tracking and segmentation, capable of robustly track cell thought the whole image sequence. We also introduced a completely new method for the detection of cell death in time-series of measured cell features using LSTM neural network. To the best of our knowledge, this is also the first study introducing CDS and cell density as a parameter typical for individual cell death subroutines.

Data availability
The Matlab code is available at GitHub (https://github.com/tomasvicar/CellDeathDetect). Annotated quantitative phase image dataset (with cell tracking masks and labels of cell death timepoints and types of death) used in the manuscript is available at Zenodo repository (www.zenodo.org), https://doi.org/10.5281/zenodo.2601562.