Research ArticleCELL BIOLOGY

Single-cell cytometry via multiplexed fluorescence prediction by label-free reflectance microscopy

See allHide authors and affiliations

Science Advances  15 Jan 2021:
Vol. 7, no. 3, eabe0431
DOI: 10.1126/sciadv.abe0431

Abstract

Traditional imaging cytometry uses fluorescence markers to identify specific structures but is limited in throughput by the labeling process. We develop a label-free technique that alleviates the physical staining and provides multiplexed readouts via a deep learning–augmented digital labeling method. We leverage the rich structural information and superior sensitivity in reflectance microscopy and show that digital labeling predicts accurate subcellular features after training on immunofluorescence images. We demonstrate up to three times improvement in the prediction accuracy over the state of the art. Beyond fluorescence prediction, we demonstrate that single cell–level structural phenotypes of cell cycles are correctly reproduced by the digital multiplexed images, including Golgi twins, Golgi haze during mitosis, and DNA synthesis. We further show that the multiplexed readouts enable accurate multiparametric single-cell profiling across a large cell population. Our method can markedly improve the throughput for imaging cytometry toward applications for phenotyping, pathology, and high-content screening.

INTRODUCTION

Cell morphology features are powerful phenotypical readouts, which have been the basis for pathology for decades. They are also the underlying mechanisms for varieties of imaging cytometry and high-content screening platforms to characterize pathological changes and responses to drug treatments (1). The most widely used approach for imaging readouts are fluorescence labels that highlight specific subcellular components or cell functions through immunofluorescence (IF), fluorescent reporter cells, or dyes. However, the throughput of these approaches is fundamentally limited by the physical process of labeling. The IF staining is labor intensive and generally requires cell fixation that does not allow kinetic observations of live cells over time. Fluorescent reporter cells are permissive for longitudinal live cell imaging. However, the process of gene editing and validation takes a substantial amount of time and can be difficult to introduce multiple markers within the same cells for multiplexed analysis. Regardless of the fluorescence labeling approaches, the overlapping of the fluorescence emission spectra further limits the multiplexing capability. To alleviate these limitations, here, we develop a label-free single-cell cytometry that is highly multiplexed and can forgo the physical staining via a deep learning (DL)–augmented digital labeling method.

Our work relies on the premise that label-free scattering-based microscopy captures rich structural information and can be effective to characterize cell morphological features (2). Bright-field, phase contrast, and differential interference contrast microscopy have been routinely used for observing and quantifying cell morphology (3). Scattering-based microscopy and tomography techniques have been increasingly used to reconstruct cellular structures (4). Of particular interest is the reflectance-mode microscopy that provides exquisite sensitivity in detecting nanoscale structural changes beyond the diffraction limit (57). By capturing backscattering signals, reflectance imaging provides access to the highest spatial frequency components in the reciprocal Fourier space and thus can provide higher structural contrast than the transmission techniques (2). Our recent work shows that backscattering signals allow resolving finer details than the transmission counterparts (8, 9). In this work, we further leverage the higher sensitivity provided by the reflectance-mode microscopy and demonstrate how enriched label-free information allows predicting highly accurate subcellular structural features.

The framework of this study is summarized in Fig. 1. The angle-dependent backscattering features are captured with dark-field oblique illumination and paired with IF images (Fig. 1A). Using the IF images as the ground truths, multiple DL models are independently trained for individual IF labels (Fig. 1B). Once all the models are trained, we perform digital multiplexing by feeding the same label-free input to each network and make different IF predictions in parallel (Fig. 1C). By doing so, multiple subcellular structures and cell states can be revealed simultaneously without physical labeling. While previous works have shown that DL models can disentangle the complex structures captured in the label-free data and make in silico fluorescence labeling with high accuracy (1012) or holistically capture “hidden” structural features that are not easily perceived or described (1322), these results are fundamentally limited by the weak structural contrast from the transmission modes that contain only forward scattering information. By exploiting the enhanced resolution and sensitivity in the backscattering data, we demonstrate a marked increase in the fluorescence prediction accuracy with up to three times improvement as compared to the current state of the art.

Fig. 1 Overview of the DL-augmented label-free cytometry technique.

(A) A multimodal LED array reflectance microscope is developed to acquire co-registered label-free reflectance and fluorescence images. Reflectance images from oblique dark-field illumination and computed drDPC contain rich morphological information and are the multichannel input to our DL model. Two-channel epifluorescence images are acquired on the same sample to serve as the ground truth for training our DL model. (B) Individual DL models are trained independently with paired label-free and IF images. The saliency map is used to reveal specific label-free features captured by the model to perform the transformation. (C) To perform digital multiplexed predictions, the same reflectance input is fed to each network and makes six different IF predictions in parallel.

One distinct contribution of our work is to advance beyond the prediction of fluorescence images and to demonstrate accurate structural phenotyping and quantitative single-cell cytometry using digitally multiplexed fluorescence images. We show that our DL model can correctly capture and predict characteristic subcellular features during the cell cycle, including morphological changes of nuclei and Golgi apparatus. We also show that our DL model can capture the structural features of cell proliferation and recapitulate the DNA duplication through the cell cycle. The label-free structural features identified in proliferating cells are not obvious by visual inspection of the raw images, demonstrating that our holistic DL model can potentially capture previously unidentified cellular attributes with high accuracy. Another distinct advantage of multiplexing several fluorescence markers is that it enables the development of multiple quantitative metrics for imaging cytometry and phenotyping across the large cell population and at the single-cell level [e.g., cell size, nuclear-cytoplasmic ratio (NCR), nuclear roughness, and Golgi eccentricity]. As other “omic” platform technologies are rapidly being developed to evaluate biological processes at various levels (e.g., genomics, transcriptomics, proteomics, and metabolomics), single cell–level structural metrics will complement these population-based studies, particularly when a contextual phenotypic shift is expected only for a subset of cells. As a demonstration, we evaluate several cellular features, including morphology and fluorescence-expressing intensity on the DL-predicted digitally multiplexed readouts.

A common criticism of DL-based methods is the “black-box” nature of these models (23). To overcome this issue, we adapt the attention mechanism (24) to elucidate on the working mechanism of our DL model. We construct the saliency map that highlights the most important subcellular features contributing to each IF prediction by the network. Our results show that the structural components in label-free reflectance input that correspond to the fluorescence labels can be correctly identified by the saliency map, and the “attention” is consistent across different cell batches when predicting all IF labels. This indicates that our network learns to extract the salient and specific structural information from the reflectance images matching the underlying subcellular components. In addition, the improved prediction accuracy is attributed to the enhanced resolution and sensitivity to subcellular structures from the backscattering information.

RESULTS

Oblique illumination-based reflectance microscopy captures rich morphological information

The imaging platform is based on our recently developed light-emitting diode (LED) array reflectance microscope for capturing co-registered label-free reflectance and fluorescence images (8). By flexibly controlling the LED patterns, this new platform enables capturing multiple angle-dependent backscattering contrasts in the dark field without any mechanical switching. On the basis of our previous work (8), we heuristically optimize the illumination strategy and implement half-annulus LED patterns along four different orientations (including top, bottom, left, and right) (see Fig. 1A). In addition, we compute the dark-field reflectance differential phase contrast (drDPC) based on the raw measurements (see Materials and Methods). The raw oblique illumination dark-field and drDPC images contain complementary structural contrasts. In particular, subcellular structures are shown with high contrast in the raw dark-field measurements, including the nuclei, nucleoli, and hyperreflective structures at the nuclear periphery. Cell membranes with sharp boundaries are highlighted in the drDPC images, with cytoplasm spreading on the substrate with thicker nuclei at the cells’ centers. The extra cell topography information exhibited in the drDPC images are found particularly useful for predicting IF labels that are sensitive to morphology of the cell boundary. These label-free images are used as the multichannel input to our DL model. On the same platform, two-channel epifluorescence images are concurrently acquired on the same sample to serve as the ground truth for training our DL models (Fig. 1A) (see Materials and Methods). The significance of this new microscopy platform is that we capture enriched label-free information by multiple contrasts in the reflectance mode. This empowers our label-free high-content cytometry technique to uncover highly sensitive and specific structural phenotypes at the single-cell level across large cell populations.

Individual fluorescence prediction achieves state-of-the-art performance

To evaluate the performance of our DL models, we take measurements on fixed HeLa cells containing, in total, six IF labels, including DNA (Hoechst), Golgi apparatus (GM130), endosome (EEA1), actin (phalloidin), proliferation [EdU (5-ethynyl-2’-deoxyuridine)], and apoptosis [terminal deoxynucleotidyl transferase–mediated deoxyuridine triphosphate nick end labeling (TUNEL)]. Specifically, five separate batches of IF staining are performed with GM130, EEA1, phalloidin, EdU, and TUNEL, each of which is costained with Hoechst (see Materials and Methods). We then train six networks for performing individual IF label predictions using paired reflectance-fluorescence image dataset (Fig. 1B). Additional details about the network implementation and the data preprocessing procedure are provided in Materials and Methods and figs. S1 and S2, respectively.

A major goal of our study is to investigate the structural contrast captured by different label-free reflectance modes and understand how they affect the DL-based fluorescence predictions. To this end, we conduct ablation studies on the drDPC input in fig. S3 and table S1. As expected, the additional drDPC channels substantially improve the actin IF label predictions by ~9% because drDPC highlights cell topography and boundaries as compared to the plain dark-field images. Notably, the additional drDPC channels also markedly improve the prediction accuracy for proliferation and apoptosis IF labels, by ~8 and ~17%, respectively. We hypothesize that this is because distinct structural features exhibited during proliferation or apoptosis can be more prominently displayed by drDPC and subsequently recognized by the DL model.

After training, we first evaluate each network’s prediction accuracy on unseen reflectance input from the same cell batch. Figure 2 shows the label-free input, the individual IF ground truth, and the prediction for all six labels. The predicted subcellular structures and cell states have excellent visual agreement with the ground truths. Characteristic morphological features are clearly recovered, including rounded nuclei, cytoplasmic endosome, spreading cell membrane (actin), and Golgi apparatus at the nuclei periphery. Selective cellular events or functions such as proliferation and apoptosis are also captured by the DL predictions. Additional examples of the prediction results are shown in fig. S4.

Fig. 2 Visualization of the results from the six IF prediction networks.

The rows show sample dark-field reflectance images from each input stack, the network’s IF prediction, and the ground-truth IF image, respectively. The columns show six IF labels covering four different subcellular features, including nuclei (DNA), endosome, actin, and Golgi apparatus, as well as two different cell states, including proliferation and apoptosis. The IF predictions have excellent visual agreement with the ground truths in all six cases.

We quantify the prediction accuracy by computing several evaluation metrics on the network’s predictions based on the underlying cytometry tasks. Specifically, we formulate the predictions of the DNA, endosome, actin, and Golgi apparatus as regression problems because they are distinct subcellular structures and thus better described in their morphologies. Accordingly, the performance is quantified by calculating the image patch-wise Pearson correlation coefficient (PCC), which measures the pixel-level similarity of morphologies between the DL prediction and the ground truth (see Materials and Methods). The distributions of the PCCs for the four IF label predictions are shown in the violin plots in Fig. 3A. Notably, all four of our regression DL models achieve higher accuracy as compared to the current state-of-the-art techniques based on transmission label-free microscopy data (10). The median PCCs on the DNA, endosome, actin, and Golgi apparatus label predictions achieve 87.25, 91.85, 92.01, and 59.82%, respectively. These results agree well with the example visualizations in Fig. 2. The cellular features in the reflectance images associated with the corresponding fluorescence label are clearly visible. Although these scattering signals are entangled with other signals in the raw label-free images, our result shows that our DL models are able to recognize and distill these salient features with high accuracy. As visualized in the violin plots and quantified by the 25 and 75% quantiles in Fig. 3A, different amounts of variations in the prediction accuracy are seen for the four labels. In general, we observe that the larger deviation of the prediction from the ground truth, as measured by the median value, is also associated with larger accuracy variations. To investigate these spatial variations, we visualize the patch-wise PCC as a spatial map in fig. S5. Low-value PCC outliers are generally observed in background regions and hence are removed by a standard algorithm (see Materials and Methods). Other than the background outliers, the PCCs are consistent across all the cell regions for all the four IF label predictions, which demonstrates the overall robustness of the DL model.

Fig. 3 Quantitative evaluation of the DL prediction.

The violin plots show the quantitative metrics for each IF label prediction. The upper and lower bounds of each gray bar represent the 25 and 75% quantiles, respectively; the center white point marked by the black dashed line denotes the median. In total, 676 testing image patches are aggregated for computing the statistics for each label. (A) Image patch-wise PCCs of the predictions for nuclei (DNA), endosome, actin, and Golgi apparatus, which evaluates the pixel-level similarity between the regression-type predictions and ground-truth subcellular features. (B) Image patch-wise AUC of the proliferation and apoptosis predictions, which assess the pixel-level detection accuracy. (C) Quantitative evaluation of the cell-level detection performance of the proliferation and apoptosis predictions.

We formulate the prediction of the proliferation and apoptosis as detection problems, which is more biologically meaningful because they are selective cell states. The proliferation labels are only present in the DNA-replicating cells. The apoptosis labels are only for cells undergoing programmed deaths. Accordingly, the performance is first quantified by the pixel-level area under the receiver operating characteristic (ROC) curve (AUC) to measure the ability of separating the positives (i.e., those expressing the fluorescence) from the negatives at each pixel (see Materials and Methods). The calculation is image batch-wise, the same way as we calculated PCC. The distributions of the AUCs over all batches for the two IF labels are shown in the violin plots in Fig. 3B. The median detection accuracies are 91.65 and 77.24% for the proliferation and apoptosis, respectively. Next, we further evaluate the cell-level detection performance to give a more direct assessment on these two label predictions. To do so, we perform single-cell segmentation on both the predicted and ground-truth IF images and then identify each prediction as one of the four possible detection outcomes, including the true positive (TP), true negative (TN), false positive (FP), and false negative (FN) from which we compute the cell-level detection metrics, including the sensitivity and specificity (see Materials and Methods). As summarized in Fig. 3C, the cell-level proliferation prediction achieves 83.55% sensitivity and 90.92% specificity; the apoptosis prediction achieves 75.91% sensitivity and 98.88% specificity. Notably, the scattering features in the proliferating cells cannot be easily perceived from the raw reflectance images, yet our DL model can capture the salient structural features with high accuracy. The balanced high sensitivity and specificity validates the reliability of our DL models for identifying these highly selective cell states/events. Overall, these individual label prediction results validate our hypothesis that the improved sensitivity and resolution in reflectance images contain rich morphological features that can be used effectively for structural phenotyping by DL.

Multiplexed prediction recovers biological accurate cellular structures

Next, we demonstrate the digital multiplexing capability by feeding the same reflectance input to each network and make six different IF predictions in parallel. By doing so, multiple subcellular structures and cell states are revealed simultaneously. In Fig. 4A, the image multiplexes the nucleus, Golgi apparatus, actin, and endosome virtual IF labels in a single wide field of view (FOV) for a large cell population. In Fig. 4B, the virtual labels for proliferating and apoptotic cells are multiplexed with the dark-field reflectance input in the same FOV as in Fig. 4A. These multiplexed predictions are performed on the cell batch under the Golgi apparatus staining condition. To further demonstrate the robustness of this digital multiplexing procedure, we show additional examples of multiplexed predictions performed on different cell batches/staining conditions in fig. S6.

Fig. 4 Multiplexed prediction on six IF labels from the same label-free input.

(A) Visualization of the full-FOV multiplexed prediction including DNA (blue), endosome (red), actin (green), and Golgi apparatus (yellow), and (B) proliferation (cyan) and apoptosis (magenta) from the same reflectance input (gray scale). (C to R) Zoomed in of DNA, Golgi apparatus, multiplexed, and proliferation predictions. White circles indicate representative cell morphology during different phases of the cell cycle, including (C to F) interphase, (G to J) metaphase, (K to N) anaphase, and (O to R) telophase.

Our results show that the DL model can correctly capture and predict characteristic subcellular features during the cell cycle. During interphase, the nuclei have a regular rounded shape with nucleoli present, and the Golgi apparatus is anchored primarily to one side of the nuclei (Fig. 4, C to E). At this stage, cells that have initiated or ongoing DNA/chromatin replications have a positive signal for proliferation (Fig. 4F). When cells enter mitosis, the chromosomes start to condense toward the centers of the cells and the nucleoli disappears. Golgi apparatus undergoes vesiculation and fragmentation, and its components are found scattered throughout the cytoplasm in the form of tiny (~50 nm) vesicles, often referred to as the “Golgi haze” (25). During metaphase, chromosomes align at the metaphase plate and the cell shape also changes markedly, bulging into a sphere (Fig. 4, G to I). Golgi haze appears rounded, with a shaded center where chromosomes are located (Fig. 4H). During anaphase, the duplicated chromosomes separate from one another and move to opposite poles of the spindle (Fig. 4, K to M). During telophase, chromosomes start to decondense and begin to take on a more interphase-like shape (Fig. 4, O to Q). In this stage, Golgi apparatus has also completed replication and reassembled into two closely spaced cell bodies, referred to as “Golgi twins” (Fig. 4P) (25). There is no DNA replication during mitosis, so the markers for proliferation (incorporation of the fluorescent nucleoside, EdU) are absent for the metaphase, anaphase, and telophase (Fig. 4, J, N, and R). As shown in Fig. 4 (C to R), these structural, subcellular, cell cycle–dependent features are accurately captured and predicted by our DL model, which validates our hypothesis that label-free reflectance imaging and DL enable structural phenotyping.

Cell profile analysis on multiplexed images allows phenotyping and quantitative cytometry

A distinct advantage of multiplexing several markers is that it enables the development of multivariant quantitative metrics for imaging cytometry and phenotyping across the large cell population and at the single-cell level. As a demonstration, we evaluate several cellular features, including cell morphology and fluorescence intensity, on the digitally multiplexed readouts. First, we generate fluorescence intensity scatter plots similar to those used in the flow cytometry, of EdU versus Hoechst for all cells from the ground truth and digitally multiplexed IF images in Fig. 5, B and C, respectively. Evaluating the scatter plots on the population level, the DL-multiplexed prediction matches well with the ground truth, both of which show the increase in EdU and doubling of Hoechst intensity in the S and G2-M phases of the cell cycle, respectively. Next, we further evaluate the results on the individual cell level by overlaying the detection outcome for proliferation of every cell onto the same scatter plot in Fig. 5D. Specifically, by comparing the predicted proliferation IF label and the ground truth, each data point is labeled as one of the four detection outcomes (TP, TN, FP, or FN) (see details in Materials and Methods). Our analysis shows that the incorrect predictions (including FP and FN) tend to cluster around the boundaries between the S phase and S1 phases, leading to confusions in the DL predictions. There are relatively fewer incorrect predictions in the G2-M phase, which are expected because of the distinctive morphological features during mitosis that can be easily captured by the DL model (e.g., metaphase, anaphase, and telophase in Fig. 4).

Fig. 5 Cell profile analysis on digital multiplexed IF staining.

(A) Illustration of the cell cycle. (B to D) Scatter plots for the whole cell-level EdU (proliferating DNA) and Hoechst (DNA) concentrations from (B) the costained ground truth, (C) the DL prediction, and (D) the detection performance for the EdU predictions quantified by TP, TN, FP, and FN, across the entire cell population under the proliferation staining condition. The total numbers of sample cells collected from the ground truth and DL predictions are 14,778 and 14,275, respectively. The numbers of TP, FP, TN, and FN in (D) are 4089, 852, 8529, and 805, respectively (in the brackets). (E to L) Comparisons of the statistics of eight single-cell profile metrics extracted from the entire cell population in the ground truth (GT) and the DL predictions (Pred), including (E) nuclear size, (F) DNA (nuclear) fluorescence intensity contrast, (G) cell (actin) size, (H) compactness of the actin, (I) NCR measured by the area ratio between the nuclei and actin, (J) endosome size, (K) eccentricity of the Golgi apparatus distribution, and (L) concentration of the Golgi apparatus. The total numbers of sample cells collected from the ground truth and DL predictions are 20,021 and 25,257 in (E) and (F), 6183 and 5151 in (G) and (H), 2246 and 1491 in (I), 15944 and 22408 in (J), 3380 and 9392 in (K), and 3380 and 3380 in (L), respectively. A.U., arbitrary units. All the single-cell profile metrics show good agreements between the predictions and the ground truths.

In Fig. 5 (E to L), we extract several biologically relevant single-cell profile metrics using the predicted IF labels and compare them with the ground truths (see Materials and Methods), as visualized in the violin plots. In particular, we show the statistics of eight different morphological and subcellular structural parameters. In Fig. 5 (E and F), we gather statistical data about the DNA label to measure the nuclear size and intensity contrast. In Fig. 5 (G and H), we evaluate the actin size (i.e., cell size) and its compactness. In Fig. 5I, we compute the NCR, an important marker for cancers, as the area ratio between the nucleus and actin. In Fig. 5J, we measure the endosome size. In Fig. 5 (K and L), we collect morphological parameters about the Golgi apparatus, including the eccentricity and concentration. Additional metrics are provided in fig. S7. For all these single-cell profile metrics, the prediction and the ground truth show excellent agreement. These results demonstrate that our DL-augmented label-free cytometry can provide comprehensive morphological quantifications with high accuracy at the single-cell level, which is the key element for phenotyping and high-content screening (26).

Saliency map reveals inner mechanism of the deep neural network

Deep neural networks have shown high expressivity for complex models but suffer from poor explainability. Many theoretical explanations for the DL model have resorted to statistical perspective while treating the overall model as a “black box.” Instead, we use the attention-based technique (24) to elucidate on the specific label-free subcellular features that contribute to the fluorescence prediction. To do so, we treat the trained DL model as a mapping function between the input and the output. We then visualize the network’s gradient with respect to the input and extract the salient features (i.e., those having the largest gradients) the network pays most attention to (see details in Materials and Methods). The resulting “saliency map” highlights the most important features contributing to the IF prediction. By doing so, the saliency map directly evaluates the specificity of the structural features extracted from the reflectance images and how they are transformed to the target fluorescence labels by our network.

Figure 6 shows the computed saliency maps for each network across different sample batches and labeling conditions. Distinct subcellular features not only are highlighted by the network’s saliency map but also have good visual correspondence to the targeted fluorescence label. By inspecting different columns, we show morphologically distinct features from different networks, indicating that different networks can learn to recognize and focus on specific features present in the label-free images. For example, the DNA saliency maps show emphasis on nuclear boundaries and some subcellular structures. The actin saliency maps show concentration over the whole cell and spreading out to the membrane boundaries. The saliency maps for Golgi apparatus generally form shapes in “partial-moon” or circular lines. By contrast, the saliency maps for proliferation and apoptosis show that the network selectively pays attention to certain features around the nuclei. In addition, the saliency maps show that our network learns to extract invariant structural features specific to the underlying fluorescence label regardless of the cell preparation and labeling processes. Across different rows, we observe consistent saliency maps under different sample batches/staining conditions for the same labeling network.

Fig. 6 Saliency maps from each network across six cell batches with different staining conditions.

The columns show the label-free input (the first dark-field reflectance channel) and the saliency maps for six different IF labels, including DNA (blue), endosome (red), actin (green), Golgi apparatus (yellow), proliferation (cyan), and apoptosis (magenta). The rows show the label-free input and the saliency maps from six cell batches under different staining conditions. The saliency maps show good consistency across different batches and highlight distinct morphological features.

DISCUSSION

In this study, we have presented a DL-augmented label-free cytometry technique that accurately predicted six fluorescence targets in parallel at the single-cell level. The accuracy has been improved up to three times in predicting subcellular structures as compared to the current state of the art. The DL model is able to accurately recognize subcellular organelles, such as Golgi apparatus reconfigure during the cycle of proliferation, as well as to distinguish subtle morphological differences between the proliferating and nonproliferating cells. These results demonstrate the data-driven model’s unique capability of holistically extracting “nonintuitive” structural features from the label-free imaging data on a large cell population. The specificity for predicting cellular features by our DL model is illuminated by the saliency maps. This analysis demonstrated the ability of the DL models in processing highly complex and entangled structural information from scattering images.

Beyond predicting IF labels, we have further demonstrated quantitative cytometry analysis based on the multiplexed digital output from our DL models. Our analysis has shown that a multitude of single-cell profile metrics can be accurately extracted from the DL predictions. The digital multiplexing enabled us to simultaneously quantify several morphological features on multiple subcellular components across a large cell population. This capability drastically improves the technique’s throughput for structural phenotyping in the application of imaging cytometry, such as high-content analysis/screening. Cell morphological features are effective phenotypes for different disease states and environmental influences. This phenomenon is well described and practiced in pathology and cell biology. Nuclear condensation, enlargement, and increased NCR are ubiquitous hallmarks of cancers (5, 6). Cell morphology is distinct for different cell types, which is often denoted in their terminology (i.e., astrocytes, macrophage, squamous, and columnar cells), and stem cells change structures along separate differentiation paths (13). It has been shown that cell morphological changes can be directly associated with changes of morphogenic gene expressions (27), and comprehensive morphological profiling can be used to detect genetic functions (28). Meanwhile, it has been shown that DL techniques can holistically capture complex structural features for classification. This has found broad applications in detecting cell types (1214), cell states (1518, 22), drug response (19), and stem cell lineage (20). By fully leveraging the label-free and high multiplexing nature of our technique, it can potentially generate dramatic impacts in imaging cytometry by offering unprecedented information content and discovering new compound morphological features necessitating multiplexed fluorescence readouts.

Label-free, DL-augmented method of cell morphology profiling is data driven and ultimately relies on the rich information content in the images. Our LED array reflectance microscopy enables multicontrast imaging (i.e., angle-dependent dark field and drDPC) by detecting the angle-dependent backscattering signals by a programmable LED array without any mechanical moving parts. The superior sensitivity in detecting subtle structures using backscattering than transmission microscopy is well documented (2). The superb sensitivity of backscattering-based method has been demonstrated in a variety of techniques, such as partial wave spectroscopy (29), confocal light absorption and scattering spectroscopic microscopy (30), confocal reflectance quantitative phase microscopy (31), and spatial-domain low-coherence quantitative phase microscopy (32). In addition, the angle-dependent measurement has been used to measure characteristic structural length scale (6) and to enable three-dimensional (3D) reconstruction of the refractive index distribution (33). Leveraging angle-dependent reflectance signal, we outperformed the state of the art for predicting multiple subcellular components. By switching to a higher numerical aperture (NA) objective lens, the prediction performance of our DL model to subcellular structures, such as those missed in our actin label predictions, may be further improved. Recently, the LED array microscopy has also been extensively explored in transmission that allows sampling the low spatial frequency components in the Fourier space (33). In addition, backscattering spectroscopic techniques further enable characterization of ultrastructural phenotypes with sensitivity down to the nanometer-length scale (34). A potential future improvement of our imaging system is to incorporate additional transmission and multispectral LED array illumination to fully exploit the angle- and wavelength-dependent scattering contrast with a single objective lens by versatile illumination engineering.

One limitation of our current work is that it is based on fixed cells that does not allow longitudinal imaging. This can be overcome by using fluorescent reporter cell lines or live cell dyes to provide the fluorescence ground truth (10) and enable dynamic observation. The additional temporal dimension may further improve the model’s sensitivity in cell phenotyping and find previously unidentified label-free features by incorporating the information about the cell dynamics (20, 22, 35). Another limitation of the DL framework we used here is that it cannot be generalized to different types of cells. Techniques based on transfer learning (36, 37) and domain adaptation (38) will be investigated in our future work to overcome this limitation.

The variations in the prediction accuracy of our DL models are currently evaluated post hoc and based on pixel-wise and cell-level metrics by comparing the DL predictions and the ground truth. For many biomedical applications, it is beneficial to understand how much error the model may make without knowing the ground truth, i.e., the confidence of the model predictions. Emerging Bayesian DL-based uncertainty quantification techniques have proved useful to provide a proxy estimate of the prediction accuracy and quantify the model confidence (39, 40), which will be adapted in our future work.

In summary, we have reported a label-free imaging cytometry technique that multiplexes six IF labels in parallel with high accuracy via DL models. We have validated the fluorescence predictions by comparing them to the ground-truth IF images. In addition, we have conducted imaging cytometry studies on several quantitative morphological metrics on subcellular structures and phenotyping of cell proliferation. Last, the specificity of the DL model is assessed by visualizing the saliency map at the single-cell level across different staining and fixation conditions. With this unique combination of new capabilities, this novel framework may find wide applications in image-based cytometry, in particular, for high-content screening and analysis.

MATERIALS AND METHODS

Cell preparation and IF staining

HeLa cells were cultured in a Dulbecco’s modified Eagle’s medium with 10% fetal bovine serum (Gibco, 10564011) and 5% penicillin-streptomycin (Gibco, 15140122). The cells were trypsinized and passaged twice a week. Two days before the staining and imaging, cells were cultured on glass-bottom petri dishes (FluoroDish FD35-100), which were first treated with 10% poly-l-lysine (Sigma-Aldrich, RNBG0769) with phosphate-buffered saline (PBS) for 10 min in an incubator. The staining and imaging were performed on the glass-bottom dishes.

We follow the standard IF staining protocols. In total, six IF stains are used to label DNA (Hoechst), actin (phalloidin, Alexa Fluor 488 Phalloidin, Invitrogen, A12379), endosome (EEA1, Santa Cruz Biotechnology, sc-137130, AFF488), Golgi apparatus (GM130, Cell Signaling Technology, 12480), proliferation (EdU, Click-iT Plus EdU Alexa Fluor 488 kit, Invitrogen), and apoptosis (TUNEL, Click-iT TUNEL Alexa Fluor 488 kit, Invitrogen). The HeLa cells were first fixed with ice-cold methanol, washed three times (10 min each) in 0.05% PBST (0.05% Triton X-100 PBS solution), and incubated for 20 min at room temperature in a blocking solution containing 0.25% Triton X-100 and 10% bovine serum albumin in PBS. Alexa Fluor 488–conjugated antibodies were diluted in the blocking solution with the recommended concentration by the manufacturers and incubated with cells to label the specific subcellular components (EEA1 for endosome and phalloidin for actin). For actin staining, cells were fixed with ice-cold acetone to preserve the structures. For Golgi staining, a secondary antibody anti-rabbit immunoglobulin (Santa Cruz Biotechnology, 4412S) was diluted in blocking solution and used to culture the cells for 1.5 hours at room temperature in the dark. To stain cell proliferation and apoptosis, we used EdU and TUNEL assays, respectively, according to the recommended protocol by Invitrogen. The apoptosis was induced by culturing the cells with 1 μM staurosporine for 24 hours. In all the above stains, cell nuclei were counterstained with 1× Hoechst 33342.

Image data acquisition

We collect the data using our custom-built multimodal reflectance microscope (8), as shown in Fig. 1A. A custom-built LED array consisting of two LED rings is used for providing controllable dark-field illumination in reflection. We use commercially available LEDs (APTF1616SEEZGQBDC, Kingbright) that can provide three independent Red, Green, and Blue (RGB) color channels (central wavelength is 460, 515, and 630 nm, respectively). All the LEDs are individually addressable using two cascaded LED drivers (TLC5955, Texas Instruments). A microcontroller (Teensy 3.2, PJRC) provides the camera trigger signal through digital input/output pins and simultaneously controls the LED illumination pattern. The LED array is mounted around the objective lens (10× 0.3 NA, UPlanFL N, Olympus, Japan) using a 3D-printed adapter. The tube lens is a commercial single-lens reflex (SLR) lens (Nikon AF DC-NIKKOR 135mm f/2D) to maximize the FOV. The microscope provides an overall ×7.5 magnification. A scientific complementary metal-oxide–semiconductor (sCMOS) camera (CS2100M-USB, Thorlabs, 1920 × 1080 pixels, 5.04-μm pixel size, 16-bit depth) is used to acquire the images. We capture four dark-field reflectance images by using half-annulus green LED patterns along different orientations, including top (ITop), bottom (IBottom), left (ILeft), and right (IRight). The exposure time is 700 ms. Two drDPC images along two orthogonal orientations are generated byIDPC1=ITopIBottomITop+IBottom,IDPC2=ILeftIRightILeft+IRight(1)

The fluorescence excitations are provided by two LED sources (M365LP1 and M470L4, Thorlabs; central wavelength, 365 and 470 nm, respectively) combined with a dichroic mirror (DMLP425R, Thorlabs). The epifluorescence illumination is introduced by a 50/50 beam splitter (CCM1-BS013, Thorlabs). The emission filters (MF460-60 and MF525-39, Thorlabs) are placed on a filter wheel (CFW6, Thorlabs) for blue and green fluorescence emissions. Two-channel fluorescence images are acquired sequentially after acquiring the reflectance images. The exposure time is 400 ms for IF imaging. Specifically, the first green channel is for one of the five IF antibodies conjugated with the green fluorophores (Alexa Fluor 488) for endosome, actin, Golgi apparatus, proliferation, and apoptosis; the second blue channel is for the costained DNA. We capture 30 image stacks for each sample batch/IF stain.

Data preprocessing procedure

The raw reflectance and fluorescence images are preprocessed before feeding into our deep neural networks for training. The preprocessing procedure consists of four steps, including flat-field cropping, image denoising, background correction, and intensity normalization. Because the fluorescence excitation illumination is not evenly distributed across the entire rectangular FOV, we first perform flat-field cropping by using only the central 1000 × 1080–pixel region for training, where the excitation is approximately uniform. Second, we perform image denoising on the measurements. We apply two denoising approaches. In the first approach, we apply an unsupervised DL-based denoising algorithm, noise2void (41), to suppress the sensor noise present in the images. To do so, each 1000 × 1080–pixel image is cropped into 256 × 256–pixel patches. Each image patch is fed to a blind-spot network to perform denoising. After denoising, the patches are then stitched back together by alpha blending. This unsupervised denoising algorithm is found to be effective in removing unstructured, signal-independent noise, including the sensor noise and isolated hot pixels, in particular, for measurements with low signal-to-noise ratios (SNRs). The whole training and inference (denoising) procedure takes ~10 hours for processing the entire dataset containing 30 images. We find that this denoising procedure is only necessary for processing the Golgi and proliferation fluorescence images, as well as for the reflectance images for the actin prediction where the sensor noise severely corrupts the images. In the second approach, when the SNRs are sufficiently high for the measurements on other cell batches, we use a computationally more efficient morphology opening operation to remove the hot pixels in the fluorescence images under the assumption that hot pixels are isolated pixels with extreme intensity values. The opening operation takes a square kernel of size 2 × 2 pixels. This hot-pixel removal procedure takes ~15 min to process the entire dataset containing 30 images. Third, we perform background correction on the fluorescence images by eliminating the potential background bias across the batches. To do so, we calculate the histogram of each fluorescence image and denote the mode value (i.e., the most frequent value) as the constant background. This background of each fluorescence image is subtracted; the negative values from the subtraction are clipped to zero. Fourth, we perform intensity normalization by normalizing the pixel values of both the input and output images to be between 0 and 1. Additional details about the data preprocessing steps are shown in fig. S2.

Our network takes 256 × 256–pixel input images. Accordingly, we split the 30 pairs of images into 256 × 256–pixel patches to generate the training and testing data. For each IF prediction, 512 training samples and 128 testing samples are randomly generated from the full-FOV image pairs. Each input stack consists of four different channels combining dark-field reflectance images from different oblique illumination patterns (Fig. 1A). In addition, we construct the two-direction drDPC images using Eq. 1 from the dark-field images (Fig. 1A), which were found particularly effective for predicting the actin, proliferation, and apoptosis labels. We provide example comparisons of the prediction results with and without the drDPC inputs for all six IF labels in fig. S3. Quantitative comparisons of the prediction results with and without the drDPC inputs for all six IF labels are provided in table S1. When making the full-FOV predictions (e.g., Fig. 4), we use the entire 1080 × 1920–pixel reflectance images because they do not suffer from the nonuniform illumination issue.

Neural network implementation

We develop a convolutional neural network (CNN) to learn the highly complex nonlinear mapping between the morphology information contained in the multichannel reflectance images and the fluorescence labels. The network structure follows the encoder-decoder “U-net” architecture and further incorporates the dense blocks and skip connections to enable high-resolution information prediction (39). The input of the preprocessed 256 × 256–pixel reflectance image stack passes through the “encoder” path consisting of four dense blocks followed by the max-pooling layers, and the bottleneck feature maps are then fed into the “decoder” path with four dense blocks followed by upsampling layers. The skip connections bridge the lower-level activation maps with higher-level activation maps and preserve the high-frequency information. More details about the network are provided in fig. S1. We use the negative PCC as the training loss (42). We train our network using ADAM optimizer with 500 epochs and 0.1 training/validation splitting. No overfitting is observed during the training.

Quantitative evaluation of network prediction

We use the PCC to evaluate the performance of the regression type of problems. Specifically, the PCC is used to quantify the prediction quality for pervasive subcellular features, including DNA, endosome, actin, and Golgi apparatus labels. It computes the statistical correlation between the predicted and ground-truth IF image patches and is able to quantify the pixel-level similarity on the fine subcellular features. The PCC between the prediction X and the ground truth Y (each image is reshaped to an N-dimensional vector) is computed asPCC=Σi=1N(XiX¯)(YiY¯)Σi=1N(XiX¯)2Σi=1N(YiY¯)2+ϵ(2)where ϵ = 10−10 is a small regularizer to prevent zero denominator, ·¯ denotes the mean, and i is the index of each vector. The value of the PCC ranges from −1 to +1, where ±1 indicates total positive or negative correlation and 0 indicates no correlation. The PCC computation is implemented by a custom code in Python. The PCCs are computed on the testing image patches, each containing a 128 × 128–pixel FOV. In total, the statistics from 676 patches from four large FOV testing images of size 908 × 908 pixels at each staining condition are aggregated and shown in the violin plots in Fig. 3B.

To evaluate the spatial variations of the prediction performance, we construct the PCC map for each label prediction. To do so, 169 consecutive image patches are obtained from each large FOV image by cropping the image with a 128 × 128–pixel sliding window and 64 × 64–pixel overlap between the neighboring patches. The patch-wise PCC map is computed and shown in fig. S5. The lower PCC patches are found to generally align with the background region where very low reliable IF readings are present. We treat these patches as outliers. Accordingly, we use a median-based outlier detection and removal algorithm (see section S1) to remove these background outliers when constructing the violin plots in Fig. 3.

We use the AUC to quantify the detection performance for identifying the selective cell events including the proliferation and apoptosis. To plot the ROC of the detection label performance, the ground-truth labels are first binarized with certain thresholds. We use the well-established Otsu’s method to reliably compute the binarization threshold values (43), which finds the optimal values that maximize the intensity variance between the signals and the background based on a histogram clustering criterion. Specifically, the threshold values are calculated on the basis of the aggregated histogram from the entire ground-truth dataset of proliferation and apoptosis, and are found to be 0.31 and 0.11, respectively. Next, each continuous-valued pixel in the predicted image is regarded as the predicted probability of expressing the IF label at this pixel. By using the binarized ground-truth images as the target, the predictor (the trained CNN) achieves different pixel-wise true-positive rate (TPR) and false-positive rate (FPR) under different detection thresholds on the predicted images. By varying the detection thresholds, the TPR and FPR as functions of the thresholds can be plotted on the ROC curve. The AUC measures the area under the ROC curve and provides an aggregated quantification of the performance across all possible detection thresholds. The AUC is computed by the built-in functions “roc_curve” and “auc” in the scikit-learn module in Python. The AUCs are computed on the testing image patches; each contains a 128 × 128–pixel FOV. In total, the statistics from 676 testing image patches from four large FOV testing images of size 908 × 908 pixels at each condition are aggregated and shown in the violin plots in Fig. 3B.

We also quantify the single cell–level detection accuracy for the proliferation and apoptosis label predictions. To do so, we develop an automatic image processing pipeline to segment and identify each prediction as TP, TN, FP, or FN in CellProfiler. Subsequently, we compute the cell-level detection metrics, including the sensitivity and specificity. The sensitivity (also known as recall) is computed from the TPR, and the specificity (also known as selectivity) is computed from the TNR asTPR=TPP=TPTP+FN,TNR=TNN=TNTN+FP(3)

The implementation details of the image processing pipeline are described in the next section and in section S2 and fig. S8.

Digital cytometry analysis

We develop a digital cytometry analysis framework for exploring the interdependencies of different fluorescence markers on the multiplexed predictions. A commonly used flow cytometry analysis is performed by displaying the scatter plot of single cell–level proliferating DNA concentration in the log scale against the DNA concentration in the linear scale. Different from flow cytometry that directly collects the integrated fluorescence intensity from each cell, our method performs imaging with subcellular resolution across a large cell population. As a result, we first perform cell segmentation and then aggregation of the fluorescence signals within each cell region to carry out the single-cell digital cytometry analysis. We perform the digital cytometry analysis to relate the cell-level Hoechst and EdU fluorescence concentration using CellProfiler (44). Because our data contain co-registered two-channel fluorescence images with costained Hoechst and EdU, we can directly compare the ground-truth cytometry scatter plot with that from our multiplexed prediction. To generate the scatter plots, we process the co-registered Hoechst and EdU images and the multiplexed predictions using CellProfiler with a standard single-cell segmentation-based pipeline and extract the paired cell-level fluorescence intensity data. A small value (10−10) is added to the fluorescence intensity of the proliferating DNA before taking the log operation to avoid the singularity at 0. Three distinct clusters representing the S, G1, and G2-M phase are shown in the digital cytometry scatter plots. Additional details are provided in section S3 and fig. S9.

We construct an image processing pipeline to quantify the single cell–level detection performance on the proliferation and apoptosis labels (denoted as the target IF in this paragraph). First, for each of the cell event labels, we perform segmentation on the costained DNA channel ground-truth images to find all the nuclear localizations. Second, we segment all the nuclei in both the target IF prediction and ground-truth images. Third, each nucleus is labeled as one of the four possible detection outcomes. Specifically, a prediction is TP if both the prediction and the ground truth express the target IF. A prediction is TN if neither the prediction nor the ground truth expresses the target. A prediction is FP if the prediction expresses the target IF, while the ground truth does not. A prediction is FN if the prediction does not express the target IF while the ground truth does. To further investigate the effect of false detections (FP/FN) on the prediction digital cytometry scatter plot, we apply the above pipeline and label each prediction as TP, TN, FP, or FN in Fig. 5D. Additional details are provided in section S2.

Cell profile analysis

We use CellProfiler (44) to generate the single-cell profiles across each fluorescence image. We feed the ground truth and the predicted IF images of DNA, actin, endosome, and Golgi apparatus to CellProfiler. After initial cell segmentation, single cell–level parameters of morphology and intensity distribution are computed automatically by different measurement modules in CellProfiler, including the fluorescence marker size (area), compactness, eccentricity, fluorescence concentration, and single cell–level fluorescence variance and contrast. In addition, we compute the compound metric, NCR, based on the multiplexed DNA and actin fluorescence labels. To compute the NCR, co-registered ground-truth/prediction images containing the DNA and actin labels are processed individually in CellProfiler. The NCR is computed as the ratio between the area of segmented nuclei and actin masks. The single-cell profiles also present the same outliers from the background regions, which are eliminated by our outlier removal algorithm (see section S1) when constructing the violin plots. Additional details are provided in section S4.

Saliency map visualization

Our network performs image-to-image translation and can be treated as a function relating the input and output images. On the other hand, the norm of the network output is a scalar function. As a result, computing and visualizing the gradient over the norm is still possible. On the basis of this notion, we compute the saliency map as the gradient of the norm of the output with respect to the given input image. The practical implementation is achieved by the automatic differentiation feature in TensorFlow, as detailed in section S5. We specified the gradient modifier as the absolute values of the gradient, which shows the regions in the input that contribute most to the change in the output regardless of the sign of the change (i.e., negative or positive). We also used the guided backpropagation to propagate only the positive gradients for positive activations to achieve a smoother visualization. The saliency map was computed for each network and on different sample input (for varying sample batches/fixation conditions). The inputs are 256 × 256–pixel image stacks randomly selected from the testing groups under the six sample conditions. The computed saliency maps are normalized to have a uniform range between 0 and 1 for visualization.

SUPPLEMENTARY MATERIALS

Supplementary material for this article is available at http://advances.sciencemag.org/cgi/content/full/7/3/eabe0431/DC1

https://creativecommons.org/licenses/by-nc/4.0/

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.

REFERENCES AND NOTES

Acknowledgments: We thank A. Matlock for discussions on cell phenotyping visualization and evaluation and Boston University Shared Computing Cluster for proving the computational resources. Funding: This work was supported by the NSF (1846784), the NIH (R01NS108464), and the Boston University Hariri Institute Research Incubation Award. Author contributions: L.T. and J.Y. conceived the idea. S.F. and Y.M.K. prepared the cell sample culturing, fixation, and staining and acquired all the imaging data. S.C., Y.L., and Y.X. conducted the image processing, network training, cell phenotype and profile analysis, quantitative evaluation, and saliency map analysis. W.S. developed the LED array reflectance microscope platform. L.T., J.Y., and S.C. further discussed the results and refined the deep learning model, cell profiling, and cytometry pipeline. All authors contributed to the writing of the manuscript. Competing interests: The authors declare that they have no competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors.

Stay Connected to Science Advances

Navigate This Article