2.1.

Cell Culture

Normal human lung fibroblasts (NHLFs, Lonza, Basel, Switzerland) were cultured in Fibrolife S2 medium (LifeLine Cell Technology, Maryland, United States) with all provided supplements and used at p6-p9. MDA-MB-231 breast tumor cells, obtained from ATCC (United States), were transfected with RFP-Puro Lentiviral Control Vector (Cell Biolabs, Inc., California, United States). Cells were cultured in DMEM (ThermoFisher, Massachusetts, United States) at 37°C with 5% ${\mathrm{CO}}_2$. Tumor spheroids were created by plating NHLFs, stained with CellTracker™ Green CMFDA (Thermofisher, Massachusetts, United States), and RFP-transfected MDA-MB-231 breast tumor cells in a 1:1 ratio at 1M cells/ml in 10 ml of DMEM with 10% FBS on Corning^® ultra-low attachment culture dishes (CLS3261 from Corning, New York, United States). These spheroids were collected after 2 days, filtered through a $70\mu \mathrm{m}$ strainer to eliminate single cells and spheroids smaller than $70\mu \mathrm{m}$ in diameter, and then mixed with fibrin gel solution ($2\mathrm{U}/\mathrm{ml}$ thrombin, $3\mathrm{mg}/\mathrm{ml}$ fibrinogen, Millipore Sigma, Missouri, United States) to obtain a concentration of 2000 spheroids/ml. They were then injected into the central gel channel of microfluidic devices (AIM Biotech, United States). Quantifying the number of spheroids was done by aliquoting $50\mu \mathrm{l}$ of spheroid suspension solution into a flat-bottom 96-well plate well and counting the number of spheroids inside a $50\mu \mathrm{l}$ drop. After fibrin polymerization, the DMEM medium was added to the media channels flanking the gel region. One day after seeding spheroids and gel solution into devices, we changed the media for the device. On day 2, devices were fixed and permeabilized by Triton-X 100 (Millipore Sigma, Missouri, United States), and nuclei were stained by DAPI. As a proof of concept, fixation was done for ease of imaging many devices at the 2-day timepoint, but for live images in future research, Hoechst can be done instead of DAPI.

2.2.

Image Acquisition

Each microfluidic device has one or multiple sample regions, each accommodating several spheroids. Each tumor spheroid was imaged using a 20× objective (Olympus) on a confocal scanning microscope (FV-1000, Olympus, Japan) with a z-step of $4\mu \mathrm{m}$. For each z position, five types of images (channels) were recorded: blue (excitation wavelength 405 nm)/emission wavelength 461 nm) for nuclei, green (excitation 473 nm/emission 520 nm) for fibroblasts, red (excitation 559 nm, emission 572 nm) for tumor cells, and far red (excitation 635 nm, emission 668 nm, no dichroic mirror) for reflectance light and transmission light from the red laser (559 nm, transmission). We took images with typical width × length dimensions of 463 to $636\mu \mathrm{m}\times 382$ to $636\mu \mathrm{m}$ (pixel size $1.25\text{pixel}/\mu \mathrm{m}$), respectively. The height of each z-stack was varied depending on spheroid size, typically between 80 and $100\mu \mathrm{m}$, with a $4\mu \mathrm{m}$ step size. Cells near the device’s top or bottom were excluded due to strong confocal reflectance signals at these interfaces.

2.3.

Single-Cell Segmentation Using Nuclear Counterstaining and Reflectance Imaging

Next, we developed a FIJI image processing plugin that segmented individual cells in each z-stack multicellular image.²⁵ First, each z-stack image was split into smaller tiles ($159\mu \mathrm{m}\times 159\mu \mathrm{m}$ or $200\times 200\text{pixels}$), which had a 20% tile overlap [Fig. 1(a)]. The overlap between tiles allowed the removal of cells that were truncated by the tiling process. Each tile was saved as a multichannel TIFF image, reducing memory usage.

Fig. 1

Overview of the image processing and CNNs training and validation procedure. (a) Tumor spheroids are generated by co-culturing tumor cells and fibroblasts within microfluidic devices, followed by 3D cell image acquisition using confocal microscopy. (b) The 3D images are segmented into individual cells using three channels: DAPI, reflection, and transmission signals. These segmented cells are employed for CNN training, validation, and testing. (c) Cell types are distinguished by specific fluorescent signals: green denotes fibroblasts (FB), and red signifies tumor cells (TC).

To segment the cells, FIJI’s “marker-controlled watershed” algorithm needed to be used with DAPI and reflectance light (Fig. S1 in the Supplementary Material).²⁶ This algorithm treats the input image as a topographic surface, with higher gray values corresponding to greater “altitude,” and simulates a flooding process from seed points. The segmentation procedure comprised the following steps: (1) signal combination: we combined the signal intensities from the DAPI channel (representing the nucleus) and the reflectance channel (highlighting the cell’s cytoskeleton) to form a composite cell representation. Subsequently, Gaussian blur was applied to produce the “blurred cells” image (Fig. S1 in the Supplementary Material, labeled “blurred cells”). (2) Input preparation: based on the blurred image, the original DAPI images, and reflectance images, three input images required for the marker-controlled watershed plugin were generated. “Seed markers” were determined as local maxima from the DAPI image, defining nucleus locations. The “flooding input image,” reflecting cell borders, was derived from a gradient operation applied to the “blurred cells” image. (3) Mask generation: a “Mask” image was generated using the Weka 3D segmentation tool, an integrated machine learning plugin in FIJI. This Mask image identified the location of the cells’ cytoskeleton while distinguishing them from the surrounding fibrin matrix. A classifier was manually trained to differentiate cells from the matrix.²⁶ (4) The “marker-controlled watershed” module processed the input images, including “seed markers,” “flooding input image,” and “Mask,” to compute the segmented image, defining individual cells. The reflection image is crucial for defining the area of a cell, as the cell segmentation cannot be done using only the area of the nuclei defined by the DAPI signal (Fig. S1 in the Supplementary Material). Only cells that did not intersect with the image borders were retained after excluding border effects. Additionally, any background signal resulting from potential reflection at the fibrin matrix-glass slide interface, erroneously identified as segmented objects by the “marker-controlled watershed” plugin, was excluded by removing all objects touching the image border, as it consistently covered parts of the image border. The 20% overlap in the tiling process enabled the removal of cells touching the tile’s borders. Subsequently, cells were filtered based on volume, excluding those too small (less than 15 pixels or $37.9\mu {\mathrm{m}}^3$) or too large (more than 60,000 pixels or $\mathrm{151,686}\mu {\mathrm{m}}^3$). Nuclei with integrated intensities lower than the background level were also excluded. 3D regions of interest (ROIs) encompassing segmented cells were calculated using the 3D manager plugin and recorded for intensity measurement and subsequent 3D reconstruction.²⁷ The reflectance, brightfield, and DAPI images of each unlabeled cell were then used as inputs to be classified by the machine learning program [Fig. 1(b)]. Each large image features one or two spheroids. Each tile is a small part of the large image and represents a group of cells within the spheroid or cells that migrate from the spheroid to the nearby matrix. The FIJI macro initially generated over 4852 single-cell 3D images total from the original multicellular images, averaging about 121 single-cell images from each of the 40 multicellular images. Every single cell 3D image had a width and length of $159\mu \mathrm{m}\times 159\mu \mathrm{m}$ and the height of the original image (80 to $100\mu \mathrm{m}$).

Next, we used DAPI images together with a green or red signal image from each tile to create the ground-truth images of fibroblasts or tumor cells, respectively [Fig. 1(c)]. This process involved the marker-controlled watershed segmentation approach (Fig. S2 in the Supplementary Material). DAPI signal local maxima were still employed to define seed points. For the flooding input image, direct use of the green or red channel image, without additional DAPI input, was preferred due to the channel specificity of the labeled cells compared to the non-specific reflectance image. The mask image [Figs. S2(a) and S2(b) in the Supplementary Material] was obtained by automatically thresholding the fibroblast or tumor cell images using Huang’s algorithm for fibroblasts and Li’s algorithm for tumor cells in FIJI.^28,29 Subsequently, segmented fibroblasts and tumor cells were obtained within the reference channels. The cells were obtained through marker-controlled watershed segmentation, using the DAPI local maxima images as markers, the gradient of the red and green images as flooding input images, and the thresholded binary images of the red or green channels as masks. These segmented cells served as ground-truth data for NN training. In some cases, tumor cells and fibroblasts formed dense clusters, making it challenging to define cell boundaries precisely and leading to slight overlaps in ground-truth labeling. We established specific criteria to classify cells as fibroblasts or tumor cells for the training dataset selection. A cell was classified as a fibroblast when its ROI exhibited a fibroblast signal above a defined threshold. Similarly, cells with an ROI with a greater tumor cell signal than the threshold were designated as tumor cells. We introduced an intensity index we designated the $F$ index, calculated as $F/(F+C)$, where $F$ and $C$ represented the binary fibroblast and tumor cell areas within the segmented ROI. An $F$ index exceeding 0.5 denoted a fibroblast, whereas a value equal to or below 0.5 indicated a tumor cell. This criterion ensured that fibroblasts contained more green pixels than red and vice versa.

However, there were instances where cells were inaccurately segmented, resulting in erroneous cells with fragments from multiple cells. For example, the comparison of original and segmented cells highlighted in yellow in Fig. 2(a) revealed the complexity of distinguishing fibroblasts from tumor cells, especially in densely packed regions. This complexity arose because each cell in the training dataset needed a clear and exclusive ground-truth label, precluding simultaneous categorization as both a tumor cell and a fibroblast. To address this, we implemented a strategy that assigned segmented cells to ground-truth labels based on both their $F$ index (as previously defined) and their overlap with the most similar ground-truth cell. First, we confirmed a minimum 50% overlap between the cell’s ROI with either the fibroblast or cancer cell staining. Second, we ensured an overlap of over 90% with the most similar ground truth cell.

Fig. 2

Overview of machine learning training, classification and reconstruction of unseen images (a) CNN architectures and parameters are fine-tuned for optimization, ultimately leading to the selection of the best-performing neural network, K1. (b) The optimized CNN is applied to segment cells in previously unseen images. (c) Image reconstruction process. CNN-classified cells are color-coded according to the cell type and placed at the original XYZ locations of the 3D image to create the recapitulation of the original stained image.

To evaluate the segmentation accuracy, we compared the number of cells within each 3D image of $159\mu \mathrm{m}\times 159\mu \mathrm{m}\times 104\mu \mathrm{m}$ containing several cells, obtained by cropping from three large images having dimensions of $636\mu \mathrm{m}\times 636\mu \mathrm{m}\times 104\mu \mathrm{m}$ and manually or automatically counting them.

2.4.

Neural Network Training Dataset Preparation

After segmentation, the 2D slices were converted to matrices through the Numpy library.³⁰ As the cells needed to be represented as 3D images, 2D matrices of signal intensity gray value of each z-plane were stacked together as 3D matrices. We standardized individual cell image size to $20\times 50\times 50\text{pixels}$ by adding rows and columns of 0s to the cell images with dimensions smaller than the standard size. Cells larger than these dimensions were removed as they were primarily misshapen, distorted cells, likely errors in the segmentation algorithm. Finally, the 2D matrices were stacked together as 3D matrices representing 3D images of individual cells.

For use in a rigorous test of the model’s generalizing ability, 408 cells from the original 4852 single-cell 3D images were set aside before image augmentation. These cells were from six randomly chosen, original, z-stack images of spheroids not used in training, validation, or prior testing, making them completely separate from the training process. The model’s accuracy was then assessed from its performance on these 408 cells set aside (Table S1 in the Supplementary Material). We next prepared the images for the training process. The remaining 4444 3D single-cell images were first curated by removing all images where the $F$ index and overlap method disagreed, as those cells had an ambiguous ground truth. Only 277 out of 4444 segmented cells (constituting 6.2% discordance and 93.8% accordance) had uncertain ground-truth labels, often due to weak fluorescence in either the red or green channel. Consequently, we successfully assigned ground-truth labels to 4167 single cells. Next, to ensure the data is not biased toward either cell type, we removed the tumor cells with the lowest confidence until there were equal numbers of fibroblasts and tumor cells. Therefore, we randomly selected 1144 tumor cells from a pool of 3023 and included all 1144 fibroblasts. The 2288 images were copied three times and rotated 90 deg, 180 deg, and 270 deg as a form of image augmentation to increase the amount of data available to use fourfold. The resulting 9152 3D images were then arbitrarily divided into training, validation, and testing cells. The training cells are the cells that the model “learns” from, optimizing the weights it uses to classify images into cell types. The validation cells are used as a metric to determine how well the model can generalize on the cells it has not seen already during training. The strict testing cells evaluate the model’s final performance on cells it did not see during training. Each 3D matrix representing a 3D cell image from each of the three data image types (Reflection, transmission, DAPI) was stacked together into a 4D matrix to feed into the machine learning program.

2.5.

Creating the CNN Models

We developed our models and trained and evaluated their performances using the Keras and TensorFlow libraries.⁷ We hypothesized that the order of the randomized image sets presented to Keras and TensorFlow could have some effect on the machine learning’s performance. Therefore, we designed a bootstrapping program that automated the training under the same parameters, each with a different set of images for their training, validation, and testing datasets (Table S1 in the Supplementary Material). For each training run, we put each image in the training, validation, and testing datasets into a different random order before starting the multiple training epochs for that dataset. We performed 25 training runs, selecting the model with the best validation loss.

The CNN model was built on a modified version of the VGG-16 architecture.¹⁸ The CNN model was a traditional 3D CNN with a batch normalization layer, an ReLU activation function layer, an Adam optimizer, and a Max Pooling 3D layer [Fig. 1(d)]. The training run for this machine learning model created our best-performing NN, which we dubbed K1. We used the “Adam” optimizer, with an initial learning rate of 0.0001, an exponential scheduled learning rate decay, a batch size of 16, a kernel size of 5, a pooling size of (2,2,2), and 16 CNN filters.

The loss function was set to the categorical cross-entropy cost function

To evaluate the performance of our classification method, we compared the model’s prediction to the ground truth cell type for each segmented cell. Then, we recapitulate a 3D image of the original multicellular image using machine-learning classification and the 3D manager plugin by assigning the classification from the NN to each ROI and placing them in the 3D coordination within the original image to recapitulate the 3D original image.

3.1.

Image Processing of Reflectance and Counterstaining Signals to Achieve Single-Cell Segmentation in Confocal 3D Images for Convolutional Neural Network Training, Validation, and Testing

The experimental workflow involved co-culturing tumor cells and fibroblasts in a low-adhesion well plate, facilitating the formation of heterogeneous spheroids. These spheroids were subsequently transferred into microfluidic chips and cultured within fibrin gel for an additional day. During this time, they began to disperse before being fixed and imaged by confocal microscopy [Fig. 1(a)]. The confocal microscopy images comprised three channels: blue DAPI for nuclear staining, reflectance, and transmission light for non-specific extracellular and intracellular matrix structures, and brightfield images capturing cell morphology. Additionally, two fluorescence channels, green and red, were utilized to distinguish fibroblast and tumor cell labels.

In order to streamline subsequent image processing tasks and reduce computational load, the original images were divided into smaller tiles with a 20% overlap. Single-cell segmentation was performed using a marker-control segmentation algorithm in ImageJ, primarily employing the DAPI and reflectance signals [as illustrated in Fig. 1(b)]. Following the identification of ROIs corresponding to individual cells, the combination of DAPI, reflectance, and transmission signals within each ROI constituted the non-specific cell image used for training the CNN. The label assigned to each cell was determined based on the relative intensity of the green and red signals.

Out of the entire dataset of 4852 cells, a subset of 408 cells (comprising 8.4% of the total dataset), within 8 multicellular images were reserved for strict testing purposes. The remaining cells were curated and subjected to rotational augmentations to construct the training, validation, and testing datasets, facilitating the optimization of the CNN model. Detailed cell counts for each group are provided in Table S1 in the Supplementary Material. Following the optimization process, our most proficient neural network model, denoted as K1, was established [as shown in Fig. 2(a)]. Detailed CNN optimization results are described in Sec. 2.2; see below.

Subsequently, K1 was employed to classify unseen cell images from the strict testing set, as illustrated in Fig. 2(b). To generate stained images based on the CNN-classified cells, we combined cell coordinates and the classification results from the CNN. This reconstruction process resulted in the recapitulation of the original stained image [Fig. 2(c)].

To assess the segmentation performance, we compare the segmented cells (in yellow) to the green and red reference channels of the same image, which display the positions of fibroblasts and tumor cells [Fig. 3(a)]. To quantify the performance of the segmentation FIJI plugin, we counted the total number of cells manually in 12 randomly selected tile images of $159\mu \mathrm{m}\times 159\mu \mathrm{m}\times 104\mu \mathrm{m}$ and compared them with the number of cells segmented in a blinded manner. We plot the correlation graphs between the manually counted number of cells versus the one obtained with the automatic segmentation protocol [Fig. 3(b)].

Fig. 3

Single-cell segmentation and ground-truth assignment. (a) Visual comparison of 3D automatic segmentation of single cells to reference cells in a cluster of cells. Cells that touch the border of the images are not considered. (b) Comparison of automatic counting with blind manual counting across various sampled 3D images. The manual counting was conducted prior to the automatic counting in a one-side blinded setup. (c) Comparison of the original staining, the ground-truth segmentation image, and the cell type assignment of the segmented cells, excluding cells that touch the border, based on the $F$ index (see Sec. 2). Segmented cells are slightly smaller than the original stained cells.

As a result, we obtain a strong correlation between the number of cells obtained by automatic and manual counting in each individual tile image ($R^2=0.95$). We found that the ratio between the number of automatically counted and manually counted cells was 0.84. We verified the segmentation accuracy by comparing the segmented cells and the labeled images and observed that the automatic segmentation method did not pick up any false-positive cases. The program lost some cells, however, especially when they were tightly aggregated due to quasi-overlapped DAPI signals of different nuclei, causing heteroskedasticity. After segmentation, we assign the ground-truth label to each segmented cell based on the intensity of the green and red channels of the image [Fig. 3(c)]. Besides the marker-controlled watershed operation, several other 3D Segmentation algorithms in Fiji, such as 3D spot segmentation and 3D watershed were tested but visually, the segmentation accuracy of these methods are less than marker-controlled watershed because the latter uses the information from the reflectance images, which contain both signals from matrix and cells that are needed for the detection of cell borders.

3.2.

Optimized NNs Achieve 67% Classification Accuracy in the Testing Dataset

The average training time per epoch was 330 s, and the median epoch of minimum validation loss was 9 (Table 1).

Table 1

Optimized NN results. The modified VGG-16 model was trained in several training runs, and the best model was selected from these runs.

We tested our best NN (K1) against a set of 3D cell images from our 1632 cells rotated from 408 cells from the strict validation data set, which was omitted from the NN training. The NN achieves a maximum accuracy of 70% and a median accuracy of 67% (Table 1).

3.3.

Reflectance Image is Essential for the Training of CNNs

To identify what data were most important to the machine learning’s predictions, we trained the machine learning on datasets where the training data were restricted and found the optimized machine learning’s corresponding classification accuracy. We compared the machine learning’s performance when the training dataset used has either all three channels (DAPI, reflectance, and transmission signal), a combination of two channels, or only a single channel. We also restricted the number of cells the machine learning was allowed to train on.

When we average the classification accuracy of all optimized NNs obtained with each training dataset, we can confirm that the accuracy of the NNs trained with three channel-dataset is among the highest, together with DAPI and reflectance combination ($\sim 70\%$ median), or at least equal to the one obtained with either one or two channel(s) [Fig. 4(a)]. Although the DAPI-reflectance dataset led to slightly higher accuracy (70%) compared to the three-channel dataset (67%), we still use the three-channel dataset for validation and testing because the brightfield image of cells contains biological information that we might be able to extract in the future.

Fig. 4

Input features and data set tests. (a) Comparison of validation accuracy between different combinations of channels DAPI, reflectance, and brightfield as model inputs on the test set. From left to right: all three image types DAPI, reflectance (R) and brightfield (BF), DAPI and R, BF and DAPI, R and BF, only BF, only DAPI, and only R. ($n=25$ to 27). (b) Impact of the training set size on classification accuracy ($n=27$, 26, and 16, respectively, one-way ANOVA with Tukey with multiple comparisons test; *, $P<0.05$, **, $P<0.01$).

To test whether increased data would significantly improve our machine learning’s accuracy, we removed the image augmentation of rotation and cut the number of cells used in training by 90% and 50%. As expected, the accuracy was much lower with smaller data sizes, stemming from insufficient examples causing overfitting and reduced generalizability. The loss of accuracy when the data is reduced by 50% is small, implying additional data will have a minor improvement in accuracy [Fig. 4(b)].

3.4.

Recapitulation of 3D Image Succeeded in Reproducing the Original Cell Position and Type

By keeping track of which cells we used from our rigorous test dataset, we recapitulated what the predictions would have looked like as an image for each 3D ROI. We then compared the nuclei, reflectance, and transmission original images, the fibroblast and tumor cell channels of the original images, and the corresponding segmented ground-truth image with these images (Fig. 5). The recapitulated 3D images demonstrate that our method combining cell segmentation and machine learning can reflect the cell distribution of the stained 3D tissue culture images.

Fig. 5

Recapitulation of the initial 3D images using machine learning cell classification. Each row represents one representative sample. From left to right: 1. Original non-labeled images are the projection of three-channel z-stack images: nuclei stained by DAPI, reflectance, and transmission images of the confocal microscope. 2. The original labeled image of cells. 3. Ground-truth image obtained by identifying cell type (either fibroblast or tumor cell) based on the intensity of the green and red channels of the original-labeled image. ROIs without a nucleus are excluded. 4. The recapitulated 3D images by machine learning.